Plot with imshow, not xarray

docsrc/notebooks/large_country/large_country.org

@@ -20,12 +20,10 @@ import ee
 import numpy as np
 from osgeo import gdal
 from geefcc import get_fcc, sum_raster_bands
-import xarray as xr
 import geopandas
 import matplotlib.pyplot as plt
 from matplotlib.colors import ListedColormap
 import matplotlib.patches as mpatches
-import cartopy.crs as ccrs
 #+end_src
 
 #+RESULTS:
@@ -48,9 +46,11 @@ ncpu
 #+end_src
 
 #+RESULTS:
-: 7
+: 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. 
+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 one degree.
+
+A buffer can be useful if we want to avoid "edge effects", while computing distance to forest edge for example. One degree tiles are used to download the data from GEE in parallel.
 
 #+begin_src python
 start_time = time.time()
@@ -75,7 +75,7 @@ print('Execution time:', round(elapsed_time, 2), 'minutes')
 #+end_src
 
 #+RESULTS:
-: Execution time: 1.15 minutes
+: Execution time: 20.15 minutes
 
 * Transform multiband fcc raster in one band raster 
 
@@ -155,42 +155,29 @@ grid = geopandas.read_file(grid_gpkg)
 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_coarsen["band_data"].plot(ax=ax, cmap=color_map, add_colorbar=False)
-grid_image = grid.boundary.plot(ax=ax, color="grey", linewidth=0.5)
-borders_image = borders.boundary.plot(ax=ax, color="black", linewidth=0.5)
-buffer_image = buffer.boundary.plot(ax=ax, color="black", linewidth=0.5)
-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 700 :align center
-#+RESULTS:
-[[file:fcc.png]]
-
-* Test to plot with gdal and imshow                                :noexport:
-
-#+begin_src python :results graphics file output :file fcc_tmp.png
 with gdal.Open("out_tmf/fcc_tmf_coarsen.tif", gdal.GA_ReadOnly) as ds:
     raster_image = ds.ReadAsArray()
+    nrow, ncol = raster_image.shape
+    xmin, xres, _, ymax, _, yres = ds.GetGeoTransform()
+    extent = [xmin, xmin + xres * ncol, ymax + yres * nrow, ymax]
 
 # Plot
 fig = plt.figure()
-ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
-plt.imshow(raster_image, cmap=color_map)
+ax = plt.subplot(111)
+ax.imshow(raster_image, cmap=color_map, extent=extent,
+          resample=False)
 grid_image = grid.boundary.plot(ax=ax, color="grey", linewidth=0.5)
 borders_image = borders.boundary.plot(ax=ax, color="black", linewidth=0.5)
 buffer_image = buffer.boundary.plot(ax=ax, color="black", linewidth=0.5)
 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_tmp.png", bbox_inches="tight", dpi=100)
+fig.savefig("fcc.png", bbox_inches="tight", dpi=200)
 #+end_src
 
+#+attr_rst: :width 700 :align center
 #+RESULTS:
-[[file:fcc_tmp.png]]
+[[file:fcc.png]]
 
+Lines in black represent country borders and the 10 km buffer. One degree tiles in grey cover the whole buffer and were used to download the data in parallel.
 
 # End

docsrc/notebooks/large_country/large_country.rst

@@ -22,12 +22,10 @@ cores, n-1 cores will be used in parallel.
     import numpy as np
     from osgeo import gdal
     from geefcc import get_fcc, sum_raster_bands
-    import xarray as xr
     import geopandas
     import matplotlib.pyplot as plt
     from matplotlib.colors import ListedColormap
     import matplotlib.patches as mpatches
-    import cartopy.crs as ccrs
 
 We initialize Google Earth Engine.
 
@@ -46,10 +44,12 @@ We can compute the number of cores used for the computation.
 
 ::
 
-    7
+    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. 
+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 one degree.
+
+A buffer can be useful if we want to avoid “edge effects”, while computing distance to forest edge for example. One degree tiles are used to download the data from GEE in parallel.
 
 .. code:: python
 
@@ -73,7 +73,7 @@ We estimate the computation time to download 159 1-degree tiles using several co
 
 ::
 
-    Execution time: 1.15 minutes
+    Execution time: 20.15 minutes
 
 Transform multiband fcc raster in one band raster
 -------------------------------------------------
@@ -148,17 +148,26 @@ We plot the forest cover change map.
 
 .. code:: python
 
+    with gdal.Open("out_tmf/fcc_tmf_coarsen.tif", gdal.GA_ReadOnly) as ds:
+        raster_image = ds.ReadAsArray()
+        nrow, ncol = raster_image.shape
+        xmin, xres, _, ymax, _, yres = ds.GetGeoTransform()
+        extent = [xmin, xmin + xres * ncol, ymax + yres * nrow, ymax]
+
     # Plot
     fig = plt.figure()
-    ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())
-    raster_image = fcc_tmf_coarsen["band_data"].plot(ax=ax, cmap=color_map, add_colorbar=False)
+    ax = plt.subplot(111)
+    ax.imshow(raster_image, cmap=color_map, extent=extent,
+              resample=False)
     grid_image = grid.boundary.plot(ax=ax, color="grey", linewidth=0.5)
     borders_image = borders.boundary.plot(ax=ax, color="black", linewidth=0.5)
     buffer_image = buffer.boundary.plot(ax=ax, color="black", linewidth=0.5)
     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)
+    fig.savefig("fcc.png", bbox_inches="tight", dpi=200)
 
 .. image:: fcc.png
     :width: 700
     :align: center
+
+Lines in black represent country borders and the 10 km buffer. One degree tiles in grey cover the whole buffer and were used to download the data in parallel.

docsrc/reference.rst

@@ -2,7 +2,11 @@ Python API
 **********
 
 get_fcc
-=========
+=======
 
 .. autofunction:: geefcc.get_fcc
 
+sum_raster_bands
+================
+
+.. autofunction:: geefcc.sum_raster_bands

geefcc/sum_raster_bands.py

@@ -14,7 +14,8 @@ def sum_raster_bands(input_file, output_file="sum.tif",
 
     :param input file: Input file with several bands.
 
-    :param output_file: Output file with several bands.
+    :param output_file: Output file with one band corresponding to the
+        sum of the input bands.
 
     :param blk_rows: Number of rows for block. This is used to break
         lage raster files in several blocks of data that can be hold

geefcc/temp.py

@@ -1,14 +1,14 @@
-            # Adjusted extent
-            xmin = aoi[0]
-            xmax = aoi[2]
-            xmax_tap = xmin + math.ceil((xmax - xmin) / scale) * scale
-            ymin = aoi[1]
-            ymax = aoi[3]
-            ymax_tap = ymin + math.ceil((ymax - ymin) / scale) * scale
+# Adjusted extent
+xmin = aoi[0]
+xmax = aoi[2]
+xmax_tap = xmin + math.ceil((xmax - xmin) / scale) * scale
+ymin = aoi[1]
+ymax = aoi[3]
+ymax_tap = ymin + math.ceil((ymax - ymin) / scale) * scale
 
-            # projWin = [ulx, uly, lrx, lry]
-            gdal.Translate(output_file, vrt_file,
-                           maskBand=None,
-                           projWin=[xmin, ymax_tap, xmax_tap, ymin],
-                           creationOptions=copts,
-                           callback=cback)
+# projWin = [ulx, uly, lrx, lry]
+gdal.Translate(output_file, vrt_file,
+               maskBand=None,
+               projWin=[xmin, ymax_tap, xmax_tap, ymin],
+               creationOptions=copts,
+               callback=cback)