""" Title: Semi-Natural Habitat (SNH) Layer Processing Script Description: This script processes geospatial raster data to generate a Semi-Natural Habitat (SNH) layer, which is used for analyzing ecological and environmental indicators. It takes input rasters for Morphological Spatial Pattern Analysis (MSPA), extensive grasslands, and CORINE land cover data, reclassifies and combines them, applies spatial masking, and produces a final output raster. Key Features: - Splits large rasters into smaller tiles for efficient processing. - Reclassifies MSPA tiles into categories (Core, Edge, Linear). - Combines MSPA and extensive grasslands data into an SNH layer. - Masks the SNH layer using a binary CORINE raster. - Merges processed tiles back into a single raster file. - Handles user-defined input and output paths dynamically. Usage: 1. Define paths to the project directory and input rasters. 2. Uncomment the desired steps in the `main` function based on the required processing. 3. Run the script in an environment with required libraries installed (GDAL, NumPy, Rasterio, etc.). 4. The script saves intermediate and final outputs to specified directories. Author: Rui Catarino Date: 16/01/25 """ import os import sys from pathlib import Path import subprocess import time import geopandas as gpd from shapely.geometry import box from osgeo import gdal, ogr import numpy as np import rasterio from scipy.ndimage import distance_transform_edt from scipy.ndimage import convolve import matplotlib.pyplot as plt import pandas as pd from rasterstats import zonal_stats from affine import Affine # User-defined inputs project_path = input("Enter the root project directory: ").strip() or "/path/to/project" output_path = os.path.join(project_path, "Output_Data") temp_path = os.path.join(output_path, "Temp_SNH") # Input files MSPA_tif = input("Enter the path to the MSPA raster file: ").strip() or "/path/to/MSPA.tif" CORINE_tif = input("Enter the path to the CORINE raster file: ").strip() or "/path/to/CORINE.tif" extGL_tif = input("Enter the path to the extensive grasslands raster file: ").strip() or "/path/to/extGL.tif" # Output files corine_bin_tif = os.path.join(output_path, "U2018_CLC2018_10bin.tif") SNH_tif = os.path.join(output_path, "SNH_2018.tif") SNHmasked_tif = os.path.join(output_path, "SNH_masked.tif") # Temporary directories MSPA_tiles_path = os.path.join(temp_path, "MSPA_tiles") recMSPA_tiles_path = os.path.join(temp_path, "MSPA_rec_tiles") extGL_tiles_path = os.path.join(temp_path, "extGL_tiles") SNH_tiles_path = os.path.join(temp_path, "SNH_tiles") # Ensure directories exist Path(MSPA_tiles_path).mkdir(parents=True, exist_ok=True) Path(recMSPA_tiles_path).mkdir(parents=True, exist_ok=True) Path(extGL_tiles_path).mkdir(parents=True, exist_ok=True) Path(SNH_tiles_path).mkdir(parents=True, exist_ok=True) # Report time function def report_time(start_time, step_name): elapsed_time = time.time() - start_time hours, remainder = divmod(elapsed_time, 3600) minutes, seconds = divmod(remainder, 60) print(f"\n{step_name} completed in {int(hours)}h {int(minutes)}m {int(seconds)}s.") # Function to check if an output file exists and is valid (not corrupted) def check_output_tif(output_raster_path): """ Check if the output file exists and whether it is corrupted. Parameters: output_raster_path (str): Path to the output raster file. Returns: bool: True if the file exists and is not corrupted, False if the file is corrupted or doesn't exist. """ if os.path.exists(output_raster_path): try: # Try opening the raster to check if it is readable ds = gdal.Open(output_raster_path) if ds is not None: print(f"Output file {output_raster_path} already exists and is valid. Skipping...") return True # File exists and is valid else: print(f"Output file {output_raster_path} is corrupted. Reprocessing...") os.remove(output_raster_path) # Delete corrupted file return False # File is corrupted except Exception as e: print(f"Error while checking {output_raster_path}: {e}") print(f"File {output_raster_path} is corrupted. Deleting it.") os.remove(output_raster_path) # Delete corrupted file return False # File is corrupted or unreadable return False # File doesn't exist # Function to save array to a GeoTIFF file def save_array_to_tiff(array, output_path, driver, mspa_ds): out_ds = driver.Create(str(output_path), mspa_ds.RasterXSize, mspa_ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetProjection(mspa_ds.GetProjection()) out_ds.SetGeoTransform(mspa_ds.GetGeoTransform()) out_ds.GetRasterBand(1).WriteArray(array) out_ds.FlushCache() out_ds = None # Close and save the file # Define function to print the extent using rasterio def print_tile_extent(tile_path): """ Function to print the extent of a given tile using rasterio. Parameters: tile_path (str): Path to the tile file. """ with rasterio.open(tile_path) as ds: bounds = ds.bounds print(f"Extent of tile {tile_path}:") print(f" X Min: {bounds.left}") print(f" X Max: {bounds.right}") print(f" Y Min: {bounds.bottom}") print(f" Y Max: {bounds.top}") # Function to Convert the NumPy array back into a GDAL dataset before using gdal.Translate. def array_to_gdal_dataset(array, affine_transform, projection, nodata_value=None): """ Converts a NumPy array to a GDAL dataset in memory. """ driver = gdal.GetDriverByName('MEM') # In-memory GDAL driver rows, cols = array.shape dataset = driver.Create('', cols, rows, 1, gdal.GDT_Float32) # Create a dataset with 1 band # Set the affine transformation and projection dataset.SetGeoTransform(affine_transform) dataset.SetProjection(projection) # Write the array to the dataset's first band band = dataset.GetRasterBand(1) band.WriteArray(array) # Set NoData value if provided if nodata_value is not None: band.SetNoDataValue(nodata_value) band.FlushCache() # Ensure changes are written to the dataset return dataset # Function to crop tiles def create_tiles(mspa_rst, gl_rst, MSPA_tiles_path, gl_tiles_path, y2): """ Cuts both mspa_rst and gl_rst into horizontal strips (tiles) without buffer and saves them to output folders. Parameters: mspa_rst (str): Path to the MSPA raster file. gl_rst (str): Path to the gl raster file. MSPA_tiles_path (str): Path to the folder where MSPA tiles will be saved. gl_tiles_path (str): Path to the folder where gl tiles will be saved. y2 (int): Height of each tile in kilometers. """ step_start_time = time.time() print("\n" + "="*40) print(f"Processing input rasters: \n - {mspa_rst} \n - {gl_rst} \n") # Load rasters mspa_ds = gdal.Open(mspa_rst) gl_ds = gdal.Open(gl_rst) # Get geotransforms (top-left corner and pixel size) and raster size for rasters mspa_geotransform = mspa_ds.GetGeoTransform() gl_geotransform = gl_ds.GetGeoTransform() # Calculate extents (width and height in pixels) mspa_x_extent = mspa_ds.RasterXSize * mspa_geotransform[1] gl_x_extent = gl_ds.RasterXSize * gl_geotransform[1] x2 = max(mspa_x_extent, gl_x_extent) raster_height = max(mspa_ds.RasterYSize, gl_ds.RasterYSize) y2 = y2 * 1000 print(f"MSPA Geotransform: {mspa_geotransform}") print(f"Grassland Geotransform: {gl_geotransform}") mspa_extent = ( mspa_geotransform[0], mspa_geotransform[3], mspa_geotransform[0] + mspa_ds.RasterXSize * mspa_geotransform[1], mspa_geotransform[3] + mspa_ds.RasterYSize * mspa_geotransform[5] ) gl_extent = ( gl_geotransform[0], gl_geotransform[3], gl_geotransform[0] + gl_ds.RasterXSize * gl_geotransform[1], gl_geotransform[3] + gl_ds.RasterYSize * gl_geotransform[5] ) print(f"MSPA Extent: {mspa_extent}") print(f"Grassland Extent: {gl_extent}") # Calculate number of vertical cuts based on y2 and raster height y_cuts = int(raster_height / y2) print(f"Calculated number of cuts (height-wise): {y_cuts}") for i in range(y_cuts): tile_id = f"tile_{i}.tif" # Set the y offsets for each tile y1 = i * y2 # File paths for MSPA and gl tiles mspa_path = os.path.join(MSPA_tiles_path, tile_id) gl_path = os.path.join(gl_tiles_path, tile_id) width_x = mspa_ds.RasterXSize process_mspa = not check_output_tif(str(mspa_path)) process_gl = not check_output_tif(str(gl_path)) if process_mspa: command_mspa = [ "gdal_translate", "-of", "GTiff", "-co", "BIGTIFF=YES", "-co", "TILED=YES", "-srcwin", "0", str(y1), str(width_x), str(y2), mspa_rst, mspa_path ] print(f"Processing tile {i + 1}/{y_cuts} for MSPA") subprocess.run(command_mspa, check=True) if process_gl: command_gl = [ "gdal_translate", "-of", "GTiff", "-co", "BIGTIFF=YES", "-srcwin", "0", str(y1), str(width_x), str(y2), gl_rst, gl_path ] print(f"Processing tile {i + 1}/{y_cuts} for Grasslands") subprocess.run(command_gl, check=True) print(f"Tiles saved to -{MSPA_tiles_path} \n -{gl_tiles_path}") report_time(step_start_time, "Tiles creation") print("="*40) # Function to reclassify MSPA tiles def recla_MSPA(MSPA_tiles_path, recMSPA_tiles_path, start_feature=0): """ Reclassify MSPA tiles into different classes (Core, Edge, Linear), adjust core scores based on distance to edges, and save the reclassified tiles. """ step_start_time = time.time() print("\n" + "="*40) # A line of "=" for separation print(f"Reclassifying MSPA tiles in classes (Core, Edge, Linear), ...") # Get the list of tiles from MSPA directory mspa_tiles = list(Path(MSPA_tiles_path).glob("*.tif")) total_tiles = len(mspa_tiles) # Initialize GDAL driver once driver = gdal.GetDriverByName('GTiff') for i, mspa_tile in enumerate(mspa_tiles): if i < start_feature: continue # Skip tiles until reaching the start_feature index tile_id = mspa_tile.stem.split("_")[-1] print(f"MSPA reclassification, processing tile {tile_id} [{i+1}/{total_tiles}] ...") # Check if the output file exists and is not corrupted output_mspa_tile_path = Path(recMSPA_tiles_path) / f"tile_{tile_id}.tif" # Set the output file path if check_output_tif(str(output_mspa_tile_path)): continue # Skip if the file exists and is valid # Open MSPA tile mspa_ds = gdal.Open(str(mspa_tile)) mspa_array = mspa_ds.GetRasterBand(1).ReadAsArray() # Ensure the array is float32 to handle the reclassification properly mspa_array = mspa_array.astype(np.float32) ### Reclassify Core class (Woody areal interior) ### core_array = np.copy(mspa_array) core_array[(core_array == 1)] = 0 core_array[(core_array == 17)] = 1 core_array[(core_array == 117)] = 1 core_array[(core_array != 1)] = 0 # All non-Core pixels set to 0 ### Reclassify Edge class (Woody areal exterior) ### edge_array = np.copy(mspa_array) edge_array[(edge_array == 2)] = 0 edge_array[(edge_array == 3)] = 2 edge_array[(edge_array == 5)] = 2 edge_array[(edge_array == 9)] = 2 edge_array[(edge_array == 35)] = 2 edge_array[(edge_array == 67)] = 2 edge_array[(edge_array == 103)] = 2 edge_array[(edge_array == 105)] = 2 edge_array[(edge_array == 119)] = 2 edge_array[(edge_array == 135)] = 2 edge_array[(edge_array == 137)] = 2 edge_array[(edge_array == 167)] = 2 edge_array[(edge_array == 169)] = 2 edge_array[(edge_array != 2)] = 0 # All non-Edge pixels set to 0 ### Reclassify Linear class (Woody linear) ### linear_array = np.copy(mspa_array) linear_array[(linear_array == 3)] = 0 linear_array[(linear_array == 1)] = 3 linear_array[(linear_array == 33)] = 3 linear_array[(linear_array == 37)] = 3 linear_array[(linear_array == 65)] = 3 linear_array[(linear_array == 69)] = 3 linear_array[(linear_array == 101)] = 3 linear_array[(linear_array == 109)] = 2 linear_array[(linear_array == 133)] = 3 linear_array[(linear_array == 165)] = 3 linear_array[(linear_array != 3)] = 0 # All non-Linear pixels set to 0 ### Combine arrays ### mspa_array = core_array + edge_array + linear_array # Save the final reclassified MSPA array save_array_to_tiff(mspa_array, output_mspa_tile_path, driver, mspa_ds) # Clean up mspa_ds = None # Report the time taken print(f"Output saved in: {SNH_tiles_path}") report_time(step_start_time, "Reclassifying MSPA Tiles in classes (Core, Edge, Linear), ") print("="*40) # A line of "=" for separation # Function to score extGL_path in Herbaceous Areal, and adding these into the reclassified MSPA tiles def SNH_tiles(recMSPA_tiles_path, extGL_path, SNH_tiles_path, start_feature=0): """ Merge tiles from recMSPA_tiles_path and extGL_path, applying a specific rule: Add the value of extGL_path only when recMSPA_tiles_path is not 1, 2, or 35. Save the merged tiles in SNH_tiles_path. """ step_start_time = time.time() print("\n" + "="*40) # A line of "=" for separation print(f"Adding Herbaceous Areal to MSPA...") # Get the list of tiles from recMSPA_tiles_path recMSPA_tiles = list(Path(recMSPA_tiles_path).glob("*.tif")) total_tiles = len(recMSPA_tiles) # Initialize GDAL driver once driver = gdal.GetDriverByName('GTiff') for i, recMSPA_tile in enumerate(recMSPA_tiles): if i < start_feature: continue # Skip tiles until reaching the start_feature index tile_id = recMSPA_tile.stem.split("_")[-1] extGL_tile = Path(extGL_path) / f"tile_{tile_id}.tif" # Check if corresponding extGL tile exists if not extGL_tile.exists(): print(f"Corresponding extGL tile not found for {tile_id}. Stopping process.") return # Stop the function if the tile is missing in extGL_path # Define output paths for merged tiles SNH_tile_path = Path(SNH_tiles_path) / f"tile_{tile_id}.tif" # Check if the output file already exists and is not corrupted if check_output_tif(str(SNH_tile_path)): continue # Skip if the file exists and is valid print(f"SNH classification. Processing tile {tile_id} [{i+1}/{total_tiles}]...") # Open recMSPA and extGL tiles recMSPA_ds = gdal.Open(str(recMSPA_tile)) extGL_ds = gdal.Open(str(extGL_tile)) recMSPA_array = recMSPA_ds.GetRasterBand(1).ReadAsArray() extGL_array = extGL_ds.GetRasterBand(1).ReadAsArray() ### Reclassify Herbaceous Areal (HA) ### extGL_array = extGL_array.astype(np.float32) extGL_array[(extGL_array == 1)] = 4 # Build the Semi Natual abitat layer. Merge logic -> add extGL value only when recMSPA value is not 1, 2, or 3 recMSPA_array_bin = np.where(recMSPA_array >0, 0, 1) SNH_array = (recMSPA_array_bin*extGL_array)+recMSPA_array # Save the final reclassified MSPA array save_array_to_tiff(SNH_array, SNH_tile_path, driver, recMSPA_ds) # Clean up recMSPA_ds, extGL_ds, recMSPA_array_bin,SNH_array = None, None, None, None print(f"Tile processed") # Report the time taken print(f"Output saved in: {recMSPA_tiles_path}") report_time(step_start_time, "Herbaceous Areal to MSPA Tiles") print("="*40) # A line of "=" for separation # Create Binary Grasslands layer def binary_CORINE(CORINE_tif, corine_bin_path, output_path, target_resolution=10): """ Reclassifies the CORINE raster to binary (sets values not in classes_to_mask to 0, and values in classes_to_mask to 1), saves it as 'corine_bin_tif', and then resamples this binary raster to the target resolution (default 50m) and saves the resampled result to 'corine_bin_path'. Parameters: - CORINE_tif (str): Path to the original CORINE raster file. - corine_bin_path (str): Path to save the resampled binary raster file. - output_path (str): Directory to save the binary raster before resampling. - target_resolution (int): The target resolution in meters for the resampled file (default is 50m). """ step_start_time = time.time() print("\n" + "="*40) print(f"Masking and resampling CORINE raster...") corine_bin_tif = os.path.join(output_path, f"U2018_CLC2018_{target_resolution}bin.tif") # Binary output file # Check if the resampled binary CORINE file already exists if os.path.exists(corine_bin_path): print(f"Resampled binary CORINE file already exists: {corine_bin_path}. Skipping reclassification and resampling...") return # Reclassify and save the binary version if not os.path.exists(corine_bin_tif): # Open the CORINE raster corine_ds = gdal.Open(CORINE_tif) corine_array = corine_ds.GetRasterBand(1).ReadAsArray() # Apply the mask: set values not in classes_to_mask to 0, and values in classes_to_mask to 1 classes_to_mask = [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22] # Classes to retain corine_mask_array = np.where(np.isin(corine_array, classes_to_mask), 1, 0) # Create the output binary raster (corine_bin_tif) driver = gdal.GetDriverByName('GTiff') output_masked = driver.Create(str(corine_bin_tif), corine_mask_array.shape[1], corine_mask_array.shape[0], 1, gdal.GDT_Byte, ["COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9", "BIGTIFF=YES"]) output_masked.SetGeoTransform(corine_ds.GetGeoTransform()) # Use the original GeoTransform output_masked.SetProjection(corine_ds.GetProjection()) # Use the original projection # Write the reclassified mask to the output GeoTIFF output_masked.GetRasterBand(1).WriteArray(corine_mask_array) # Save the reclassified array output_masked.GetRasterBand(1).SetNoDataValue(0) # Set NoData value output_masked.FlushCache() # Ensure everything is written to disk output_masked = None # Close the dataset print(f"Binary CORINE file saved as: {corine_bin_tif}") else: print(f"Binary CORINE file already exists: {corine_bin_tif}. Skipping binary creation...") # Step 2: Resample the binary version to the target resolution gdal.Warp(corine_bin_path, corine_bin_tif, xRes=target_resolution, yRes=target_resolution, resampleAlg="near", options=["COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9", "BIGTIFF=YES"]) print(f"Resampled CORINE file saved as: {corine_bin_path}") # Report the time taken report_time(step_start_time, "CORINE binary and resampling") print("="*40) # Merge the processed strips back together def merge_strips(input_dir, output_file): """ Merges all the processed strips (tiles) in the input_dir and saves the result as a single file to output_file. Parameters: input_dir (str): Directory containing the processed strips (tiles) to merge. output_file (str): Output file path where the merged raster will be saved. """ step_start_time = time.time() print("\n" + "="*40) # A line of "=" for separation print(f"Merging processed strips from {input_dir} \n into {output_file}...") # Collect all .tif files in the input directory strips_pattern = list(Path(input_dir).glob("*.tif")) # Convert PosixPath objects to strings and join them into a single space-separated string strips_pattern_str = ' '.join([str(tile) for tile in strips_pattern]) # GDAL merge command as a single string command_merge = f"gdal_merge.py -o {output_file} " \ f"-co COMPRESS=LZW " \ f"-co PREDICTOR=2 " \ f"-co BIGTIFF=YES " \ f"-co TILED=YES " \ f"-co BLOCKXSIZE=256 -co BLOCKYSIZE=256 " \ f"{strips_pattern_str}" # Run the merge command subprocess.run(command_merge, shell=True, check=True) print(f"Merged raster saved to {output_file}.") report_time(step_start_time, "Merging Strips") print("="*40) # A line of "=" for separation # Automatically create the compressed output file name by adding "_comp" before ".tif" output_file_comp = output_file.replace(".tif", "_comp.tif") # Apply additional compression with gdal_translate command_translate = f"gdal_translate -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 " \ f"-co ZLEVEL=9 -co BIGTIFF=YES -co TILED=YES -co BLOCKXSIZE=256 " \ f"-co BLOCKYSIZE=256 -a_nodata 0 {output_file} {output_file_comp}" # Run the gdal_translate command for further compression subprocess.run(command_translate, shell=True, check=True) print(f"Compressed raster saved to {output_file_comp}.") # Report the time taken report_time(step_start_time, "Merging Strips") print("="*40) # A line of "=" for separation # Mask SNH tiles with CORINE raster def mask_SNH(corine_bin_tif, SNH_tif, output_tif): """ Mask the values of SNH_tif where corine_bin_tif is not equal to 1, normalize to [0, 100], and save the result to output_tif (compressed). """ step_start_time = time.time() print("\n" + "="*40) print(f"Processing SNH with CORINE mask...") # Open SNH and CORINE binary files SNH_ds = gdal.Open(SNH_tif) corine_ds = gdal.Open(corine_bin_tif) SNH_array = SNH_ds.GetRasterBand(1).ReadAsArray() corine_array = corine_ds.GetRasterBand(1).ReadAsArray() # Mask the SNH array using the CORINE binary array masked_array = np.where(corine_array == 1, SNH_array, np.nan) # Mask SNH where CORINE != 1 # Save the result to output_tif with compression driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(output_tif, SNH_ds.RasterXSize, SNH_ds.RasterYSize, 1, gdal.GDT_Float32, options=["COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=9", "BIGTIFF=YES"]) out_ds.SetGeoTransform(SNH_ds.GetGeoTransform()) out_ds.SetProjection(SNH_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(masked_array) out_ds.GetRasterBand(1).SetNoDataValue(np.nan) out_ds.FlushCache() # Clean up SNH_ds, corine_ds, out_ds = None, None, None print(f"Masked SNH saved successfully to: {output_tif}") # Report the time taken report_time(step_start_time, "Masking and Normalizing SNH with CORINE binary") print("="*40) # Main function to run all steps def main(): print("Starting the process...") start_time = time.time() # 1. Create tiles for MSPA and Grasslands create_tiles(MSPA_tif, extGL_tif, MSPA_tiles_path, extGL_tiles_path, y2=2) # 2. Reclassify MSPA tiles into categories recla_MSPA(MSPA_tiles_path, recMSPA_tiles_path, start_feature=0) # 3. Combine MSPA and Grassland tiles into SNH tiles SNH_tiles(recMSPA_tiles_path, extGL_tiles_path, SNH_tiles_path, start_feature=0) # 4. Create binary CORINE raster binary_CORINE(CORINE_tif, corine_bin_tif, output_path, target_resolution=10) # 5. Merge SNH tiles back together into a single raster merge_strips(SNH_tiles_path, SNH_tif) # 6. Mask SNH raster using CORINE binary layer mask_SNH(corine_bin_tif, SNH_tif, SNHmasked_tif) # 7. Delete temporary folder shutil.rmtree(os.path.join(output_path,"Temp_recMSPAs"), ignore_errors=True) # Calculate and report total elapsed time print("="*40) # A line of "=" for separation report_time(start_time, "Script completed ") print("="*40) # A line of "=" for separation # Run the main function if __name__ == "__main__": main()