class DatasetCollection
# TODO: merge ReprojectDataset and ProjectRaster they are almost the same
# TODO: still needs to be tested
@staticmethod
def to_epsg(
src: gdal.DatasetCollection,
to_epsg: int = 3857,
cell_size: int = [],
method: str = "Nearest",
) -> gdal.DatasetCollection:
"""to_epsg.
- to_epsg reprojects and resamples a folder of rasters to any projection
(default the WGS84 web mercator projection, without resampling)
Parameters
----------
src: [gdal dataset]
gdal dataset object (src=gdal.Open("dem.tif"))
to_epsg: [integer]
reference number to the new projection (https://epsg.io/)
(default 3857 the reference no of WGS84 web mercator )
cell_size: [integer]
number to resample the raster cell size to a new cell size
(default empty so raster will not be resampled)
method: [String]
resampling technique default is "Nearest"
https://gisgeography.com/raster-resampling/
"Nearest" for nearest neighbor,"cubic" for cubic convolution,
"bilinear" for bilinear
Returns
-------
raster: [gdal DatasetCollection]
a GDAL in-memory file object, where you can ReadAsArray etc.
"""
if not isinstance(src, gdal.DatasetCollection):
raise TypeError(
"src should be read using gdal (gdal dataset please read it using gdal"
f" library) given {type(src)}"
)
if not isinstance(to_epsg, int):
raise TypeError(
"please enter correct integer number for to_epsg more information "
f"https://epsg.io/, given {type(to_epsg)}"
)
if not isinstance(method, str):
raise TypeError(
"please enter correct method more information see " "docmentation "
)
if cell_size:
assert isinstance(cell_size, int) or isinstance(
cell_size, float
), "please enter an integer or float cell size"
if method == "Nearest":
method = gdal.GRA_NearestNeighbour
elif method == "cubic":
method = gdal.GRA_Cubic
elif method == "bilinear":
method = gdal.GRA_Bilinear
src_proj = src.GetProjection()
src_gt = src.GetGeoTransform()
src_x = src.RasterXSize
src_y = src.RasterYSize
dtype = src.GetRasterBand(1).DataType
spatial ref
src_sr = osr.SpatialReference(wkt=src_proj)
src_epsg = src_sr.GetAttrValue("AUTHORITY", 1)
distination
spatial ref
dst_epsg = osr.SpatialReference()
dst_epsg.ImportFromEPSG(to_epsg)
transformation factors
tx = osr.CoordinateTransformation(src_sr, dst_epsg)
incase the source crs is GCS and longitude is in the west hemisphere gdal
reads longitude fron 0 to 360 and transformation factor wont work with valeus
greater than 180
if src_epsg == "4326" and src_gt[0] > 180:
lng_new = src_gt[0] - 360
transform the right upper corner point
(ulx, uly, ulz) = tx.TransformPoint(lng_new, src_gt[3])
transform the right lower corner point
(lrx, lry, lrz) = tx.TransformPoint(
lng_new + src_gt[1] * src_x, src_gt[3] + src_gt[5] * src_y
)
else:
transform the right upper corner point
(ulx, uly, ulz) = tx.TransformPoint(src_gt[0], src_gt[3])
transform the right lower corner point
(lrx, lry, lrz) = tx.TransformPoint(
src_gt[0] + src_gt[1] * src_x, src_gt[3] + src_gt[5] * src_y
)
if not cell_size:
the result raster has the same pixcel size as the source
check if the coordinate system is GCS convert the distance from angular to metric
if src_epsg == "4326":
coords_1 = (src_gt[3], src_gt[0])
coords_2 = (src_gt[3], src_gt[0] + src_gt[1])
pixel_spacing=geopy.distance.vincenty(coords_1, coords_2).m
pixel_spacing = FeatureCollection.GCSDistance(coords_1, coords_2)
else:
pixel_spacing = src_gt[1]
else:
if src_epsg.GetAttrValue('AUTHORITY', 1) != "4326":
assert (cell_size > 1), "please enter cell size greater than 1"
if the user input a cell size resample the raster
pixel_spacing = cell_size
create a new raster
cols = int(np.round(abs(lrx - ulx) / pixel_spacing))
rows = int(np.round(abs(uly - lry) / pixel_spacing))
dst = Dataset._create_dataset(cols, rows, 1, dtype, driver="MEM")
new geotransform
new_geo = (ulx, pixel_spacing, src_gt[2], uly, src_gt[4], -pixel_spacing)
set the geotransform
dst.SetGeoTransform(new_geo)
set the projection
dst.SetProjection(dst_epsg.ExportToWkt())
set the no data value
no_data_value = src.GetRasterBand(1).GetNoDataValue()
dst = Dataset._set_no_data_value(dst, no_data_value)
perform the projection & resampling
gdal.ReprojectImage(
src, dst, src_sr.ExportToWkt(), dst_epsg.ExportToWkt(), method
)
return dst