Adding compute_gadm.py

docsrc/python_api/forestatrisk.data.rst

@@ -14,6 +14,7 @@ Compute
 .. autofunction:: forestatrisk.data.compute.compute_biomass_avitabile
 .. autofunction:: forestatrisk.data.compute.compute_biomass_whrc
 .. autofunction:: forestatrisk.data.compute.compute_distance
+.. autofunction:: forestatrisk.data.compute.compute_gadm	
 .. autofunction:: forestatrisk.data.compute.compute_osm
 .. autofunction:: forestatrisk.data.compute.compute_srtm
 .. autofunction:: forestatrisk.data.compute.compute_wdpa

forestatrisk/data/compute/__init__.py

@@ -3,6 +3,7 @@
 from .compute_biomass_avitabile import compute_biomass_avitabile
 from .compute_biomass_whrc import compute_biomass_whrc
 from .compute_distance import compute_distance
+from .compute_gadm import compute_gadm
 from .compute_osm import compute_osm
 from .compute_srtm import compute_srtm
 from .compute_wdpa import compute_wdpa

forestatrisk/data/compute/compute_gadm.py

@@ -0,0 +1,78 @@
+"""Processing GADM data."""
+
+import os
+import math
+
+from osgeo import gdal
+
+from ..get_vector_extent import get_vector_extent
+
+opj = os.path.join
+opd = os.path.dirname
+
+
+def compute_gadm(ifile, ofile, proj, verbose=False):
+    """Processing GADM data.
+
+    Extract layers to a new GPKG file called "aoi_latlong.gpkg" and
+    reproject this file to output file. Layer names for country border
+    and subjurisdictions are "aoi" (for "area of interest") and "subj"
+    respectively. The extent of the area of interest is returned with
+    a buffer of 5 km.
+
+    :param input_file: Path to input GADM GPKG file.
+    :param output_file: Path to GPKG output file.
+    :param proj: Projection definition (EPSG, PROJ.4, WKT) as in
+        GDAL/OGR. Used for reprojecting data.
+    :param verbose: Logical. Whether to print messages or not. Default
+        to ``False``.
+
+    """
+
+    aoi_latlon_file = opj(opd(ofile), "aoi_latlon.gpkg")
+
+    cb = gdal.TermProgress if verbose else 0
+
+    gdal.VectorTranslate(
+        aoi_latlon_file, ifile,
+        format="GPKG",
+        layers="ADM_ADM_0",
+        layerName="aoi",
+        callback=cb,
+    )
+
+    # Use update access mode to not overwrite the file.
+    gdal.VectorTranslate(
+        aoi_latlon_file, ifile,
+        format="GPKG",
+        layers="ADM_ADM_1",
+        layerName="subj",
+        accessMode="update",
+        callback=cb,
+    )
+
+    # Reproject AOI
+    param = gdal.VectorTranslateOptions(
+        accessMode="overwrite",
+        srcSRS="EPSG:4326",
+        dstSRS=proj,
+        reproject=True,
+        format="GPKG",
+        callback=cb,
+    )
+    gdal.VectorTranslate(ofile, aoi_latlon_file,
+                         options=param)
+
+    # Compute extent
+    extent_proj = get_vector_extent(ofile)
+
+    # Region with buffer of 5km
+    xmin_reg = math.floor(extent_proj[0] - 5000)
+    ymin_reg = math.floor(extent_proj[1] - 5000)
+    xmax_reg = math.ceil(extent_proj[2] + 5000)
+    ymax_reg = math.ceil(extent_proj[3] + 5000)
+    extent_reg = (xmin_reg, ymin_reg, xmax_reg, ymax_reg)
+
+    return extent_reg
+
+# End

forestatrisk/data/country_compute.py

@@ -3,15 +3,14 @@
 import os
 from glob import glob
 from shutil import rmtree, copy2
-import math
-
-from osgeo import gdal
 
 from ..misc import make_dir
 from .compute import compute_forest
-from .compute import compute_srtm, compute_wdpa, compute_osm
+from .compute import (
+    compute_gadm, compute_srtm,
+    compute_wdpa, compute_osm
+)
 from .compute import compute_biomass_avitabile
-from .get_vector_extent import get_vector_extent
 
 
 def country_compute(
@@ -55,28 +54,10 @@ def country_compute(
     # Create output directory
     make_dir(output_dir)
 
-    # Reproject GADM
-    ofile = os.path.join(temp_dir, "ctry_PROJ.gpkg")
+    # Compute aoi file and extent from GADM file
     ifile = os.path.join(temp_dir, "gadm41_" + iso3 + "_0.gpkg")
-    param = gdal.VectorTranslateOptions(
-        accessMode="overwrite",
-        srcSRS="EPSG:4326",
-        dstSRS=proj,
-        reproject=True,
-        format="GPKG",
-    )
-    gdal.VectorTranslate(ofile, ifile, options=param)
-
-    # Compute extent
-    ifile = os.path.join(temp_dir, "ctry_PROJ.gpkg")
-    extent_proj = get_vector_extent(ifile)
-
-    # Region with buffer of 5km
-    xmin_reg = math.floor(extent_proj[0] - 5000)
-    ymin_reg = math.floor(extent_proj[1] - 5000)
-    xmax_reg = math.ceil(extent_proj[2] + 5000)
-    ymax_reg = math.ceil(extent_proj[3] + 5000)
-    extent_reg = (xmin_reg, ymin_reg, xmax_reg, ymax_reg)
+    ofile = os.path.join(temp_dir, "aoi_proj.gpkg")
+    extent_reg = compute_gadm(ifile, ofile, proj)
 
     # Computing country data
     if data_country:
@@ -90,7 +71,7 @@ def country_compute(
         compute_biomass_avitabile(proj, extent_reg)
         # Moving created files
         dist_files = [f for f in glob("dist_*.tif") if f[-6] != "t"]
-        proj_files = [f for f in glob("*_PROJ.*") if f != "ctry_PROJ.tif"]
+        proj_files = [f for f in glob("*_proj.*") if f != "aoi_proj.tif"]
         other_files = ["altitude.tif", "slope.tif", "pa.tif", "AGB.tif"]
         ifiles = dist_files + proj_files + other_files
         for ifile in ifiles: