Feat(rasters.py): Add raster intersection and resampling features to flopy.utils

.travis.yml

@@ -53,7 +53,10 @@ install:
 - export CXX="g++"
 - if [[ $TRAVIS_PYTHON_VERSION == 2.7 ]]; then pip install -r requirements27.travis.txt;
   echo requests version1; python -c "import requests; print(requests.__version__)";
-  else pip install -r requirements.travis.txt; fi
+  else pip install -r requirements.travis.txt;
+  pip install --no-binary rasterio rasterio;
+  pip install --upgrade numpy;
+  fi
 - pip install https://github.com/modflowpy/pymake/zipball/master
 - pip install --upgrade jupyter
 - pip install shapely[vectorize]

autotest/t065_test_gridintersect.py

@@ -1027,3 +1027,70 @@ def test_polygon_offset_rot_structured_grid_shapely():
     result = ix.intersect_polygon(p)
     # assert len(result) == 3.
     return result
+
+def test_rasters():
+    from flopy.utils import Raster
+    import os
+    import flopy as fp
+
+    ws = os.path.join("..", "examples", "data", "options")
+    raster_name = "dem.img"
+
+    try:
+        rio = Raster.load(os.path.join(ws, "dem", raster_name))
+    except ImportError:
+        return
+
+    ml = fp.modflow.Modflow.load("sagehen.nam", version="mfnwt",
+                                 model_ws=os.path.join(ws, 'sagehen'))
+    xoff = 214110
+    yoff = 4366620
+    ml.modelgrid.set_coord_info(xoff, yoff)
+
+    # test sampling points and polygons
+    val = rio.sample_point(xoff + 2000, yoff + 2000, band=1)
+    print(val - 2336.3965)
+    if abs(val - 2336.3965) > 1e-4:
+        raise AssertionError
+
+    x0, x1, y0, y1 = rio.bounds
+
+    x0 += 1000
+    y0 += 1000
+    x1 -= 1000
+    y1 -= 1000
+    shape = np.array([(x0, y0), (x0, y1), (x1, y1), (x1, y0), (x0, y0)])
+
+    data = rio.sample_polygon(shape, band=rio.bands[0])
+    if data.size != 267050:
+        raise AssertionError
+    if abs(np.min(data) - 1942.1735) > 1e-4:
+        raise AssertionError
+    if (np.max(data) - 2608.557) > 1e-4:
+        raise AssertionError
+
+    rio.crop(shape)
+    data = rio.get_array(band=rio.bands[0], masked=True)
+    if data.size != 267050:
+        raise AssertionError
+    if abs(np.min(data) - 1942.1735) > 1e-4:
+        raise AssertionError
+    if (np.max(data) - 2608.557) > 1e-4:
+        raise AssertionError
+
+    data = rio.resample_to_grid(ml.modelgrid.xcellcenters,
+                                ml.modelgrid.ycellcenters,
+                                band=rio.bands[0],
+                                method="nearest")
+    if data.size != 5913:
+        raise AssertionError
+    if abs(np.min(data) - 1942.1735) > 1e-4:
+        raise AssertionError
+    if abs(np.max(data) - 2605.6204) > 1e-4:
+        raise AssertionError
+
+    del rio
+
+
+if __name__ == "__main__":
+    test_rasters()

examples/data/options/dem/dem.img.aux.xml

