fix #55, actually buffer neighbours by 2km, default to shapefile outp…

…ut, tweaks to remove py.test warnings

becmodel/config.py

@@ -9,7 +9,7 @@
     "becmaster": None,
     "dem": None,
     "temp_folder": tempfile.mkdtemp(prefix="becmodel-"),
-    "out_file": "becmodel.gpkg",
+    "out_file": "becmodel.shp",
     "out_layer": "becmodel",
     "cell_size_metres": 50,
     "cell_connectivity": 1,

becmodel/main.py

@@ -43,8 +43,7 @@ class ConfigValueError(Exception):
 
 
 class BECModel(object):
-    """A class to hold a model's config, data and methods
-    """
+    """A class to hold a model's config, data and methods"""
 
     def __init__(self, config_file=None):
         LOG.info("Initializing BEC model v{}".format(becmodel.__version__))
@@ -66,8 +65,7 @@ def __init__(self, config_file=None):
         self.start_time = datetime.now()
 
     def read_config(self, config_file):
-        """Read provided config file, overwriting default config values
-        """
+        """Read provided config file, overwriting default config values"""
         LOG.info("Loading config from file: %s", config_file)
         cfg = configparser.ConfigParser()
         cfg.read(config_file)
@@ -93,8 +91,7 @@ def read_config(self, config_file):
         self.config["wksp"] = self.config["temp_folder"]
 
     def update_config(self, update_dict, reload=False):
-        """Update config dictionary, reloading source data if specified
-        """
+        """Update config dictionary, reloading source data if specified"""
         self.config.update(update_dict)
         # set config temp_folder to wksp for brevity
         if "temp_folder" in update_dict.keys():
@@ -104,8 +101,7 @@ def update_config(self, update_dict, reload=False):
             self.data = util.load_tables(self.config)
 
     def validate_config(self):
-        """Validate provided config and add aspect temp zone definitions
-        """
+        """Validate provided config and add aspect temp zone definitions"""
         # validate that required paths exist
         for key in ["rulepolys_file", "elevation"]:
             if not os.path.exists(self.config[key]):
@@ -179,6 +175,16 @@ def validate_config(self):
                         self.config["dem"]
                     )
                 )
+        if Path(self.config["out_file"]).suffix == ".gpkg":
+            self.config["output_driver"] = "GPKG"
+        elif Path(self.config["out_file"]).suffix == ".shp":
+            self.config["output_driver"] = "ESRI Shapefile"
+        else:
+            raise ConfigValueError(
+                "out_file {} specified in config invalid, output must be .shp or .gpkg".format(
+                    self.config["out_file"]
+                )
+            )
 
     def write_config_log(self):
         """dump configs to file"""
@@ -347,8 +353,7 @@ def high_elevation_merges(self):
 
     @property
     def high_elevation_types(self):
-        """Create a list of high elevation types found in the entire project
-        """
+        """Create a list of high elevation types found in the entire project"""
         return list(set([k["type"] for k in self.high_elevation_merges]))
 
     @property
@@ -408,8 +413,7 @@ def high_elevation_dissolves(self):
         return high_elevation_dissolves
 
     def load(self, overwrite=False):
-        """ Load input data, do all model calculations and filters
-        """
+        """Load input data, do all model calculations and filters"""
         # shortcuts
         config = self.config
         data = self.data
@@ -472,18 +476,18 @@ def load(self, overwrite=False):
 
         # do bounds extend outside of BC?
         bounds_ll = transform_bounds("EPSG:3005", "EPSG:4326", *data["bounds"])
-        bounds_gdf = util.bbox2gdf(bounds_ll)
+        bounds_gdf = util.bbox2gdf(bounds_ll).set_crs("EPSG:4326")
 
         # load neighbours
         # Note that the natural earth dataset is only 1:10m,
-        # buffer it by 2k to be sure it captures the edge of the province
+        # buffer it by 2km to be sure it captures the edge of the province
         nbr = (
             gpd.read_file(
                 os.path.join(os.path.dirname(__file__), "data/neighbours.geojson")
             )
             .dissolve(by="scalerank")
-            .buffer(2000)
         )
