Adding crop_to_aoi argument

docsrc/notebooks/get_started/get_started.org

@@ -41,12 +41,15 @@ if not os.path.isfile("out_tmf/forest_tmf.tif"):
         aoi="REU",
         years=[2000, 2010, 2020],
         source="tmf",
+        parallel=False,
+        crop_to_aoi=True,
         tile_size=0.5,
         output_file="out_tmf/forest_tmf.tif",
     )
 #+end_src
 
 #+RESULTS:
+: get_fcc running, 3 tiles .
 
 #+begin_src python :results value
 # Load data
@@ -137,12 +140,15 @@ if not os.path.isfile("out_gfc_50/forest_gfc_50.tif"):
         years=[2001, 2010, 2020],  # Here, first year must be 2001 (1st Jan)
         source="gfc",
         perc=50,
+        parallel=False,
+        crop_to_aoi=True,
         tile_size=0.5,
         output_file="out_gfc_50/forest_gfc_50.tif",
     )
 #+end_src
 
 #+RESULTS:
+: get_fcc running, 3 tiles ....
 
 #+begin_src python
 # Load data

docsrc/notebooks/get_started/get_started.rst

@@ -39,10 +39,17 @@ from the Tropical Moist Forest product. We will use the Reunion Island
             aoi="REU",
             years=[2000, 2010, 2020],
             source="tmf",
+            parallel=False,
+            crop_to_aoi=True,
             tile_size=0.5,
             output_file="out_tmf/forest_tmf.tif",
         )
 
+::
+
+    get_fcc running, 3 tiles .
+
+
 .. code:: python
 
     # Load data
@@ -98,8 +105,8 @@ from the Tropical Moist Forest product. We will use the Reunion Island
     # Labels
     labels = {0: "non-forest in 2000", 1:"deforestation 2000-2009",
               2:"deforestation 2010-2019", 3:"forest in 2020"}
-    patches =[mpatches.Patch(facecolor=col, edgecolor="black",
-                             label=labels[i]) for (i, col) in enumerate(colors)]
+    patches = [mpatches.Patch(facecolor=col, edgecolor="black",
+                              label=labels[i]) for (i, col) in enumerate(colors)]
 
 .. code:: python
 
@@ -129,10 +136,17 @@ Compare with forest cover change from GFC
             years=[2001, 2010, 2020],  # Here, first year must be 2001 (1st Jan)
             source="gfc",
             perc=50,
+            parallel=False,
+            crop_to_aoi=True,
             tile_size=0.5,
             output_file="out_gfc_50/forest_gfc_50.tif",
         )
 
+::
+
+    get_fcc running, 3 tiles ....
+
+
 .. code:: python
 
     # Load data
@@ -148,8 +162,8 @@ Compare with forest cover change from GFC
     plt.title("Forest cover change 2001-2010-2020, GFC")
     labels = {0: "non-forest in 2001", 1:"deforestation 2001-2009",
               2:"deforestation 2010-2019", 3:"forest in 2020"}
-    patches =[mpatches.Patch(facecolor=col, edgecolor="black",
-                             label=labels[i]) for (i, col) in enumerate(colors)]
+    patches = [mpatches.Patch(facecolor=col, edgecolor="black",
+                              label=labels[i]) for (i, col) in enumerate(colors)]
     plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
     fig.savefig("gfc.png", bbox_inches="tight", dpi=100)
 
@@ -185,7 +199,7 @@ Comparing forest cover in 2020 between TMF and GFC
     # Labels
     labels = {0: "non-forest tmf, non-forest gfc", 1:"non-forest tmf / forest gfc",
               2:"forest tmf / forest gfc", 3:"forest tmf, non-forest gfc"}
