Source code for forestatrisk.data.compute.compute_gadm

"""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
opb = os.path.basename


[docs] def compute_gadm(ifile, ofile, proj, buff, verbose=False): """Processing GADM data. Extract layers to a new GPKG file called "aoi_latlong.gpkg" (if it doesn't exist yet) 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 ifile: Path to input GADM or aoi GPKG file. :param ofile: Path to aoi GPKG output file. :param proj: Projection definition (EPSG, PROJ.4, WKT) as in GDAL/OGR. Used for reprojecting data. :param buff: Buffer in meter (m) used to extend the extent. :param verbose: Logical. Whether to print messages or not. Default to ``False``. """ # Lat/lon AOI file aoi_latlon_file = opj(opd(ofile), "aoi_latlon.gpkg") cb = gdal.TermProgress_nocb if verbose else 0 # Change layer and file names # if using file downloaded from GADM if opb(ifile)[:4] == "gadm": # Remove intermediate output file # if it exists if os.path.isfile(aoi_latlon_file): os.remove(aoi_latlon_file) # Change layer ADM_ADM_0 and file name gdal.VectorTranslate( aoi_latlon_file, ifile, format="GPKG", layers="ADM_ADM_0", layerName="aoi", callback=cb, ) # Change layer ADM_ADM_1 # NB: 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 xmin_reg = math.floor(extent_proj[0] - buff) ymin_reg = math.floor(extent_proj[1] - buff) xmax_reg = math.ceil(extent_proj[2] + buff) ymax_reg = math.ceil(extent_proj[3] + buff) extent_reg = (xmin_reg, ymin_reg, xmax_reg, ymax_reg) return extent_reg
# End