First functioning get_ffc.py function

.gitignore

@@ -1,4 +1,5 @@
 # Local website
+tmp_outputs
 docs
 
 # Byte-compiled / optimized / DLL files

geefcc/download_gadm.py

@@ -0,0 +1,28 @@
+"""Download GADM data."""
+
+import os
+from urllib.request import urlretrieve
+
+
+def download_gadm(iso3, output_file):
+    """Download GADM data for a country.
+
+    Download GADM (Global Administrative Areas) for a specific
+    country. See `<https://gadm.org>`_\\ .
+
+    :param iso3: Country ISO 3166-1 alpha-3 code.
+
+    :param output_file: Path to output GPKG file.
+
+    """
+
+    # Check for existing file
+    if not os.path.isfile(output_file):
+
+        # Download the zipfile from gadm.org
+        url = ("https://geodata.ucdavis.edu/gadm/gadm4.1/"
+               f"gpkg/gadm41_{iso3}.gpkg")
+        urlretrieve(url, output_file)
+
+
+# End

geefcc/ee_tmf.py

@@ -0,0 +1,77 @@
+"""Compute fcc on GEE using the TMF product."""
+
+import ee
+
+
+def ee_tmf(ee_project, years):
+    """Compute fcc on GEE using the TMF product.
+
+    :param ee_project: Earth Engine project name.
+
+    :param years: List of years defining time-periods for estimating
+        forest cover change. Years for computing forest cover change
+        can be in the interval 2001--2024 for GFC (GFC does not
+        provide loss for the year 2000) and 2000--2023 for TMF.
+
+    :return: An image collection for forest where each image
+        correspond to a year.
+
+    """
+
+    # Initialize Earth Engine
+    ee.Initialize(project=ee_project,
+                  opt_url=("https://earthengine-highvolume"
+                           ".googleapis.com"))
+
+    # Get annual product
+    annual_product = ee.ImageCollection(
+        "projects/JRC/TMF/"
+        "v1_2022/AnnualChanges")
+    annual_product = annual_product.mosaic().toByte()
+
+    # ap_all_year: forest if Y = 1 or 2.
+    ap_forest = annual_product.where(annual_product.eq(2), 1)
+    ap_all_year = ap_forest.where(ap_forest.neq(1), 0)
+
+    forest_list = []
+    band_names = []
+
+    for year in years:
+        id_year = year - 1990 - 1
+        ap = ap_all_year.select(list(range(id_year, 33)))
+        forest_yr = ap.reduce(ee.Reducer.sum()).gte(1)
+        forest_yr = forest_yr.set(
+            "system:time_start",
+            ee.Date.fromYMD(year, 1, 1).millis())
+        forest_list.append(forest_yr)
+        band_names.append(f"forest{year}")
+
+    forest_ic = ee.ImageCollection(forest_list)
+
+    def get_date(image):
+        """Get formatted date."""
+        date = ee.Image(image).date().format("YYYY-MM-dd")
+        return date
+
+    def filter_and_mosaic(d):
+        """Create mosaic for one date."""
+        d = ee.Date(d)
+        im = (forest_ic
+              .filterDate(d, d.advance(1, "day"))
+              .mosaic().toByte())
+        im = im.set("system:time_start", d.millis(),
+                    "system:id", d.format("YYYY-MM-dd"))
+        return im
+
+    def mosaic_by_date(img_list):
+        """Mosaic by date."""
+        unique_dates = img_list.map(get_date).distinct()
+        mosaic_list = unique_dates.map(filter_and_mosaic)
+        return ee.ImageCollection(mosaic_list)
+
+    forest = mosaic_by_date(ee.List(forest_list))
+
+    return forest
+
+
+# End

geefcc/geeic2geotiff.py

