Update large country notebook

docsrc/notebooks/get_started/get_started.org

@@ -104,8 +104,8 @@ color_map = ListedColormap(colors)
 # 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)]
 #+end_src
 
 #+RESULTS:
@@ -160,8 +160,8 @@ raster_image = fcc_gfc.plot(ax=ax, cmap=color_map, add_colorbar=False)
 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)
 #+end_src
@@ -201,7 +201,7 @@ color_map = ListedColormap(colors)
 # 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)]
 #+end_src
 

docsrc/notebooks/large_country/large_country.org

@@ -1,44 +1,169 @@
 #+title: Large countries
-#+options: toc:nil title:t num:nil author:nil
+#+options: toc:nil title:t num:nil author:nil ^:{}
+#+property: header-args:python :results output :session :exports both
+#+property: header-args :eval never-export
+#+export_select_tags: export
+#+export_exclude_tags: noexport
 
 We can use =geefcc= to download forest cover change for large countries,
 for example Perou (iso code "PER"). The country will be divided into
-several tiles which are processed in parallel. If your computer has 8
-cores, 7 cores will be used in parallel.
+several tiles which are processed in parallel. If your computer has n
+cores, n-1 cores will be used in parallel.
 
 #+begin_src python
 import os
 import time
 
 import ee
+import numpy as np
 from geefcc import get_fcc
+import xarray as xr
+import matplotlib.pyplot as plt
+from matplotlib.colors import ListedColormap
+import matplotlib.patches as mpatches
+import cartopy.crs as ccrs
 #+end_src
 
 #+RESULTS:
-: None
+
+We initialize Google Earth Engine.
 
 #+begin_src python
 # Initialize GEE
-ee.Initialize(project="forestatrisk", opt_url="https://earthengine-highvolume.googleapis.com")
+ee.Initialize(project="forestatrisk",
+              opt_url="https://earthengine-highvolume.googleapis.com")
+#+end_src
+
+#+RESULTS:
+
+We can compute the number of cores used for the computation.
+
+#+begin_src python :results value
+ncpu = os.cpu_count() - 1
+ncpu
 #+end_src
 
+#+RESULTS:
+: 3
+
+We download the forest cover change data from GEE for Peru for years 2000, 2010 and 2020, using a buffer of about 10 km around the border (0.089... decimal degrees) and a tile size of 1 degree. 
+
 #+begin_src python
 start_time = time.time()
 get_fcc(
-    aoi="MTQ",
+    aoi="PER",
     buff=0.08983152841195216,
     years=[2000, 2010, 2020],
     source="tmf",
     tile_size=1,
-    output_file="out_tmf/fcc_tmf.tif",
+    output_file="out_tmf/forest_tmf.tif",
 )
 end_time = time.time()
 #+end_src
 
+#+RESULTS:
+
+We estimate the computation time to download 159 1-degree tiles on 3 cores. 
 
 #+begin_src python
-# Computation time
 elapsed_time = (end_time - start_time) / 60
