GDAL to 3.3.0

Dockerfile

@@ -1,40 +1,41 @@
 # Basin Setup requires GDAL
-FROM geographica/gdal2:2.3.2
+FROM osgeo/gdal:ubuntu-full-3.3.0
 
 ENV TAUDEM_HASH f927ca639a1834565a76cb3df5acbcd2909d6d0d
 
+WORKDIR /code
+
 ADD https://github.com/dtarb/TauDEM/archive/${TAUDEM_HASH}.tar.gz .
 
 RUN apt-get -y update \
-    && apt-get install -y git mpich wget zip unzip python3-pip \
+    && apt-get install -y mpich git python3-pip libnetcdf-dev \
     && apt-get autoremove
 
 # Build taudem
-RUN tar -xzf ${TAUDEM_HASH}.tar.gz -C /usr/src \
-    && mkdir /usr/src/TauDEM-${TAUDEM_HASH}/bin \
-    && cd /usr/src/TauDEM-${TAUDEM_HASH}/src \
-    && make
-
-RUN mkdir -p /usr/local/taudem/bin \
-    && cp /usr/src/TauDEM-${TAUDEM_HASH}/bin/* /usr/local/taudem/bin/
+RUN tar -xzf ${TAUDEM_HASH}.tar.gz \
+    && mkdir /code/TauDEM-${TAUDEM_HASH}/bin \
+    && cd /code/TauDEM-${TAUDEM_HASH}/src \
+    && make \
+    && mkdir -p /usr/local/taudem/bin \
+    && cp /code/TauDEM-${TAUDEM_HASH}/bin/* /usr/local/taudem/bin/ \
+    && rm -rf /code/TauDEM-${TAUDEM_HASH}
 
 ENV PATH /usr/local/taudem/bin:$PATH
 
 # Copy over code base
-RUN mkdir -p /code \
-    && mkdir -p /code/basin_setup
+# RUN mkdir -p /code \
+#     && mkdir -p /code/basin_setup
 
 COPY . /code/basin_setup/
 
 RUN cd /code/basin_setup \
-    && python3 -m pip install --upgrade pip \
-    && python3 -m pip install setuptools wheel \
-    && python3 -m pip install -r requirements.txt \
+    && python3 -m pip install --no-cache-dir setuptools wheel \
+    && python3 -m pip install --no-cache-dir -r requirements.txt \
     && python3 setup.py install \
     && rm -rf /var/lib/apt/lists/*
 
 # Setup the shared volume and add it as the starting place
 VOLUME /data
 WORKDIR /data
-RUN echo "umask 0002" >> /etc/bash.bashrc
+
 ENTRYPOINT ["generate_topo"]

basin_setup/generate_topo/main.py

@@ -105,7 +105,7 @@ def load_dem(self):
         gdal.gdalwarp(
             self.config['dem_file'],
             self.images['dem'],
-            self.crs['init'],
+            self.crs,
             self.extents,
             self.cell_size,
             resample='bilinear',
@@ -140,7 +140,7 @@ def load_vegetation(self):
             elif self.config['vegetation_dataset'] == 'landfire_2.0.0':
                 veg = vegetation.Landfire200(self.config)
 
-            veg.reproject(self.extents, self.cell_size, self.crs['init'])
+            veg.reproject(self.extents, self.cell_size, self.crs)
             veg.load_clipped_images()
             veg.calculate_tau_and_k()
             veg.calculate_height()
@@ -199,9 +199,6 @@ def create_netcdf(self):
             'standard_name': 'projection_y_coordinate'
         }
 
-        for key in list(output.keys()):
-            output[key].attrs["grid_mapping"] = "projection"
-
         # Global attributes
         now = datetime.now().strftime('%Y-%m-%d %H:%M:%S')
         output.attrs = {
@@ -222,10 +219,12 @@ def create_netcdf(self):
             encoding={
                 "x": {"dtype": "f4"},
                 "y": {"dtype": "f4"},
-                "veg_type": {"dtype": 'u2'},
-                "veg_height": {"dtype": "f4"},
-                "veg_tau": {"dtype": "f4"},
-                "veg_k": {"dtype": "f4"},
+                "dem": {"dtype": "f4", "grid_mapping": "projection"},
+                "mask": {"grid_mapping": "projection"},
+                "veg_type": {"dtype": 'u2', "grid_mapping": "projection"},
+                "veg_height": {"dtype": "f4", "grid_mapping": "projection"},
+                "veg_tau": {"dtype": "f4", "grid_mapping": "projection"},
+                "veg_k": {"dtype": "f4", "grid_mapping": "projection"},
             }
         )
 

basin_setup/generate_topo/shapefile.py

@@ -10,7 +10,12 @@ def __init__(self, file_name) -> None:
 
     @property
     def crs(self):
-        return self.polygon.crs
+        """The EPSG code for the shapefile
+
+        Returns:
+            str: EPSG code, ex 'epsg:32611'
+        """
+        return self.polygon.crs.srs
 
     def mask(self, nx, ny, transform):
         """Create a raster mask from the shapefile using rasterio.features.rasterize

requirements.txt

@@ -1,14 +1,12 @@
 colorama==0.3.9
-numpy==1.16.1
-pandas>=0.24.0,<1.0.0
-pyshp==1.2.12
-requests==2.20.0
-geopandas==0.4.0
-utm==0.4.2
+numpy<1.20
+pandas>=1.0
+requests>=2.20.0
+geopandas>=0.9.0,<0.10
 coloredlogs==10.0
 spatialnc>=0.3.0,<0.4.0
 inicheck>=0.9,<0.10.0
-netcdf4==1.4.1
+netcdf4>=1.4.1
 cftime<1.1.0
 xarray>=0.16.2
 rioxarray>=0.3.1

scripts/gold_difference

@@ -23,12 +23,12 @@ def argument_parser():
     )
 
     parser.add_argument(
-        '--basin', '-b',
-        metavar='BASIN',
+        '--veg', '-v',
+        metavar='VEGETATION',
         type=str,
-        default='Lakes',
-        choices=['Lakes'],
-        help='Test basin name (Lakes)'
+        default='landfire_140',
+        choices=['landfire_140', 'landfire_200'],
+        help='Vegetation dataset to compare'
     )
 
     parser.add_argument(
@@ -58,8 +58,9 @@ def main(args):
             os.path.join(
                 repo_path,
                 'tests',
-                args.basin,
+                'Lakes',
                 'gold',
+                args.veg,
                 "{}.nc".format(gf)
             )
         )
@@ -78,4 +79,4 @@ def main(args):
 
 if __name__ == '__main__':
     script_arguments = argument_parser().parse_args()
-    main(script_arguments)
+    main(script_arguments)

tests/Lakes/data/pour_points.bna

@@ -1,2 +1,2 @@
 "Lakes","",1
-322895.0,4165837.8
+322895.0, 4165837.8

tests/Lakes/data/pour_points.geojson

@@ -0,0 +1,19 @@
+{
+  "type": "FeatureCollection",
+  "name": "pour_points",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "Primary ID": "Lakes",
+      },
+      "geometry": {
+        "type": "Point",
+        "coordinates": [
+          322895,
+          4165837.8
+        ]
+      }
+    }
+  ]
+}

tests/Lakes/gold/basin_outline.prj

@@ -1 +1 @@
-PROJCS["WGS_1984_UTM_Zone_11N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]
+PROJCS["WGS_1984_UTM_Zone_11N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],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]]

tests/Lakes/gold/landfire_140/gold_config.ini

@@ -1,11 +1,11 @@
 [generate_topo]
 
-basin_shapefile:        ./basin_outline.shp
-dem_file:               ../data/dem_epsg_32611_100m.tif
+basin_shapefile:        ../basin_outline.shp
+dem_file:               ../../data/dem_epsg_32611_100m.tif
 output_folder:          ./
 basin_name:             Lakes
 
-vegetation_folder:      ../../../landfire/landfire_140
+vegetation_folder:      ../../../../landfire/landfire_140
 vegetation_dataset:     landfire_1.4.0
 
 pad_domain:             5, 5, 5, 5

tests/Lakes/lakes_test_case.py

@@ -15,12 +15,13 @@
 class BasinSetupLakes(unittest.TestCase):
 
     BASE_INI_FILE_NAME = 'config.ini'
-    CRS = dict({'init': 'epsg:32611'})
+    CRS = 'epsg:32611'
 
     test_dir = Path(basin_setup.__file__).parent.parent.joinpath('tests')
     basin_dir = test_dir.joinpath('Lakes')
     gold_dir = test_dir.joinpath(basin_dir, 'gold')
     config_file = os.path.join(basin_dir, BASE_INI_FILE_NAME)
+    output_topo = os.path.join(basin_dir, 'output', 'topo.nc')
 
     @classmethod
     def base_config_copy(cls):

tests/generate_topo/test_generate_topo.py

@@ -1,15 +1,12 @@
-import os
 from unittest.mock import patch
 
-import numpy as np
 import xarray as xr
 from inicheck.config import UserConfig
-from inicheck.tools import cast_all_variables
 from rasterio import Affine
 
 from basin_setup.generate_topo import GenerateTopo
 from basin_setup.generate_topo.shapefile import Shapefile
-from basin_setup.generate_topo.vegetation import Landfire140, Landfire200
+from basin_setup.generate_topo.vegetation import Landfire140
 from basin_setup.utils import domain_extent
 from tests.Lakes.lakes_test_case import BasinSetupLakes
 
@@ -55,7 +52,7 @@ def test_load_basin_shapefiles(self):
 
         self.assertTrue(len(self.subject.basin_shapefiles) == 1)
         self.assertIsInstance(self.subject.basin_shapefiles[0], Shapefile)
-        self.assertDictEqual(self.subject.basin_shapefiles[0].crs, self.CRS)
+        self.assertEqual(self.subject.basin_shapefiles[0].crs, self.CRS)
 
     def test_load_dem(self):
         self.subject.crs = self.CRS
@@ -98,7 +95,7 @@ def test_run(self, mock_veg):
         gt = GenerateTopo(config_file=self.config_file)
         gt.run()
 
-        ds = xr.open_dataset(os.path.join(self.basin_dir, 'output', 'topo.nc'))
+        ds = xr.open_dataset(self.output_topo, cache=False)
 
         self.assertCountEqual(
             list(ds.coords.keys()),
@@ -109,98 +106,6 @@ def test_run(self, mock_veg):
             ['dem', 'mask', 'veg_height', 'veg_k',
                 'veg_tau', 'veg_type', 'projection']
         )
+        ds.close()
 
         self.compare_netcdf_files('landfire_140/topo.nc', 'topo.nc')
-
-
-@patch.object(Landfire140, 'veg_height_csv',
-              new='tests/Lakes/data/landfire_1.4.0/LF_140EVH_05092014.csv')
-@patch.object(Landfire140, 'clipped_images', new={
-    'veg_type': 'tests/Lakes/data/landfire_1.4.0/clipped_veg_type.tif',
-    'veg_height': 'tests/Lakes/data/landfire_1.4.0/clipped_veg_height.tif'
-})
-@patch.object(Landfire200, 'veg_height_csv',
-              new='tests/Lakes/data/landfire_2.0.0/LF16_EVH_200.csv')
-@patch.object(Landfire200, 'clipped_images', new={
-    'veg_type': 'tests/Lakes/data/landfire_2.0.0/clipped_veg_type.tif',
-    'veg_height': 'tests/Lakes/data/landfire_2.0.0/clipped_veg_height.tif'
-})
-class TestVegetationOptions(BasinSetupLakes):
-
-    @patch.object(Landfire140, 'reproject', return_value=True)
-    def test_landfire_140(self, mock_reproject):
-        gt = GenerateTopo(config_file=self.config_file)
-        gt.run()
-        self.assertTrue(mock_reproject.called)
-        self.assertTrue(mock_reproject.call_count == 1)
-
-        ds = xr.open_dataset(os.path.join(self.basin_dir, 'output', 'topo.nc'))
-
-        self.assertCountEqual(
-            list(ds.coords.keys()),
-            ['y', 'x']
-        )
-        self.assertCountEqual(
-            list(ds.keys()),
-            ['dem', 'mask', 'veg_height', 'veg_k',
-                'veg_tau', 'veg_type', 'projection']
-        )
-
-        self.compare_netcdf_files('landfire_140/topo.nc', 'topo.nc')
-
-    @patch.object(Landfire200, 'reproject', return_value=True)
-    def test_landfire_200(self, mock_reproject):
-        config = self.base_config_copy()
-
-        config.raw_cfg['generate_topo']['vegetation_dataset'] = 'landfire_2.0.0'  # noqa
-        config.raw_cfg['generate_topo']['vegetation_folder'] = '../../landfire/landfire_200'  # noqa
-
-        config.apply_recipes()
-        config = cast_all_variables(config, config.mcfg)
-
-        gt = GenerateTopo(config_file=config)
-        gt.run()
-        self.assertTrue(mock_reproject.called)
-        self.assertTrue(mock_reproject.call_count == 1)
-
-        ds = xr.open_dataset(os.path.join(self.basin_dir, 'output', 'topo.nc'))
-
-        self.assertCountEqual(
-            list(ds.coords.keys()),
-            ['y', 'x']
-        )
-        self.assertCountEqual(
-            list(ds.keys()),
-            ['dem', 'mask', 'veg_height', 'veg_k',
-                'veg_tau', 'veg_type', 'projection']
-        )
-
-        self.compare_netcdf_files('landfire_200/topo.nc', 'topo.nc')
-
-    @patch.object(Landfire140, 'reproject', return_value=True)
-    def test_no_veg(self, mock_reproject):
-        gt = GenerateTopo(config_file=self.config_file)
-        gt.config['vegetation_dataset'] = None
-        gt.run()
-
-        self.assertFalse(mock_reproject.called)
-        self.assertTrue(mock_reproject.call_count == 0)
-
-        ds = xr.open_dataset(os.path.join(self.basin_dir, 'output', 'topo.nc'))
-
-        self.assertCountEqual(
-            list(ds.coords.keys()),
-            ['y', 'x']
-        )
-        self.assertCountEqual(
-            list(ds.keys()),
-            ['dem', 'mask', 'veg_height', 'veg_k',
-                'veg_tau', 'veg_type', 'projection']
-        )
-
-        for image in gt.veg.VEG_IMAGES:
-            if image == 'veg_type':
-                # all veg type are 0 values
-                self.assertTrue(np.sum(ds[image].isnull().values) == 0)
-            else:
-                self.assertTrue(np.all(ds[image].isnull().values))

tests/generate_topo/test_shapefile.py

@@ -21,7 +21,7 @@ def test_load(self):
         self.assertTrue(len(self.shape.polygon) == 1)
 
     def test_crs(self):
-        self.assertDictEqual(self.shape.crs, self.CRS)
+        self.assertEqual(self.shape.crs, self.CRS)
 
     def test_mask(self):
         transform, x, y = domain_extent.affine_transform_from_extents(