postprocessing script CMD

sample/geslotenheid/postprocessing/openheid2indictorv2.py

@@ -41,64 +41,137 @@ def classify(x, arr, class_names):
         return class_names[np.digitize(x, arr)]
 
 
-# read data. Merk op dat hier de VERRASTERDE Viewscape output wordt gelezen. In principe zou de ruwe POINT-SHAPEFILE output
-# ook gebruikt kunnen worden. Dan zou je elk kunt kunnen classificeren en vervolgens toekennen aan een CBS 500m rapportage cel
-grid_src = r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\CBS500m_grids\versieWENR\CBS_500m_grid.shp'
-grid = gp.read_file(grid_src)
-yearA = '2018'
-yearB = '2020'
-openA = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20211213_pd20181231\4tif\monitoringsbron_openheid_2018.tif')
-openB = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\d_viewscape_output\20211208_pd20201231\4tif\monitoringsbron_openheid_2020.tif')
-
-# Output specs
-out_file_name = 'verschilgrid_openheid_{0}_{1}_V01.shp'.format(yearA, yearB)
-out_dir = r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\e_indicator\vergelijking{0}{1}'.format(yearA, yearB)
-if not os.path.isdir(out_dir):
-    os.mkdir(out_dir)
-
-# calculate gemiddelde opheid in year A voor elke 500m cell: avgopen_A
-zonal_stats = rast.zonal_stats(grid, openA.read(1), stats=['mean', 'count'], affine=openA.transform, nodata=openA.nodata)
-gemA_ha = pd.Series(data=[x['mean'] for x in zonal_stats], index=grid.index)
-countA = pd.Series(data=[x['count'] for x in zonal_stats], index=grid.index)
-print('Zonal stats {} done'.format(yearA))
-
-# calculate gemiddelde opheid in year B voor elke 500m cell: calculate avgopen_B
-zonal_stats = rast.zonal_stats(grid, openB.read(1), stats=['mean', 'count'], affine=openB.transform, nodata=openB.nodata)
-gemB_ha = pd.Series(data=[x['mean'] for x in zonal_stats], index=grid.index)
-countB = pd.Series(data=[x['count'] for x in zonal_stats], index=grid.index)
-print('Zonal stats {} done'.format(yearB))
-
-# Calculate trend between year A and year B
-diff_perc = np.subtract(np.multiply(gemB_ha.divide(gemA_ha), 100), 100)
-
-# Absolute difference
-diff_ha = gemB_ha.subtract(gemA_ha)
-
-# Classify difference in hectare to classes
-arr_vals = np.array([-100, -50, -25, 25, 50, 100])
-classes = ['1. Veel geslotener (< -100 ha)',
-           '2. Aanzienlijk geslotener (-50 - -100 ha',
-           '3. Beperkt geslotener (-25 - -50 ha)',
-           '4. Geen of nauwelijks verandering (-25 - 25 ha)',
-           '5. Beperkt opener (25 - 50 ha)',
-           '6. Aanzienlijk opener (50 - 100 ha)',
-           '7. Veel opener (>= 100 ha)']
-
-diff_class = diff_ha.apply(classify, arr=arr_vals, class_names=classes)
-
-# Append to grid shapefile and write to file
-grid['avg{}'.format(yearA)] = gemA_ha
-grid['avg{}'.format(yearB)] = gemB_ha
-grid['count{}'.format(yearA)] = countA
-grid['count{}'.format(yearB)] = countB
-grid['diff_perc'] = diff_perc
-grid['diff_ha'] = diff_ha
-grid['indicator'] = diff_class  # copy van bovenstaande, maar nu 
-grid['src_{}'.format(yearA)] = os.path.basename(openA.name)
-grid['src_{}'.format(yearB)] = os.path.basename(openB.name)
-grid['periode'] = yearB
-grid['GRID_ID'] = grid.C28992R500  # To create compliant column with cell ID
-
-# Write to file. Drop provincie column, this is added by ESRI later
-grid.drop(labels=['provincien'], axis=1).set_crs('epsg:28992', allow_override=True).to_file(os.path.join(out_dir, out_file_name))
+def to_indicator(
+    year_a: int,
+    year_b: int,
+    year_a_src: str,
+    year_b_src: str,
+    grid_src: str,
+    out_dir: str,
+    sample: bool = False,
+):
+    """
+    Take two openheid rasters and assess difference. Write to 500m grid.
+    Assumes year_b is LATER than year_a. Eg: year_b: 2020, year_a: 2018
+    """
 