@@ -0,0 +1,38 @@
+<PAMDataset>
+  <Metadata>
+    <MDI key="DataType">Generic</MDI>
+  </Metadata>
+  <Metadata domain="Esri">
+    <MDI key="PyramidResamplingType">NEAREST</MDI>
+  </Metadata>
+  <PAMRasterBand band="1">
+    <Description>Layer_1</Description>
+    <Histograms>
+      <HistItem>
+        <HistMin>1920.546752929688</HistMin>
+        <HistMax>2663.162109375</HistMax>
+        <BucketCount>256</BucketCount>
+        <IncludeOutOfRange>1</IncludeOutOfRange>
+        <Approximate>0</Approximate>
+        <HistCounts>134|325|397|465|560|738|795|788|834|858|950|819|791|953|973|883|933|1206|1766|1540|1483|1461|1789|1987|1950|1847|2074|2260|2240|2130|2414|2741|3205|2829|3080|3565|3518|3373|3543|4214|4294|3621|3424|3940|4128|3819|4117|4430|3958|3729|3208|3599|3915|3870|3599|3621|3763|3992|3597|3658|4233|4680|4572|4695|5549|5802|6150|29362|10180|9782|7344|6358|6630|6571|6265|5675|6457|7049|6732|6270|5502|5644|5135|4659|4283|4429|4117|3369|3205|3382|3696|3228|2760|2672|2858|2734|2611|2676|2778|2566|2301|2249|2428|2368|2216|2126|2258|2281|2088|2102|2202|2154|2052|1973|2002|2079|2012|1813|1885|1921|2004|1885|1896|1872|1964|1817|1695|1727|1715|1753|1756|1775|1730|1644|1645|1687|1798|1881|1773|1785|1747|1819|1751|1888|2018|1795|1522|1329|1196|1012|942|909|965|980|929|837|800|799|829|695|730|735|698|662|653|641|607|590|568|587|591|570|559|595|579|559|539|543|559|548|552|614|611|568|555|535|599|630|567|582|564|547|556|579|565|587|536|529|511|510|472|476|521|527|526|521|586|629|632|611|652|774|601|500|472|497|539|545|470|376|394|394|394|389|455|493|379|295|288|297|290|344|317|314|286|235|239|279|260|206|174|205|188|173|84|86|97|66|53|40|37|28|23|21|28|1</HistCounts>
+      </HistItem>
+    </Histograms>
+    <Metadata domain="IMAGE_STRUCTURE">
+      <MDI key="COMPRESSION">RLE</MDI>
+    </Metadata>
+    <Metadata>
+      <MDI key="LAYER_TYPE">athematic</MDI>
+      <MDI key="RepresentationType">ATHEMATIC</MDI>
+      <MDI key="STATISTICS_COVARIANCES">18624.56042831697</MDI>
+      <MDI key="STATISTICS_EXCLUDEDVALUES"></MDI>
+      <MDI key="STATISTICS_MAXIMUM">2663.162109375</MDI>
+      <MDI key="STATISTICS_MEAN">2169.7721716653</MDI>
+      <MDI key="STATISTICS_MEDIAN">0</MDI>
+      <MDI key="STATISTICS_MINIMUM">1920.5467529297</MDI>
+      <MDI key="STATISTICS_MODE">0</MDI>
+      <MDI key="STATISTICS_SKIPFACTORX">1</MDI>
+      <MDI key="STATISTICS_SKIPFACTORY">1</MDI>
+      <MDI key="STATISTICS_STDDEV">136.47183016402</MDI>
+    </Metadata>
+  </PAMRasterBand>
+</PAMDataset>