-    patches =[mpatches.Patch(facecolor=col, edgecolor="black",
+    patches = [mpatches.Patch(facecolor=col, edgecolor="black",
                              label=labels[i]) for (i, col) in enumerate(colors)]
 
 .. code:: python

geefcc/geotiff_from_tiles.py

@@ -2,19 +2,21 @@
 
 import os
 from glob import glob
-import math
 
 from osgeo import gdal
 
 opj = os.path.join
 opd = os.path.dirname
 
 
-def geotiff_from_tiles(output_file):
+def geotiff_from_tiles(crop_to_aoi, extent, output_file):
     """Make geotiff from tiles.
 
-    :param extent_latlong: Extent in lat/long.
-    :param scale: Resolution.
+    :param crop_to_aoi: Crop the raster using aoi extent.
+
+    :param extent: Result of ``get_extent_from_aoi()`` function (a
+        dictionary).
+
     :param output_file: Output file.
 
     """
@@ -44,9 +46,29 @@ def geotiff_from_tiles(output_file):
     # VRT to GeoTIFF
     # Creation options
     copts = ["COMPRESS=DEFLATE", "BIGTIFF=YES"]
-    gdal.Translate(output_file, vrt_file,
-                   maskBand=None,
-                   creationOptions=copts,
-                   callback=cback)
+    aoi_isfile = extent["aoi_isfile"]
+    borders_gpkg = extent["borders_gpkg"]
+    extent_latlong = extent["extent_latlong"]
+    if crop_to_aoi:
+        if aoi_isfile:
+            gdal.Warp(output_file, vrt_file,
+                      cropToCutline=True,
+                      warpOptions=["CUTLINE_ALL_TOUCHED=TRUE"],
+                      cutlineDSName=borders_gpkg,
+                      creationOptions=copts,
+                      callback=cback)
+        else:
+            xmin, ymin, xmax, ymax = extent_latlong
+            ulx_uly_lrx_lry = [xmin, ymax, xmax, ymin]
+            gdal.Translate(output_file, vrt_file,
+                           projWin=ulx_uly_lrx_lry,
+                           maskBand=None,
+                           creationOptions=copts,
+                           callback=cback)
+    else:
+        gdal.Translate(output_file, vrt_file,
+                       maskBand=None,
+                       creationOptions=copts,
+                       callback=cback)
 
 # End

geefcc/get_fcc.py

@@ -21,7 +21,7 @@ def get_fcc(aoi,
             source="tmf",
             perc=75,
             tile_size=1,
-            # crop_to_aoi=False,
+            crop_to_aoi=False,
             parallel=False,
             ncpu=None,
             output_file="fcc.tif"):
@@ -52,17 +52,17 @@ def get_fcc(aoi,
 
     :param tile_size: Tile size for parallel computing.
 
+    :param crop_to_aoi: Crop the raster GeoTIFF file to **aoi with
+        buffer**. If ``False``, the output file will match the
+        **grid** covering the aoi with buffer.
+
     :param parallel: Logical. Parallel (if ``True``) or sequential (if
         ``False``) computing. Default to ``False``.
 
     :param ncpu: Number of CPU to use for parallel computing. If None,
         it will be set to the number of cores on the computer minus
         one.
 
-    :param crop_to_aoi: Crop the raster GeoTIFF file to aoi. If False,
-        the output file will match the grid covering the aoi with
-        buffer.
-
     :param output_file: Path to output GeoTIFF file. If directories in
         path do not exist they will be created.
 
@@ -78,10 +78,10 @@ def get_fcc(aoi,
     scale = 0.000269494585235856472  # in dd, ~30 m
 
     # Get aoi
-    ext = get_extent_from_aoi(aoi, buff, out_dir)
-    aoi_isfile = ext["aoi_isfile"]
-    borders_gpkg = ext["borders_gpkg"]
-    extent_latlong = ext["extent_latlong"]
+    extent = get_extent_from_aoi(aoi, buff, out_dir)
+    aoi_isfile = extent["aoi_isfile"]
+    borders_gpkg = extent["borders_gpkg"]
+    extent_latlong = extent["extent_latlong"]
 
     # Make minimal grid
     grid_gpkg = opj(out_dir, "grid.gpkg")
@@ -136,6 +136,6 @@ def get_fcc(aoi,
             pool.join()
 
     # Geotiff from tiles
-    geotiff_from_tiles(output_file)
+    geotiff_from_tiles(crop_to_aoi, extent, output_file)
 
 # End