@@ -0,0 +1,186 @@
+"""Saves an xarray dataset to a Cloud Optimized GeoTIFF (COG).
+
+Adapted from:
+https://gist.github.com/GerardoLopez/35123d4a15aa31f3ea4b01efb5b26d4d
+"""
+
+import os
+
+from osgeo import gdal, osr, gdal_array
+import xarray as xr
+
+
+def get_dst_dataset(dst_img, cols, rows, layers, dtype, proj, gt):
+    """Create a GDAL data set in Cloud Optimized GeoTIFF (COG) format.
+
+    :param dst_img: Output filenane full path
+    :param cols: Number of columns
+    :param rows: Number of rows
+    :param layers: Number of layers
+    :param dtype: GDAL type code
+    :param proj: Projection information in WKT format
+    :param gt: GeoTransform tupple
+
+    :return dst_ds: GDAL destination dataset object
+
+    """
+
+    gdal.UseExceptions()
+    try:
+        # Default driver options to create a COG
+        driver = gdal.GetDriverByName('GTiff')
+        driver_options = ['COMPRESS=DEFLATE',
+                          'PREDICTOR=1',
+                          'BIGTIFF=YES',
+                          'TILED=YES',
+                          'COPY_SRC_OVERVIEWS=YES']
+
+        # Create driver
+        dst_ds = driver.Create(dst_img, cols, rows, layers,
+                               dtype, driver_options)
+
+        # Set cartographic projection
+        dst_ds.SetProjection(proj)
+        dst_ds.SetGeoTransform(gt)
+
+    except Exception as err:
+        if err.err_level >= gdal.CE_Warning:
+            # print('Cannot write dataset: %s' % self.input.value)
+            # Stop using GDAL exceptions
+            gdal.DontUseExceptions()
+            raise RuntimeError(err.err_level, err.err_no, err.err_msg)
+
+    gdal.DontUseExceptions()
+    return dst_ds
+
+
+def get_resolution_from_xarray(xarray):
+    """Method to create a tuple (x resolution, y resolution) in x and y
+    dimensions.
+
+    :param xarray: xarray with latitude and longitude variables.
+
+    :return: tuple with x and y resolutions
+    """
+
+    x_res = xarray.longitude.values[1] - xarray.longitude.values[0]
+    y_res = xarray.latitude.values[0] - xarray.latitude.values[1]
+
+    return (x_res, y_res)
+
+
+def xarray2geotiff(xarray, data_var, out_dir, index):
+    """Saves an xarray dataset to a Cloud Optimized GeoTIFF (COG).
+
+    :param xarray: xarray Dataset.
+    :param data_var: Data variable in the xarray dataset.
+    :param out_dir: Output directory.
+    :param index: Tile index.
+
+    """
+
+    # Create GeoTransform - perhaps the user requested a
+    # spatial subset, therefore is compulsory to update it
+
+    # GeoTransform -- case of a "north up" image without
+    #                 any rotation or shearing
+    #  GeoTransform[0] top left x
+    #  GeoTransform[1] w-e pixel resolution
+    #  GeoTransform[2] 0
+    #  GeoTransform[3] top left y
+    #  GeoTransform[4] 0
+    #  GeoTransform[5] n-s pixel resolution (negative value)
+
+    # Reorganize the data
+    xarray = xarray.transpose("time", "latitude", "longitude")
+
+    # Create tmp xarray DataArray
+    _xarray = getattr(xarray, data_var)
+
+    x_res, y_res = get_resolution_from_xarray(_xarray)
+
+    gt = (_xarray.longitude.data[0] - (x_res / 2.),
+          x_res, 0.0,
+          _xarray.latitude.data[-1] - (y_res / 2.),
+          0.0, y_res)
+
+    # Coordinate Reference System (CRS) in a PROJ4 string to a
+    # Spatial Reference System Well Known Text (WKT)
+    crs = xarray.attrs['crs']
+    crs = int(crs[5:9])
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(crs)
+    proj = srs.ExportToWkt()
+
+    # Get GDAL datatype from NumPy datatype
+    if _xarray.dtype == 'bool':
+        dtype = gdal.GDT_Byte
+    else:
+        dtype = gdal_array.NumericTypeCodeToGDALTypeCode(_xarray.dtype)
+
+    # File name with index
+    fname = os.path.join(out_dir, f"forest_{index}.tif")
+
+    # Dimensions
+    layers, rows, cols = _xarray.shape
+
+    # Create destination dataset
+    dst_ds = get_dst_dataset(
+        dst_img=fname, cols=cols, rows=rows,
+        layers=layers, dtype=dtype, proj=proj, gt=gt)
+
+    for layer in range(layers):
+        dst_band = dst_ds.GetRasterBand(layer + 1)
+
+        # Date
+        if 'time' in _xarray.dims:
+            dst_band.SetMetadataItem(
+                'time',
+                _xarray.time.data[layer].astype(str))
+
+        # Data variable name
+        dst_band.SetMetadataItem('data_var', data_var)
+
+        # Data
+        data_npa = _xarray[layer].data
+        index = list(reversed(range(rows)))
+        data_npa = data_npa[index]
+        dst_band.WriteArray(data_npa)
+
+    dst_band = None
+    del dst_ds
+
+
+def geeic2geotiff(index, extent, forest, proj, scale, out_dir):
+    """Write a GEE image collection to a geotiff.
+
+    :param index: Tile index.
+    :param extent: Tile extent.
+    :param forest: Forest image collection from GEE.
+    :param proj: Projection.
+    :param scale: Scale.
+    :param output_dir: Output directory.
+
+    """
+
+    # 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"})
+        .rename({"sum": "forest_cover"})
+    )
+    ds = ds.load()
+
+    # Load and write data to geotiff
+    xarray2geotiff(ds, "forest_cover", out_dir, index)
+
+
+# End Of File