examples/data/options/dem/dem.img.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<metadata xml:lang="en"><Esri><CreaDate>20181204</CreaDate><CreaTime>08092900</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>FALSE</SyncOnce><DataProperties><itemProps><itemName Sync="TRUE">dem.img</itemName><itemLocation><linkage Sync="TRUE">file://C:\Users\jlarsen\Desktop\GSFLOW_class\trunk\exercises\saghen_prms\examples\sagehen\hru_params\dem_rasters\dem.img</linkage><protocol Sync="TRUE">Local Area Network</protocol></itemLocation><imsContentType Sync="TRUE">002</imsContentType><nativeExtBox><westBL Sync="TRUE">214270.000000</westBL><eastBL Sync="TRUE">221720.000000</eastBL><southBL Sync="TRUE">4366610.000000</southBL><northBL Sync="TRUE">4373510.000000</northBL><exTypeCode Sync="TRUE">1</exTypeCode></nativeExtBox></itemProps><coordRef><type Sync="TRUE">Projected</type><geogcsn Sync="TRUE">GCS_North_American_1983</geogcsn><csUnits Sync="TRUE">Linear Unit: Meter (1.000000)</csUnits><projcsn Sync="TRUE">NAD_1983_UTM_Zone_11N</projcsn><peXml Sync="TRUE">&lt;ProjectedCoordinateSystem xsi:type='typens:ProjectedCoordinateSystem' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/10.5'&gt;&lt;WKT&gt;PROJCS[&amp;quot;NAD_1983_UTM_Zone_11N&amp;quot;,GEOGCS[&amp;quot;GCS_North_American_1983&amp;quot;,DATUM[&amp;quot;D_North_American_1983&amp;quot;,SPHEROID[&amp;quot;GRS_1980&amp;quot;,6378137.0,298.257222101]],PRIMEM[&amp;quot;Greenwich&amp;quot;,0.0],UNIT[&amp;quot;Degree&amp;quot;,0.0174532925199433]],PROJECTION[&amp;quot;Transverse_Mercator&amp;quot;],PARAMETER[&amp;quot;False_Easting&amp;quot;,500000.0],PARAMETER[&amp;quot;False_Northing&amp;quot;,0.0],PARAMETER[&amp;quot;Central_Meridian&amp;quot;,-117.0],PARAMETER[&amp;quot;Scale_Factor&amp;quot;,0.9996],PARAMETER[&amp;quot;Latitude_Of_Origin&amp;quot;,0.0],UNIT[&amp;quot;Meter&amp;quot;,1.0],AUTHORITY[&amp;quot;EPSG&amp;quot;,26911]]&lt;/WKT&gt;&lt;XOrigin&gt;-5120900&lt;/XOrigin&gt;&lt;YOrigin&gt;-9998100&lt;/YOrigin&gt;&lt;XYScale&gt;450445547.3910538&lt;/XYScale&gt;&lt;ZOrigin&gt;-100000&lt;/ZOrigin&gt;&lt;ZScale&gt;10000&lt;/ZScale&gt;&lt;MOrigin&gt;-100000&lt;/MOrigin&gt;&lt;MScale&gt;10000&lt;/MScale&gt;&lt;XYTolerance&gt;0.001&lt;/XYTolerance&gt;&lt;ZTolerance&gt;0.001&lt;/ZTolerance&gt;&lt;MTolerance&gt;0.001&lt;/MTolerance&gt;&lt;HighPrecision&gt;true&lt;/HighPrecision&gt;&lt;WKID&gt;26911&lt;/WKID&gt;&lt;LatestWKID&gt;26911&lt;/LatestWKID&gt;&lt;/ProjectedCoordinateSystem&gt;</peXml></coordRef><RasterProperties><General><PixelDepth Sync="TRUE">32</PixelDepth><HasColormap Sync="TRUE">FALSE</HasColormap><CompressionType Sync="TRUE">RLE</CompressionType><NumBands Sync="TRUE">1</NumBands><Format Sync="TRUE">IMAGINE Image</Format><HasPyramids Sync="TRUE">TRUE</HasPyramids><SourceType Sync="TRUE">continuous</SourceType><PixelType Sync="TRUE">floating point</PixelType><NoDataValue Sync="TRUE">-3.4028235e+038</NoDataValue></General></RasterProperties><lineage><Process ToolSource="c:\arcgis\desktop10.5\ArcToolbox\Toolboxes\Data Management Tools.tbx\Clip" Date="20181204" Time="080929">Clip ..\examples\sagehen\dem\imgn40w121_13_filled.img "214063.834030762 4366573.51459876 221714.946388238 4373594.09248124" ..\examples\sagehen\hru_params\dem_rasters\dem_clip.img # # NONE NO_MAINTAIN_EXTENT</Process><Process ToolSource="c:\arcgis\desktop10.5\ArcToolbox\Toolboxes\Data Management Tools.tbx\ProjectRaster" Date="20181204" Time="080933">ProjectRaster ..\examples\sagehen\hru_params\dem_rasters\dem_clip.img ..\examples\sagehen\hru_params\dem_rasters\dem.img PROJCS['NAD_1983_UTM_Zone_11N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-117.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]] BILINEAR "10 10" # "214110 4366620" PROJCS['NAD_1983_UTM_Zone_11N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-117.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]</Process></lineage></DataProperties><SyncDate>20181204</SyncDate><SyncTime>08093300</SyncTime><ModDate>20181204</ModDate><ModTime>08093300</ModTime></Esri><dataIdInfo><envirDesc Sync="TRUE"> Version 6.2 (Build 9200) ; Esri ArcGIS 10.5.1.7333</envirDesc><dataLang><languageCode value="eng" Sync="TRUE"></languageCode><countryCode value="USA" Sync="TRUE"></countryCode></dataLang><idCitation><resTitle Sync="TRUE">dem.img</resTitle><presForm><PresFormCd value="005" Sync="TRUE"></PresFormCd></presForm></idCitation><spatRpType><SpatRepTypCd value="002" Sync="TRUE"></SpatRepTypCd></spatRpType><dataExt><geoEle><GeoBndBox esriExtentType="search"><exTypeCode Sync="TRUE">1</exTypeCode><westBL Sync="TRUE">-120.321168</westBL><eastBL Sync="TRUE">-120.231832</eastBL><northBL Sync="TRUE">39.466206</northBL><southBL Sync="TRUE">39.401697</southBL></GeoBndBox></geoEle></dataExt></dataIdInfo><mdLang><languageCode value="eng" Sync="TRUE"></languageCode><countryCode value="USA" Sync="TRUE"></countryCode></mdLang><mdChar><CharSetCd value="004" Sync="TRUE"></CharSetCd></mdChar><distInfo><distFormat><formatName Sync="TRUE">Raster Dataset</formatName></distFormat></distInfo><mdHrLv><ScopeCd value="005" Sync="TRUE"></ScopeCd></mdHrLv><mdHrLvName Sync="TRUE">dataset</mdHrLvName><refSysInfo><RefSystem><refSysID><identCode code="26911" Sync="TRUE"></identCode><idCodeSpace Sync="TRUE">EPSG</idCodeSpace><idVersion Sync="TRUE">6.13(3.0.1)</idVersion></refSysID></RefSystem></refSysInfo><spatRepInfo><Georect><cellGeo><CellGeoCd Sync="TRUE" value="002"></CellGeoCd></cellGeo><numDims Sync="TRUE">2</numDims><tranParaAv Sync="TRUE">1</tranParaAv><chkPtAv Sync="TRUE">0</chkPtAv><cornerPts><pos Sync="TRUE">214270.000000 4366610.000000</pos></cornerPts><cornerPts><pos Sync="TRUE">214270.000000 4373510.000000</pos></cornerPts><cornerPts><pos Sync="TRUE">221720.000000 4373510.000000</pos></cornerPts><cornerPts><pos Sync="TRUE">221720.000000 4366610.000000</pos></cornerPts><centerPt><pos Sync="TRUE">217995.000000 4370060.000000</pos></centerPt><axisDimension type="002"><dimSize Sync="TRUE">745</dimSize><dimResol><value Sync="TRUE" uom="m">10.000000</value></dimResol></axisDimension><axisDimension type="001"><dimSize Sync="TRUE">690</dimSize><dimResol><value Sync="TRUE" uom="m">10.000000</value></dimResol></axisDimension><ptInPixel><PixOrientCd Sync="TRUE" value="001"></PixOrientCd></ptInPixel></Georect></spatRepInfo><contInfo><ImgDesc><contentTyp><ContentTypCd Sync="TRUE" value="001"></ContentTypCd></contentTyp><covDim><Band><dimDescrp Sync="TRUE">Layer_1</dimDescrp><maxVal Sync="TRUE">2663.162109</maxVal><minVal Sync="TRUE">1920.546753</minVal><bitsPerVal Sync="TRUE">32</bitsPerVal><valUnit><UOM type="length"></UOM></valUnit></Band></covDim></ImgDesc></contInfo><mdDateSt Sync="TRUE">20181204</mdDateSt></metadata>