+        nbr = nbr.to_crs("EPSG:3005").buffer(2000).to_crs("EPSG:4326")
         neighbours = (
             gpd.GeoDataFrame(nbr)
             .rename(columns={0: "geometry"})
@@ -545,7 +549,9 @@ def load(self, overwrite=False):
                     terraincache_path = os.environ["TERRAINCACHE"]
                 else:
                     terraincache_path = os.path.join(config["wksp"], "terrain-tiles")
-                LOG.info("Study area bounding box extends outside of BC, using MapZen terrain tiles to fill gaps")
+                LOG.info(
+                    "Study area bounding box extends outside of BC, using MapZen terrain tiles to fill gaps"
+                )
                 tt = TerrainTiles(
                     data["bounds"],
                     11,
@@ -721,7 +727,7 @@ def model(self):
         self.data = data
 
     def postfilter(self):
-        """ Tidy the output bec zones by applying several filters:
+        """Tidy the output bec zones by applying several filters:
         - majority
         - noise
         - area closing (fill in 0 areas created by noise filter)
@@ -766,13 +772,13 @@ def postfilter(self):
             majority(
                 data["becinit_grouped"],
                 morphology.rectangle(
-                    width=self.filtersize_low, height=self.filtersize_low
+                    nrows=self.filtersize_low, ncols=self.filtersize_low
                 ),
             ),
             majority(
                 data["becinit_grouped"],
                 morphology.rectangle(
-                    width=self.filtersize_steep, height=self.filtersize_steep
+                    nrows=self.filtersize_steep, ncols=self.filtersize_steep
                 ),
             ),
         )
@@ -931,38 +937,35 @@ def postfilter(self):
         )
 
         # set crs
-        data["becvalue_polys"].crs = {"init": "epsg:3005"}
+        data["becvalue_polys"].crs = "EPSG:3005"
 
         # clip to aggregated rule polygons
         # (buffer the dissolved rules out and in to ensure no small holes
         # are created by dissolve due to precision errors)
         data["rulepolys"]["rules"] = 1
-        X = data["rulepolys"].dissolve(by="rules").buffer(.01).buffer(-.01)
+        X = data["rulepolys"].dissolve(by="rules").buffer(0.01).buffer(-0.01)
         Y = gpd.GeoDataFrame(X).rename(columns={0: "geometry"}).set_geometry("geometry")
         data["becvalue_polys"] = gpd.overlay(
             data["becvalue_polys"], Y, how="intersection"
         )
 
         # add area_ha column
-        data["becvalue_polys"]["AREA_HECTARES"] = (
+        data["becvalue_polys"]["AREA_HA"] = (
             data["becvalue_polys"]["geometry"].area / 10000
         )
 
         # round to 1 decimal place
-        data["becvalue_polys"].AREA_HECTARES = data[
-            "becvalue_polys"
-        ].AREA_HECTARES.round(1)
+        data["becvalue_polys"].AREA_HA = data["becvalue_polys"].AREA_HA.round(1)
 
         # remove rulepoly fields
         data["becvalue_polys"] = data["becvalue_polys"][
-            ["BGC_LABEL", "AREA_HECTARES", "becvalue", "geometry"]
+            ["BGC_LABEL", "AREA_HA", "becvalue", "geometry"]
         ]
 
         self.data = data
 
     def write(self, discard_temp=False):
-        """ Write outputs to disk
-        """
+        """Write outputs to disk"""
 
         # if not specified otherwise, dump all intermediate raster data to file
         if not discard_temp:
@@ -1000,7 +1003,7 @@ def write(self, discard_temp=False):
         # define output schema
         schema = {
             "geometry": "MultiPolygon",
-            "properties": {"BGC_LABEL": "str:9", "AREA_HECTARES": "float:16"},
+            "properties": {"BGC_LABEL": "str:9", "AREA_HA": "float:16"},
         }
 
         # cast all features to multipolygon so that they match schema above
@@ -1011,11 +1014,14 @@ def write(self, discard_temp=False):
         ]
 
         # write output vectors to file
+        # Supported formats are shapefile or geopackage, indicated by file
+        # extension in config[out_file].
+        # **note that config[out_layer] is ignored if writing to shapefile**
         self.data["becvalue_polys"].to_file(
             self.config["out_file"],
             layer=self.config["out_layer"],
             schema=schema,
-            driver="GPKG",
+            driver=self.config["output_driver"],
         )
 
         # dump config settings to file

becmodel/util.py

@@ -1,6 +1,5 @@
 
 from math import trunc
-import os
 from pathlib import Path
 
 import logging