geefcc/geotiff_from_tiles.py

@@ -0,0 +1,63 @@
+"""Make geotiff from tiles."""
+
+import os
+from glob import glob
+import math
+
+from osgeo import gdal
+
+opj = os.path.join
+opd = os.path.dirname
+
+
+def geotiff_from_tiles(extent_latlong, scale, out_dir):
+    """Make geotiff from tiles.
+
+    :param extent_latlong: Extent in lat/long.
+    :param scale: Resolution.
+    :param out_dir: Output directory.
+
+    """
+
+    # Dir for forest tiles
+    out_dir_tiles = opj(out_dir, "forest_tiles")
+
+    # Make vrt
+    tif_forest_files = glob(opj(out_dir_tiles, "forest_*.tif"))
+    # Callback
+    verbose = False
+    cback = gdal.TermProgress if verbose else 0
+    forest_vrt = gdal.BuildVRT(
+        opj(out_dir, "forest.vrt"),
+        tif_forest_files,
+        callback=cback)
+    # Flush cache
+    # https://gis.stackexchange.com/questions/306664/
+    # gdal-buildvrt-not-creating-any-output
+    # https://gis.stackexchange.com/questions/44003/
+    # python-equivalent-of-gdalbuildvrt
+    forest_vrt.FlushCache()
+    forest_vrt = None
+
+    # VRT to GeoTIFF
+    # Creation options
+    copts = ["COMPRESS=DEFLATE", "BIGTIFF=YES"]
+    # Adjusted extent
+    xmin = extent_latlong[0]
+    xmax = extent_latlong[2]
+    xmax_tap = xmin + math.ceil((xmax - xmin) / scale) * scale
+    ymin = extent_latlong[1]
+    ymax = extent_latlong[3]
+    ymax_tap = ymin + math.ceil((ymax - ymin) / scale) * scale
+
+    # Call to gdal_translate
+    ifile = opj(out_dir, "forest.vrt")
+    ofile = opj(out_dir, "fcc.tif")
+    # projWin = [ulx, uly, lrx, lry]
+    gdal.Translate(ofile, ifile,
+                   maskBand=None,
+                   projWin=[xmin, ymax_tap, xmax_tap, ymin],
+                   creationOptions=copts,
+                   callback=cback)
+
+# End