+    # read data. Merk op dat hier de VERRASTERDE Viewscape output wordt gelezen. In principe zou de ruwe POINT-SHAPEFILE output
+    # ook gebruikt kunnen worden. Dan zou je elk kunt kunnen classificeren en vervolgens toekennen aan een CBS 500m rapportage cel
+
+    grid = gp.read_file(grid_src)
+    if sample:
+        grid = grid.sample(25)
+
+    yearA = str(year_a)
+    yearB = str(year_b)
+    openA = rio.open(year_a_src)
+    openB = rio.open(year_b_src)
+
+    # Output specs
+    out_file_name = "verschilgrid_openheid_{0}_{1}_V01.shp".format(yearA, yearB)
+    out_dir = os.path.join(out_dir, f"vergelijking{yearA}{yearB}")
+    if not os.path.isdir(out_dir):
+        os.mkdir(out_dir)
+
+    # calculate gemiddelde opheid in year A voor elke 500m cell: avgopen_A
+    zonal_stats = rast.zonal_stats(
+        grid,
+        openA.read(1),
+        stats=["mean", "count"],
+        affine=openA.transform,
+        nodata=openA.nodata,
+    )
+    gemA_ha = pd.Series(data=[x["mean"] for x in zonal_stats], index=grid.index)
+    countA = pd.Series(data=[x["count"] for x in zonal_stats], index=grid.index)
+    print("Zonal stats {} done".format(yearA))
+
+    # calculate gemiddelde opheid in year B voor elke 500m cell: calculate avgopen_B
+    zonal_stats = rast.zonal_stats(
+        grid,
+        openB.read(1),
+        stats=["mean", "count"],
+        affine=openB.transform,
+        nodata=openB.nodata,
+    )
+    gemB_ha = pd.Series(data=[x["mean"] for x in zonal_stats], index=grid.index)
+    countB = pd.Series(data=[x["count"] for x in zonal_stats], index=grid.index)
+    print("Zonal stats {} done".format(yearB))
+
+    # Calculate trend between year A and year B
+    diff_perc = np.subtract(np.multiply(gemB_ha.divide(gemA_ha), 100), 100)
+
+    # Absolute difference
+    diff_ha = gemB_ha.subtract(gemA_ha)
+
+    # Classify difference in hectare to classes
+    bins = np.array([-100, -50, -25, 25, 50, 100])
+    classes = [
+        "1. Veel geslotener (< -100 ha)",
+        "2. Aanzienlijk geslotener (-50 - -100 ha",
+        "3. Beperkt geslotener (-25 - -50 ha)",
+        "4. Geen of nauwelijks verandering (-25 - 25 ha)",
+        "5. Beperkt opener (25 - 50 ha)",
+        "6. Aanzienlijk opener (50 - 100 ha)",
+        "7. Veel opener (>= 100 ha)",
+    ]
+
+    diff_class = diff_ha.apply(classify, arr=bins, class_names=classes)
+
+    # Append to grid shapefile and write to file
+    grid["avg{}".format(yearA)] = gemA_ha
+    grid["avg{}".format(yearB)] = gemB_ha
+    grid["count{}".format(yearA)] = countA
+    grid["count{}".format(yearB)] = countB
+    grid["diff_perc"] = diff_perc
+    grid["diff_ha"] = diff_ha
+    grid["indicator"] = diff_class  # copy van bovenstaande, maar nu
+    grid["src_{}".format(yearA)] = os.path.basename(openA.name)
+    grid["src_{}".format(yearB)] = os.path.basename(openB.name)
+    grid["periode"] = year_b
+    grid["GRID_ID"] = grid.C28992R500  # To create compliant column with cell ID
+
+    # Write to file. Drop provincie column, this is added by ESRI later
+    grid.drop(labels=["provincien"], axis=1).set_crs(
+        "epsg:28992", allow_override=True
+    ).to_file(os.path.join(out_dir, out_file_name))
+
+
+if __name__ == "__main__":
+
+    def create_indicator(**kwargs):
+        to_indicator(
+            year_a=kwargs.get("year_a"),
+            year_a_src=kwargs.get("src_a"),
+            year_b=kwargs.get("year_b"),
+            year_b_src=kwargs.get("src_b"),
+            grid_src=kwargs.get("grid_src"),
+            out_dir=kwargs.get("out_dir"),
+            sample=kwargs.get("sample"),
+        )
+
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("year_a", help="earliest year", type=int)
+    parser.add_argument("src_a", help="path to openheid raster in year a", type=str)
+    parser.add_argument("year_b", help="latest year", type=int)
+    parser.add_argument("src_b", help="path to openheid raster in year b", type=str)
+    parser.add_argument(
+        "--grid_src",
+        help="path to grid shapefile",
+        type=str,
+        default=r"w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\CBS500m_grids\versieWENR\CBS_500m_grid.shp",
+    )
+    parser.add_argument(
+        "--out_dir",
+        help="output directory",
+        type=str,
+        default=r"w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\e_indicator",
+    )
+    parser.add_argument(
+        "--sample", help="use sample for testing purposes", action="store_true"
+    )
+    parser.set_defaults(func=create_indicator)
+    args = parser.parse_args()
+    args.func(**vars(args))
+    print("done")