From bad483b3910218dc828c993863d540793111090d Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Wed, 17 Jul 2024 09:44:13 -0700 Subject: [PATCH] update(Raster): add new methods and checks (#2267) * update(Raster): add new methods and checks * add feature `to_crs()` for re-projecting rasters * add static method `raster_from_array()` to allow users to make rasters from data * update `resample_to_grid()` - removed multithread and thread_pool kwargs - added initial check for raster/modelgrid intersection * add testing for raster improvements * Updates for raster_intersection_example.py * Linting * linting part 2 * Catch point with GeoSpatialUtil() in add_region() --- .../Notebooks/raster_intersection_example.py | 27 +- autotest/test_gridintersect.py | 108 ++++++++ flopy/utils/rasters.py | 256 ++++++++++++++++-- flopy/utils/triangle.py | 1 + 4 files changed, 374 insertions(+), 18 deletions(-) diff --git a/.docs/Notebooks/raster_intersection_example.py b/.docs/Notebooks/raster_intersection_example.py index 4c5eecbcd..447581faa 100644 --- a/.docs/Notebooks/raster_intersection_example.py +++ b/.docs/Notebooks/raster_intersection_example.py @@ -1,11 +1,12 @@ # --- # jupyter: # jupytext: +# notebook_metadata_filter: metadata # text_representation: # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.14.5 +# jupytext_version: 1.14.4 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -570,6 +571,30 @@ # The `ibound` array and the `top` array can be used to build or edit the BAS and DIS file objects in FloPy +# ## Raster re-projection +# +# The `Raster` class has a built in `to_crs()` method that allows for raster reprojection. The `to_crs()` method has two possible parameters that can be used to define reprojection and one additional parameter for in place reprojection: +# +# - `crs`: the crs parameter can take many different formats of coordinate refence systems (WKT string, epsg code, pyproj.CRS, rasterio.CRS, proj4 string, epsg string, etc...) +# - `epsg`: integer epsg number +# - `inplace`: bool, default False creates a new raster object, True modifies the existing Raster object +# +# Here's example usage: + +cur_crs = rio.crs +print(cur_crs) +print(rio.transform) + +rio_reproj = rio.to_crs(crs="EPSG:4326") # WGS84 dec. lat lon +print(rio_reproj.crs) +print(rio_reproj.transform) + +# Reproject as an inplace operation + +rio.to_crs(epsg=4326, inplace=True) +print(rio.crs) +print(rio.transform) + # ## Future development # # Potential features that draw on this functionality could include: diff --git a/autotest/test_gridintersect.py b/autotest/test_gridintersect.py index 4497cc3a2..268de09f9 100644 --- a/autotest/test_gridintersect.py +++ b/autotest/test_gridintersect.py @@ -1449,6 +1449,114 @@ def test_raster_sampling_methods(example_data_path): ) +@requires_pkg("rasterio") +def test_raster_reprojection(example_data_path): + ws = example_data_path / "options" / "dem" + raster_name = "dem.img" + + wgs_epsg = 4326 + wgs_xmin = -120.32116799649168 + wgs_ymax = 39.46620605907534 + + raster = Raster.load(ws / raster_name) + + print(raster.crs.to_epsg()) + wgs_raster = raster.to_crs(crs=f"EPSG:{wgs_epsg}") + + if not wgs_raster.crs.to_epsg() == wgs_epsg: + raise AssertionError(f"Raster not converted to EPSG {wgs_epsg}") + + transform = wgs_raster._meta["transform"] + if not np.isclose(transform.c, wgs_xmin) and not np.isclose( + transform.f, wgs_ymax + ): + raise AssertionError(f"Raster not reprojected to EPSG {wgs_epsg}") + + raster.to_crs(epsg=wgs_epsg, inplace=True) + transform2 = raster._meta["transform"] + for ix, val in enumerate(transform): + if not np.isclose(val, transform2[ix]): + raise AssertionError("In place reprojection not working") + + +@requires_pkg("rasterio") +def test_create_raster_from_array_modelgrid(example_data_path): + ws = example_data_path / "options" / "dem" + raster_name = "dem.img" + + raster = Raster.load(ws / raster_name) + + xsize = 200 + ysize = 100 + xmin, xmax, ymin, ymax = raster.bounds + + nbands = 5 + nlay = 1 + nrow = int(np.floor((ymax - ymin) / ysize)) + ncol = int(np.floor((xmax - xmin) / xsize)) + + delc = np.full((nrow,), ysize) + delr = np.full((ncol,), xsize) + + grid = flopy.discretization.StructuredGrid( + delc=delc, + delr=delr, + top=np.ones((nrow, ncol)), + botm=np.zeros((nlay, nrow, ncol)), + idomain=np.ones((nlay, nrow, ncol), dtype=int), + xoff=xmin, + yoff=ymin, + crs=raster.crs, + ) + + array = np.random.random((grid.ncpl * nbands,)) * 100 + robj = Raster.raster_from_array(array, grid) + + if nbands != len(robj.bands): + raise AssertionError("Number of raster bands is incorrect") + + array = array.reshape((nbands, nrow, ncol)) + for band in robj.bands: + ra = robj.get_array(band) + np.testing.assert_allclose( + array[band - 1], + ra, + err_msg="Array not properly reshaped or converted to raster", + ) + + +@requires_pkg("rasterio", "affine") +def test_create_raster_from_array_transform(example_data_path): + import affine + + ws = example_data_path / "options" / "dem" + raster_name = "dem.img" + + raster = Raster.load(ws / raster_name) + + transform = raster._meta["transform"] + array = raster.get_array(band=raster.bands[0]) + + array = np.expand_dims(array, axis=0) + # same location but shrink raster by factor 2 + new_transform = affine.Affine( + transform.a / 2, 0, transform.c, 0, transform.e / 2, transform.f + ) + + robj = Raster.raster_from_array( + array, crs=raster.crs, transform=new_transform + ) + + rxmin, rxmax, rymin, rymax = robj.bounds + xmin, xmax, ymin, ymax = raster.bounds + + if ( + not ((xmax - xmin) / (rxmax - rxmin)) == 2 + or not ((ymax - ymin) / (rymax - rymin)) == 2 + ): + raise AssertionError("Transform based raster not working properly") + + if __name__ == "__main__": sgr = get_rect_grid(angrot=45.0, xyoffset=10.0) ls = LineString([(5, 10.0 + np.sqrt(200.0)), (15, 10.0 + np.sqrt(200.0))]) diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py index e56fee2cd..fe9155121 100644 --- a/flopy/utils/rasters.py +++ b/flopy/utils/rasters.py @@ -95,10 +95,8 @@ def __init__( if isinstance(crs, CRS): pass - elif isinstance(crs, int): - crs = CRS.from_epsg(crs) - elif isinstance(crs, str): - crs = CRS.from_string(crs) + elif crs is not None: + crs = CRS.from_user_input(crs) else: TypeError("crs type not understood, provide an epsg or proj4") @@ -126,6 +124,13 @@ def __init__( if isinstance(rio_ds, rasterio.io.DatasetReader): self._dataset = rio_ds + @property + def crs(self): + """ + Returns a rasterio CRS object + """ + return self._meta["crs"] + @property def bounds(self): """ @@ -140,6 +145,13 @@ def bounds(self): return xmin, xmax, ymin, ymax + @property + def transform(self): + """ + Returns the affine transform for the raster + """ + return self._meta["transform"] + @property def bands(self): """ @@ -184,6 +196,126 @@ def ycenters(self): self.__xycenters() return self.__ycenters + def to_crs(self, crs=None, epsg=None, inplace=False): + """ + Method for re-projecting rasters from an existing CRS to a + new CRS + + Parameters + ---------- + crs : CRS user input of many different kinds + epsg : int + epsg code input that defines the coordinate system + inplace : bool + Boolean flag to indicate if the operation takes place "in place" + which reprojects the raster within the current object or the + default (False) to_crs() returns a reprojected Raster object + + Returns + ------- + Raster or None: returns a reprojected raster object if + inplace=False, otherwise the reprojected information + overwrites the current Raster object + + """ + from rasterio.crs import CRS + + if self.crs is None: + raise ValueError( + "Cannot transform naive geometries. " + "Please set a crs on the object first." + ) + if crs is not None: + dst_crs = CRS.from_user_input(crs) + elif epsg is not None: + dst_crs = CRS.from_epsg(epsg) + else: + raise ValueError("Must pass either crs or epsg.") + + # skip if the input CRS and output CRS are the exact same + if self.crs.to_epsg() == dst_crs.to_epsg(): + return self + + return self.__transform(dst_crs=dst_crs, inplace=inplace) + + def __transform(self, dst_crs, inplace): + """ + + Parameters + ---------- + dst_crs : rasterio.CRS object + inplace : bool + + Returns + ------- + Raster or None: returns a reprojected raster object if + inplace=False, otherwise the reprojected information + overwrites the current Raster object + + """ + import rasterio + from rasterio.io import MemoryFile + from rasterio.warp import ( + Resampling, + calculate_default_transform, + reproject, + ) + + height = self._meta["height"] + width = self._meta["width"] + xmin, xmax, ymin, ymax = self.bounds + + transform, width, height = calculate_default_transform( + self.crs, dst_crs, width, height, xmin, ymin, xmax, ymax + ) + + kwargs = { + "transform": transform, + "width": width, + "height": height, + "crs": dst_crs, + "nodata": self.nodatavals[0], + "driver": self._meta["driver"], + "count": self._meta["count"], + "dtype": self._meta["dtype"], + } + + with MemoryFile() as memfile: + with memfile.open(**kwargs) as dst: + for band in self.bands: + reproject( + source=self.get_array(band), + destination=rasterio.band(dst, band), + src_transform=self.transform, + src_crs=self.crs, + dst_transform=transform, + dst_crs=dst_crs, + resampling=Resampling.nearest, + ) + with memfile.open() as dataset: + array = dataset.read() + bands = dataset.indexes + meta = dataset.meta + + if inplace: + for ix, band in enumerate(bands): + self.__arr_dict[band] = array[ix] + + self.__xcenters = None + self.__ycenters = None + self._meta.update({k: v for k, v in kwargs.items()}) + self._dataset = None + + else: + return Raster( + array, + bands, + meta["crs"], + meta["transform"], + meta["nodata"], + meta["driver"], + ) + def __xycenters(self): """ Method to create np.arrays of the xy-cell centers @@ -327,8 +459,6 @@ def resample_to_grid( modelgrid, band, method="nearest", - multithread=False, - thread_pool=2, extrapolate_edges=False, ): """ @@ -363,11 +493,6 @@ def resample_to_grid( `'mode'` for majority sampling - multithread : bool - DEPRECATED boolean flag indicating if multithreading should be - used with the ``mean`` and ``median`` sampling methods - thread_pool : int - DEPRECATED number of threads to use for mean and median sampling extrapolate_edges : bool boolean flag indicating if areas without data should be filled using the ``nearest`` interpolation method. This option @@ -377,16 +502,18 @@ def resample_to_grid( ------- np.array """ - if multithread: - warnings.warn( - "multithread option has been deprecated and will be removed " - "in flopy version 3.3.8" - ) - import_optional_dependency("scipy") rasterstats = import_optional_dependency("rasterstats") from scipy.interpolate import griddata + xmin, xmax, ymin, ymax = modelgrid.extent + rxmin, rxmax, rymin, rymax = self.bounds + if any([rxmax < xmin, rxmin > xmax, rymax < ymin, rymin > ymax]): + raise AssertionError( + "Raster and model grid do not intersect. Check that the grid " + "and raster are in the same coordinate reference system" + ) + method = method.lower() if method in ("linear", "nearest", "cubic"): xc = modelgrid.xcellcenters @@ -770,6 +897,101 @@ def load(raster: Union[str, os.PathLike]): meta["driver"], ) + @staticmethod + def raster_from_array( + array, + modelgrid=None, + nodataval=1e-10, + crs=None, + transform=None, + ): + """ + Method to create a raster from an array. When using a modelgrid to + define the transform, delc and delr must be uniform in each dimension. + Otherwise, the user can define their own transform using the affine + package. + + Parameters + ---------- + array : np.ndarray + array of (n-bands, nrows, ncols) for the raster + modelgrid : flopy.discretization.StructuredGrid + StructuredGrid object (optional), but transform must be defined + if a StructuredGrid is not supplied + nodataval : (int, float) + Null value + crs : coordinate reference system input of many types + transform : affine.Affine + optional affine transform that defines the spatial parameters + of the raster. This must be supplied if a modelgrid is not + used to define the transform + + Returns + ------- + Raster object + """ + from affine import Affine + + if not isinstance(array, np.ndarray): + array = np.array(array) + + if modelgrid is not None: + if crs is None: + if modelgrid.crs is None: + raise ValueError( + "Cannot create a raster from a grid without a " + "coordinate reference system, please provide a crs " + "using crs=" + ) + crs = modelgrid.crs + + if modelgrid.grid_type != "structured": + raise TypeError( + f"{type(modelgrid)} discretizations are not supported" + ) + + if not np.all(modelgrid.delc == modelgrid.delc[0]): + raise AssertionError("DELC must have a uniform spacing") + + if not np.all(modelgrid.delr == modelgrid.delr[0]): + raise AssertionError("DELR must have a uniform spacing") + + yul = modelgrid.yvertices[0, 0] + xul = modelgrid.xvertices[0, 0] + angrot = modelgrid.angrot + transform = Affine( + modelgrid.delr[0], 0, xul, 0, -modelgrid.delc[0], yul + ) + + if angrot != 0: + transform *= Affine.rotation(angrot) + + if array.size % modelgrid.ncpl != 0: + raise AssertionError( + f"Array size {array.size} is not a multiple of the " + f"number of cells per layer in the model grid " + f"{modelgrid.ncpl}" + ) + + array = array.reshape((-1, modelgrid.nrow, modelgrid.ncol)) + + if transform is not None: + if crs is None: + raise ValueError( + "Cannot create a raster without a coordinate reference " + "system, please use crs= to provide a coordinate reference" + ) + + bands, height, width = array.shape + + return Raster( + array, + bands=list(range(1, bands + 1)), + crs=crs, + transform=transform, + nodataval=nodataval, + ) + def plot(self, ax=None, contour=False, **kwargs): """ Method to plot raster layers or contours. diff --git a/flopy/utils/triangle.py b/flopy/utils/triangle.py index d001bd73c..e378fb318 100644 --- a/flopy/utils/triangle.py +++ b/flopy/utils/triangle.py @@ -131,6 +131,7 @@ def add_region(self, point, attribute=0, maximum_area=None): None """ + point = GeoSpatialUtil(point, shapetype="point").points self._regions.append([point, attribute, maximum_area]) def build(self, verbose=False):