examples/data/options/dem/model_boundary.CPG

@@ -0,0 +1 @@
+UTF-8

examples/data/options/dem/model_boundary.prj

@@ -0,0 +1 @@
+PROJCS["NAD_1983_UTM_Zone_11N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-117.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

examples/data/options/dem/model_boundary.shp.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<metadata xml:lang="en"><Esri><CreaDate>20190816</CreaDate><CreaTime>10383000</CreaTime><ArcGISFormat>1.0</ArcGISFormat><SyncOnce>TRUE</SyncOnce><DataProperties><itemProps><itemLocation><linkage Sync="TRUE">file://\\IGSWCGWWLT4175\C$\Users\jlarsen\Desktop\dev-geospatial\dem\model_boundary</linkage><protocol Sync="TRUE">Local Area Network</protocol></itemLocation></itemProps><copyHistory><copy source="C:\Users\jlarsen\Desktop\dev-geospatial\model_boundary" dest="\\IGSWCGWWLT4175\C$\Users\jlarsen\Desktop\dev-geospatial\dem\model_boundary" date="20190816" time="10383000"></copy></copyHistory></DataProperties></Esri></metadata>

flopy/utils/__init__.py

@@ -44,3 +44,4 @@
 from .recarray_utils import create_empty_recarray, ra_slice
 from .mtlistfile import MtListBudget
 from .optionblock import OptionBlock
+from .rasters import Raster