diff --git a/.gitignore b/.gitignore index f2e5287..54801d1 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ # Local website +tmp_outputs docs # Byte-compiled / optimized / DLL files diff --git a/geefcc/download_gadm.py b/geefcc/download_gadm.py new file mode 100644 index 0000000..58e22a2 --- /dev/null +++ b/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 ``_\\ . + + :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 diff --git a/geefcc/ee_tmf.py b/geefcc/ee_tmf.py new file mode 100644 index 0000000..54bd0bb --- /dev/null +++ b/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 diff --git a/geefcc/geeic2geotiff.py b/geefcc/geeic2geotiff.py new file mode 100644 index 0000000..bc2653a --- /dev/null +++ b/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 diff --git a/geefcc/geotiff_from_tiles.py b/geefcc/geotiff_from_tiles.py new file mode 100644 index 0000000..d681902 --- /dev/null +++ b/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 diff --git a/geefcc/get_extent_from_aoi.py b/geefcc/get_extent_from_aoi.py new file mode 100644 index 0000000..ec31c1c --- /dev/null +++ b/geefcc/get_extent_from_aoi.py @@ -0,0 +1,87 @@ +"""Get extent from aoi.""" + +import os + +from download_gadm import download_gadm +from make_grid import create_buffer +from get_vector_extent import get_vector_extent + +opj = os.path.join +opd = os.path.dirname + + +def get_extent_from_aoi(aoi, buff, out_dir): + """Get extent from aoi. + + :param aoi: Area of interest defined either by a country iso code + (three letters), a vector file, or an extent in lat/long + (tuple with (xmin, ymin, xmax, ymax)). + + :param buff: Buffer around the aoi. In decimal degrees + (e.g. 0.08983152841195216 correspond to ~10 km at the + equator). + + :out_dir: Output directory. + + :return: A dictionary including extent and border vector file. + """ + + # Set aoi_isfile + aoi_isfile = True + + # aoi = country iso code + if isinstance(aoi, str) and len(aoi) == 3: + # Download borders + iso = aoi + borders_gpkg = opj(out_dir, f"gadm41_{iso}_0.gpkg") + download_gadm(iso, output_file=borders_gpkg) + # Buffer around borders + if buff > 0: + buff_file = opj( + out_dir, + f"gadm41_{iso}_buffer.gpkg") + create_buffer(input_file=borders_gpkg, + output_file=buff_file, + buffer_dist=buff) + borders_gpkg = buff_file + # Extent + extent_latlong = get_vector_extent(borders_gpkg) + + # aoi = extent + elif isinstance(aoi, tuple) and len(aoi) == 4: + aoi_isfile = False + # nb: We could create a vector file here... + borders_gpkg = None + if buff > 0: + extent_latlong = (aoi[0] - buff, aoi[1] - buff, + aoi[2] + buff, aoi[3] + buff) + else: + extent_latlong = aoi + + # aoi = gpkg file + elif os.path.isfile(aoi) and aoi[-5:] == ".gpkg": + # Buffer around borders + if buff > 0: + buff_file = opj( + out_dir, + "borders_buffer.gpkg") + create_buffer(input_file=aoi, + output_file=buff_file, + buffer_dist=buff) + borders_gpkg = buff_file + else: + borders_gpkg = aoi + # Extent + extent_latlong = get_vector_extent(borders_gpkg) + + # Else raise error + else: + raise ValueError("aoi must be either a country iso code, " + "an extent, or a gpkg file") + + return {"extent_latlong": extent_latlong, + "borders_gpkg": borders_gpkg, + "aoi_isfile": aoi_isfile} + + +# End diff --git a/geefcc/get_fcc.py b/geefcc/get_fcc.py index 994a7cc..aeddf9d 100644 --- a/geefcc/get_fcc.py +++ b/geefcc/get_fcc.py @@ -1,11 +1,113 @@ """Get forest cover change data.""" -# Import -# import os +import os +import multiprocessing as mp +from get_extent_from_aoi import get_extent_from_aoi +from make_dir import make_dir +from make_grid import make_grid, grid_intersection +from geeic2geotiff import geeic2geotiff +from ee_tmf import ee_tmf +from geotiff_from_tiles import geotiff_from_tiles -def get_fcc(): - """Get forest cover change data.""" - return "get_fcc" +opj = os.path.join +opd = os.path.dirname + + +def get_fcc(ee_project, + aoi, + buff=0, + years=[2000, 2010, 2020], + source="tmf", + perc=75, + tile_size=1, + output_file="fcc.tiff"): + """Get forest cover change data. + + :param ee_project: Earth Engine project name. + + :param aoi: Area of interest defined either by a country iso code + (three letters), a vector file, or an extent in lat/long + (tuple with (xmin, ymin, xmax, ymax)). + + :param buff: Buffer around the aoi. In decimal degrees + (e.g. 0.08983152841195216 correspond to ~10 km at the + equator). + + :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. + + :param source: Either "gfc" for Global Forest Change or "tmf" for + Tropical Moist Forest. If "gfc", the tree cover threshold + defining the forest must be specified with parameter `perc`. + + :param perc: Tree cover threshold defining the forest. + + :param tile_size: Tile size for parallel computing. + + :output_file: Path to output GeoTIFF file. One band for each + year. Value 1 for forest and 0 for non-forest. + + """ + + # Output dir + out_dir = opd(output_file) + + # Variables + proj = "EPSG:4326" + epsg_code = 4326 + 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"] + + # Make minimal grid + grid_gpkg = opj(out_dir, "grid.gpkg") + grid = make_grid(extent_latlong, buff=0, tile_size=tile_size, + scale=scale, proj=epsg_code, ofile=grid_gpkg) + if aoi_isfile: + min_grid = opj(out_dir, "min_grid.gpkg") + grid_i = grid_intersection(grid, grid_gpkg, min_grid, + borders_gpkg) + # Update grid and file + grid = grid_i + grid_gpkg = min_grid + + # Forest image collection + if source == "tmf": + forest = ee_tmf(ee_project, years) + + # Create dir for forest tiles + out_dir_tiles = opj(out_dir, "forest_tiles") + make_dir(out_dir_tiles) + + # Write tiles in parallel + ncpu = os.cpu_count() - 1 + pool = mp.Pool(processes=ncpu) + args = [(i, ext, forest, proj, scale, out_dir_tiles) + for (i, ext) in enumerate(grid)] + pool.starmap_async(geeic2geotiff, args).get() + pool.close() + pool.join() + + # Geotiff from tiles + geotiff_from_tiles(extent_latlong, scale, out_dir) + + + +# # Test +# ee_project = "forestatrisk" +# aoi = "REU" +# buff = 0.08983152841195216 +# years = [2000, 2010, 2020] +# source = "tmf" +# perc = 75 +# tile_size = 0.5 +# output_file = "tmp_outputs/fcc.tiff" # End diff --git a/geefcc/get_vector_extent.py b/geefcc/get_vector_extent.py new file mode 100644 index 0000000..1aaa403 --- /dev/null +++ b/geefcc/get_vector_extent.py @@ -0,0 +1,29 @@ +"""Get the extent of a shapefile.""" + +from osgeo import ogr + + +def get_vector_extent(input_file): + """Compute the extent of a vector file. + + This function computes the extent (xmin, ymin, xmax, ymax) of a + shapefile. + + :param input_file: Path to the input vector file. + + :return: The extent as a tuple (xmin, ymin, xmax, ymax). + + """ + + in_data_dource = ogr.Open(input_file) + in_layer = in_data_dource.GetLayer() + extent = in_layer.GetExtent() + extent = (extent[0], extent[2], extent[1], extent[3]) + + return extent # (xmin, ymin, xmax, ymax) + + +# Function alias for compatibility with old versions of forestatrisk. +extent_shp = get_vector_extent + +# End diff --git a/geefcc/make_dir.py b/geefcc/make_dir.py new file mode 100644 index 0000000..eb375cf --- /dev/null +++ b/geefcc/make_dir.py @@ -0,0 +1,29 @@ +"""Make new directory.""" + +import os + + +def make_dir(newdir): + """Make new directory. + + * Already exists, silently complete. + * Regular file in the way, raise an exception. + * Parent directory(ies) does not exist, make them as well. + + :param newdir: Directory path to create. + + """ + if os.path.isdir(newdir): + pass + elif os.path.isfile(newdir): + raise OSError( + ("A file with the same name as the desired " + f"dir, '{newdir}', already exists.") + ) + else: + head, tail = os.path.split(newdir) + if head and not os.path.isdir(head): + make_dir(head) + # print "_mkdir %s" % repr(newdir) + if tail: + os.mkdir(newdir) diff --git a/geefcc/make_grid.py b/geefcc/make_grid.py new file mode 100644 index 0000000..e356f1e --- /dev/null +++ b/geefcc/make_grid.py @@ -0,0 +1,192 @@ +"""Make minimal grid with buffer around polygons.""" + +import os + +import numpy as np +from osgeo import ogr, osr + + +def create_buffer(input_file, output_file, buffer_dist): + """Create buffer. + + Make buffers around features of a layer and saves them to a new + layer. + + Source: https://pcjericks.github.io/py-gdalogr-cookbook/ + vector_layers.html#create-buffer + + :param input_file: Input filename. + :param output_file: Output filename (`.gpkg`). + :param buffer_dist: Buffer distance (in unit of CRS). + + """ + input_ds = ogr.Open(input_file) + input_lyr = input_ds.GetLayer() + + driver = ogr.GetDriverByName("GPKG") + if os.path.exists(output_file): + driver.DeleteDataSource(output_file) + ds = driver.CreateDataSource(output_file) + lyr = ds.CreateLayer( + "buffer", + geom_type=ogr.wkbPolygon) + feature_defn = lyr.GetLayerDefn() + + for feature in input_lyr: + in_geom = feature.GetGeometryRef() + geom_buffer = in_geom.Buffer(buffer_dist) + + out_feature = ogr.Feature(feature_defn) + out_feature.SetGeometry(geom_buffer) + lyr.CreateFeature(out_feature) + out_feature = None + + +def gpkg_from_grid(grid, proj=4326, ofile="grid.gpkg"): + """Make vector file from grid. + + Source: https://pcjericks.github.io/py-gdalogr-cookbook/ + vector_layers.html#create-fishnet-grid + + :param grid: List of extents for each grid cell. + :param proj: Projection as EPSG code, default to 4326. + :param ofile: Output file, default to `grid.gpkg`. + + """ + + # Set up the shapefile driver + driver = ogr.GetDriverByName("GPKG") + + # Create the data source + if os.path.exists(ofile): + os.remove(ofile) + ds = driver.CreateDataSource(ofile) + + # Create the spatial reference system, WGS84 + srs = osr.SpatialReference() + srs.ImportFromEPSG(proj) + + # Create one layer + layer = ds.CreateLayer("grid", srs, ogr.wkbPolygon) + + # Add an ID field + id_field = ogr.FieldDefn("id", ogr.OFTInteger) + layer.CreateField(id_field) + + # Feature definition + feature_def = layer.GetLayerDefn() + + # Create grid cells + for (i, coords) in enumerate(grid): + # Get coordinates + xmin = coords[0] + ymin = coords[1] + xmax = coords[2] + ymax = coords[3] + # Create geometry + ring = ogr.Geometry(ogr.wkbLinearRing) + ring.AddPoint(xmin, ymax) + ring.AddPoint(xmax, ymax) + ring.AddPoint(xmax, ymin) + ring.AddPoint(xmin, ymin) + ring.AddPoint(xmin, ymax) + poly = ogr.Geometry(ogr.wkbPolygon) + poly.AddGeometry(ring) + # Add geometry to layer + feature = ogr.Feature(feature_def) + feature.SetGeometry(poly) + feature.SetField("id", i) + layer.CreateFeature(feature) + feature = None + + # Return + ds = None + + +def make_grid(extent, buff, tile_size, scale, proj=4326, + ofile="grid.gpkg"): + """Make overlapping grid from an extent and resolution. + + :param extent: Extent (xmin, ymin, xmax, ymax). + :param buff: Buffer (same unit as extent). + :param tile_size: Tile size (same unit as extent). + :param scale: Resolution (same unit as extent). + :param proj: Projection as EPSG code, default to 4326. + :param ofile: Output file, default to `grid.gpkg`. + + :return: List of extents for each grid cell. + + """ + + # Buffer around extent + xmin = extent[0] - buff + ymin = extent[1] - buff + xmax = extent[2] + buff + ymax = extent[3] + buff + # Adapt tile_size to scale + tile_size = int(np.round(tile_size / scale)) * scale + # List of x coordinates + xlist = list(np.arange(xmin, xmax + tile_size, tile_size)) + nx = len(xlist) + # List of y coordinates + ylist = list(np.arange(ymin, ymax + tile_size, tile_size)) + ny = len(ylist) + # Grid: list of extents + grid = [(xlist[i], ylist[j], xlist[i + 1], ylist[j + 1]) + for i in range(nx - 1) for j in range(ny - 1)] + # Create vector file from grid + gpkg_from_grid(grid, proj, ofile) + return grid + + +def grid_intersection(grid, input_grid, output_grid, borders_gpkg): + """Grid intersection. + + :param grid: List of extents for grid cells. + :param input_grid: Input grid vector file. + :param output_grid: Output grid vector file. + :param borders_gpkg: Border vector file. + + :return: List of extents for intersecting grid cells. + + """ + # Grid + dr_g = ogr.GetDriverByName("GPKG") + ds_g = dr_g.Open(input_grid) + lay_g = ds_g.GetLayer() + # Borders + dr_b = ogr.GetDriverByName("GPKG") + ds_b = dr_b.Open(borders_gpkg) + lay_b = ds_b.GetLayer() + # New grid + grid_i = [] + if os.path.exists(output_grid): + os.remove(output_grid) + ds = dr_g.CreateDataSource(output_grid) + wkt = lay_g.GetSpatialRef().ExportToWkt() + srs = osr.SpatialReference() + srs.ImportFromWkt(wkt) + layer = ds.CreateLayer("grid_i", srs, ogr.wkbPolygon) + defn = lay_g.GetLayerDefn() + for i in range(defn.GetFieldCount()): + layer.CreateField(defn.GetFieldDefn(i)) + # Loop on features + for (ext, feat_g) in zip(grid, lay_g): + geom_g = feat_g.GetGeometryRef() + for feat_b in lay_b: + geom_b = feat_b.GetGeometryRef() + if geom_g.Intersects(geom_b): + grid_i.append(ext) + layer.CreateFeature(feat_g) + # Reset reading so that features of lay_b + # are accessible again + lay_b.ResetReading() + break + # Clean + ds = None + ds_b = None + ds_g = None + # Return + return grid_i + +# End Of File