-print('Execution time:', elapsed_time, 'minutes')
+print('Execution time:', round(elapsed_time, 2), 'minutes')
+#+end_src
+
+#+RESULTS:
+: Execution time: 1.15 minutes
+
+We "load" the data. The data can be to big to fit in memory. In fact, Xarray only load the metadata.
+
+#+begin_src python :results value
+#forest_tmf = rioxarray.open_rasterio("out_tmf/forest_tmf.tif", chunks=True)
+forest_tmf = xr.open_dataset("out_tmf/forest_tmf.tif", engine="rasterio", chunks={"band": -1, "y": 1000, "x": 1000}).astype("byte")
+forest_tmf
+#+end_src
+
+#+RESULTS:
+: <xarray.Dataset> Size: 10GB
+: Dimensions:      (band: 3, y: 68620, x: 47713)
+: Coordinates:
+:   * band         (band) int64 24B 1 2 3
+:   * x            (x) float64 382kB -81.42 -81.42 -81.42 ... -68.56 -68.56 -68.56
+:   * y            (y) float64 549kB 0.05095 0.05068 0.05041 ... -18.44 -18.44
+:     spatial_ref  int64 8B ...
+: Data variables:
+:     band_data    (band, y, x) int8 10GB dask.array<chunksize=(3, 1000, 1000), meta=np.ndarray>
+
+We transform the data to have only one band describing the forest cover change with 0 for non-forest, 1 for deforestation on the period 2000--2009, 2 for deforestation on the period 2010--2019, and 3 for the remaining forest in 2020. This is done lazily by dask.
+
+#+begin_src python :results value
+fcc_tmf = forest_tmf.sum(dim="band").astype("byte").rename({"band_data": "fcc"})
+fcc_tmf
+#+end_src
+
+#+RESULTS:
+: <xarray.Dataset> Size: 3GB
+: Dimensions:      (y: 68620, x: 47713)
+: Coordinates:
+:   * x            (x) float64 382kB -81.42 -81.42 -81.42 ... -68.56 -68.56 -68.56
+:   * y            (y) float64 549kB 0.05095 0.05068 0.05041 ... -18.44 -18.44
+:     spatial_ref  int64 8B ...
+: Data variables:
+:     fcc          (y, x) int8 3GB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
+
+We define a function to compute the mode.
+
+#+begin_src python
+def find_mode(window, axis=None, **kwargs):
+    """Find the mode over all axes."""
+    uniq = np.unique(window, return_counts=True)
+    ret = uniq[0][np.argmax(uniq[1])]
+    return [[ret]]
+#+end_src
+
+#+RESULTS:
+
+We resample at lower resolution.
+
+#+begin_src python :results value
+fcc_tmf_coarse = fcc_tmf.coarsen(x=20, y=20, boundary="pad").reduce(find_mode)
 #+end_src
 
+#+RESULTS:
+
+We prepare the colors for the map.
+
+#+begin_src python
+# Colors
+cols=[(255, 165, 0, 255), (227, 26, 28, 255), (34, 139, 34, 255)]
+colors = [(1, 1, 1, 0)]  # transparent white for 0
+cmax = 255.0  # float for division
+for col in cols:
+    col_class = tuple([i / cmax for i in col])
+    colors.append(col_class)
+color_map = ListedColormap(colors)
+
+# 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)]
+#+end_src
+
+#+RESULTS:
+
+We plot the forest cover change map.
+
+#+begin_src python :results graphics file output :file fcc.png
+# Plot
+fig = plt.figure()
+ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
+raster_image = fcc_tmf_coarse["fcc"].plot(ax=ax, cmap=color_map, add_colorbar=False)
+plt.title("Forest cover change 2000-2010-2020, TMF")
+plt.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
+fig.savefig("fcc.png", bbox_inches="tight", dpi=100)
+#+end_src
+
+#+attr_rst: :width 800 :align center
+#+RESULTS:
+[[file:fcc.png]]
+
+# End

geefcc/geeic2geotiff.py

@@ -179,25 +179,27 @@ def geeic2geotiff(index, extent, ntiles, forest, proj, scale, out_dir):
 
     """
 
-    # Open dataset
-    ds = (
-        xr.open_dataset(
-            forest,
-            engine="ee",
-            crs=proj,
-            scale=scale,
-            geometry=extent,
-            chunks={"time": -1},
+    ofile = os.path.join(out_dir, f"forest_{index}.tif")
+    if not os.path.isfile(ofile):
+        # Open dataset
+        ds = (
+            xr.open_dataset(
+                forest,
+                engine="ee",
+                crs=proj,
+                scale=scale,
+                geometry=extent,
+                chunks={"time": -1},
+            )
+            .astype("b")
+            .rename({"lon": "longitude", "lat": "latitude"})
         )
-        .astype("b")
-        .rename({"lon": "longitude", "lat": "latitude"})
-    )
-    ds = ds.load()
+        ds = ds.load()
 
-    # Load and write data to geotiff
-    xarray2geotiff(ds, "forest_cover", out_dir, index)
+        # Load and write data to geotiff
+        xarray2geotiff(ds, "forest_cover", out_dir, index)
 
-    # Progress bar
-    progress_bar_async(index, ntiles)
+        # Progress bar
+        progress_bar_async(index, ntiles)
 
 # End Of File