@@ -56,7 +55,7 @@ def load_tables(config):
         data["elevation"].rename(columns=elevation_column_remap, inplace=True)
         data["elevation"].astype(
             {
-                "beclabel": np.str,
+                "beclabel": str,
                 "cool_low": np.int16,
                 "cool_high": np.int16,
                 "neutral_low": np.int16,
@@ -158,7 +157,7 @@ def load_tables(config):
             LOG.info(
                 "Input data is not specified as BC Albers, attempting to reproject"
             )
-            data["rulepolys"] = data["rulepolys"].to_crs({"init": "EPSG:3005"})
+            data["rulepolys"] = data["rulepolys"].to_crs("EPSG:3005")
         data["rulepolys"].rename(columns=str.lower, inplace=True)
         data["rulepolys"].rename(columns=rules_column_remap, inplace=True)
     except:

sample_config.cfg

@@ -8,7 +8,7 @@ rulepolys_layer = rule_polys
 elevation       = elevation.csv
 becmaster       = becmaster.csv
 temp_folder     = becmodel_tempdata
-out_file        = becmodel.gpkg
+out_file        = becmodel.shp
 out_layer       = becmodel
 
 

tests/test.cfg

@@ -5,7 +5,7 @@ rulepolys_layer = rule_polys
 elevation = tests/data/elevation.csv
 temp_folder = tests/tempdata
 cell_size_metres = 50
-out_file = tests/bectest.gpkg
+out_file = tests/bectest.shp
 high_elevation_removal_threshold_alpine = BAFA,CMA,IMA
 high_elevation_removal_threshold_parkland = p,s
 high_elevation_removal_threshold_woodland = w

tests/test_main.py

@@ -171,9 +171,9 @@ def test_local_dem(tmpdir):
     BM.load()
 
 
-def test_run(tmpdir):
+def test_run_gpkg(tmpdir):
     """
-    Check that outputs are created, properly structured and consistent
+    Check that gpkg outputs are created, properly structured and consistent
     (not necessarily correct!)
     """
     BM = BECModel(TESTCONFIG)
@@ -191,14 +191,43 @@ def test_run(tmpdir):
     with fiona.open(os.path.join(tmpdir, "bectest.gpkg")) as output:
         assert list(output.schema["properties"].keys()) == [
             "BGC_LABEL",
-            "AREA_HECTARES",
+            "AREA_HA",
         ]
     # check outputs
     df = gpd.read_file(str(os.path.join(tmpdir, "bectest.gpkg")))
-    areas = df.groupby(["BGC_LABEL"])["AREA_HECTARES"].sum().round()
+    areas = df.groupby(["BGC_LABEL"])["AREA_HA"].sum().round()
     assert list(areas) == [5156.0, 553.0, 3617.0, 7550.0, 1511.0, 5049.0]
 
 
+def test_run_shp(tmpdir):
+    """
+    Check that shapefile outputs are created, properly structured and consistent
+    (not necessarily correct!)
+    """
+    BM = BECModel(TESTCONFIG)
+    BM.update_config({"temp_folder": str(tmpdir)})
+    BM.update_config({"out_file": str(os.path.join(tmpdir, "bectest.shp"))})
+    BM.update_config({"dem": "tests/data/dem_ok.tif"})
+    BM.load()
+    BM.model()
+    BM.postfilter()
+    BM.write()
+    assert os.path.exists(tmpdir.join("00_dem.tif"))
+    assert os.path.exists(tmpdir.join("02_aspect.tif"))
+    assert os.path.exists(tmpdir.join("bectest.shp"))
+    assert fiona.listlayers(os.path.join(tmpdir, "bectest.shp")) == ["bectest"]
+    with fiona.open(os.path.join(tmpdir, "bectest.shp")) as output:
+        assert list(output.schema["properties"].keys()) == [
+            "BGC_LABEL",
+            "AREA_HA",
+        ]
+    # check outputs
+    df = gpd.read_file(str(os.path.join(tmpdir, "bectest.shp")))
+    areas = df.groupby(["BGC_LABEL"])["AREA_HA"].sum().round()
+    # note output areas are not quite the same as when using gpkg above!
+    assert list(areas) == [5156.0, 553.0, 3619.0, 7550.0, 1510.0, 5049.0]
+
+
 def test_nowoodland(tmpdir):
     BM = BECModel(TESTCONFIG)
     BM.update_config({"temp_folder": str(tmpdir)})