From c1e8a4b3b499233027dc14e50a50c97cf0af402f Mon Sep 17 00:00:00 2001
From: Emma Oceguera-Conchas <113361856+E-O-Conchas@users.noreply.github.com>
Date: Wed, 8 Nov 2023 10:33:57 +0100
Subject: [PATCH] The script has been updated for a use friendly usage and the
multiproccesing for efficient raster generation has been added
---
...mbine_mask_and_fragmentation_tiles_gdal.py | 90 +++++++++++++++++++
1 file changed, 90 insertions(+)
create mode 100644 2.combine_mask_and_fragmentation_tiles_gdal.py
diff --git a/2.combine_mask_and_fragmentation_tiles_gdal.py b/2.combine_mask_and_fragmentation_tiles_gdal.py
new file mode 100644
index 0000000..2c72067
--- /dev/null
+++ b/2.combine_mask_and_fragmentation_tiles_gdal.py
@@ -0,0 +1,90 @@
+import os
+from glob import glob
+import numpy as np
+import re
+from osgeo import gdal
+from multiprocessing import Pool
+
+# Define file paths for fragmentation and land rasters
+frag_path = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\frag_raster"
+land_path = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\land_raster"
+out_path = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\frag_land_combined"
+
+# Define a regular expression pattern to match tile filenames
+pattern = r'frag_tile_x(\d+)_y(\d+)_1.tif'
+
+# Get a list of tile filenames from the fragmentation raster folder
+tiles = []
+for file_path in glob(frag_path + os.sep + '*.tif'):
+ match = re.search(pattern, file_path)
+ if match:
+ x_tile = match.group(1)
+ y_tile = match.group(2)
+ tiles.append(f"x{x_tile}_y{y_tile}_1.tif")
+ else:
+ print(f"Skipping file: {file_path}")
+
+# Initialize GDAL
+gdal.AllRegister()
+
+def process_tile(t):
+ # Define the file paths for the current tile's fragmentation and land rasters
+ frag_p = os.path.join(frag_path, 'frag_tile_' + t)
+ land_p = os.path.join(land_path, 'mask_tile_' + t)
+
+ # Open the fragmentation and land rasters using GDAL
+ frag_ds = gdal.Open(frag_p)
+ land_ds = gdal.Open(land_p)
+
+ # Get the band objects
+ frag_band = frag_ds.GetRasterBand(1)
+ land_band = land_ds.GetRasterBand(1)
+
+ # Read the raster data as numpy arrays
+ frag_data = frag_band.ReadAsArray()
+ land_data = land_band.ReadAsArray()
+
+ # Combine the two rasters
+ out_data = np.where((frag_data == 9) & (land_data == 1), 9, land_data)
+
+ # Create the output raster
+ driver = gdal.GetDriverByName('GTiff')
+ out_ds = driver.Create(
+ os.path.join(out_path, 'lf_comb_' + t),
+ frag_ds.RasterXSize,
+ frag_ds.RasterYSize,
+ 1,
+ gdal.GDT_Byte,
+ ['COMPRESS=DEFLATE']
+ )
+
+ # Set the projection and geotransform
+ out_ds.SetProjection('EPSG:3035')
+ out_ds.SetGeoTransform(frag_ds.GetGeoTransform())
+
+ # Write the output data
+ out_band = out_ds.GetRasterBand(1)
+ out_band.SetNoDataValue(0)
+ out_band.WriteArray(out_data)
+
+ # Close the datasets
+ out_band = None
+ out_ds = None
+ frag_ds = None
+ land_ds = None
+
+# Process the tiles using multiprocessing to utilize multiple cores
+if __name__ == '__main__':
+ with Pool(processes=os.cpu_count()) as pool:
+ pool.map(process_tile, tiles)
+
+
+
+
+
+
+
+
+
+
+
--
GitLab