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()

.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:

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))])