From d607dc2fd8febd919e1ceed1b9bd61d796ac4834 Mon Sep 17 00:00:00 2001
From: Emma Oceguera-Conchas <113361856+E-O-Conchas@users.noreply.github.com>
Date: Fri, 3 Nov 2023 10:00:31 +0100
Subject: [PATCH] Added comprehensive Python script for rasterization of
shapefile data from PostgreSQL
---
1.rasterize_mask_and_fragmentation_5m.py | 165 +++++++++++++++++++++++
1 file changed, 165 insertions(+)
create mode 100644 1.rasterize_mask_and_fragmentation_5m.py
diff --git a/1.rasterize_mask_and_fragmentation_5m.py b/1.rasterize_mask_and_fragmentation_5m.py
new file mode 100644
index 0000000..42dcaa6
--- /dev/null
+++ b/1.rasterize_mask_and_fragmentation_5m.py
@@ -0,0 +1,165 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 19 10:05:00 2023
+
+@author: no67wuwu
+
+work:
+This script is used to rasterize a shapefile that has been imported to a PostGIS database.
+It produces two output rasters: a mask and a fragmentation raster.
+
+"""
+
+# Import necessary modules
+import sys
+sys.path.insert(0, r"I:\biocon\Emmanuel_Oceguera\projects\2023_03_NaturaConnect\tools")
+import MacPyver as mp
+
+
+
+# Import additional libraries
+import os
+import subprocess
+from multiprocessing import Pool
+
+
+def get_stdout(cmd, verbose=False):
+
+ """
+ Send a command to the command line and fetch the return.
+
+ :param cmd: The command to be executed
+ :param verbose: Whether to print the command's output
+ :return: The command's standard output
+
+ from http://blog.kagesenshi.org/2008/02/teeing-python-subprocesspopen-output.html
+ """
+ p = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
+ stdout = []
+ while True:
+ line = p.stdout.readline()
+ line = line.decode('utf-8') # Decode the bytes into a string
+ line = line.replace('\n','').replace('\r','')
+ if line != '':
+ stdout.append(line)
+ if verbose:
+ print(line),
+ if line == '' and p.poll() != None: #p.poll is checking if the process is still running
+ #is it is running it returns None
+ break
+ return stdout
+
+
+# Define a function to execute a GDAL command and log the result
+def execute_gdal_command(cmd, description, verbose=False):
+ """
+ Execute a GDAL command and log the result.
+
+ :param cmd: The GDAL command to be executed
+ :param description: A description of the command
+ :param verbose: Whether to print the command's output
+ :return: The command's standard output
+ """
+ print(f"Executing {description}...")
+ result = get_stdout(cmd, verbose)
+ print(f"{description} completed.")
+ return result
+
+
+# Define the main function
+def main(sbbox):
+ """
+ Main function to rasterize a sub-bounding box.
+
+ :param sbbox: The sub-bounding box coordinates (xmin, ymin, xmax, ymax)
+ """
+ # Extract the extent coordinates
+ xmin, ymin, xmax, ymax = sbbox
+
+ # Define output file names for mask and fragmentation rasters
+ outname_mask = outpath_mask + '\\mask_tile_x{0}_y{1}_1.tif'.format(int(xmin/10000), int(ymin/10000))
+ outname_frag = outpath_frag + '\\frag_tile_x{0}_y{1}_1.tif'.format(int(xmin/10000), int(ymin/10000))
+
+ # Define GDAL command for fragmentation rasterization
+ cmd_frag = """gdal_rasterize -tr 5 5 -co "COMPRESS=DEFLATE" -burn 9 -ot Byte -te {0} {1} {2} {3} -a_srs EPSG:3035 \
+ -sql "select st_intersects(f.geom, st_geomfromtext('POLYGON (({0} {1}, {2} {1}, {2} {3}, {0} {3}, {0} {1}))', 3035)), \
+ f.frag_code, f.id from fragmentation.forest_clc2018_v2020_clip_ms f where \
+ st_intersects(f.geom, st_geomfromtext('POLYGON (({0} {1}, {2} {1}, {2} {3}, {0} {3}, {0} {1}))', 3035)) \
+ and f.frag_code = 999"\
+ PG:"host=localhost user=postgres dbname=NaturaConnect-Connectivity password=07089452" {4}""".format(xmin, ymin, xmax, ymax, outname_frag)
+
+ # Define GDAL command for mask rasterization
+ cmd_mask = """gdal_rasterize -tr 5 5 -co "COMPRESS=DEFLATE" -burn 1 -ot Byte -te {0} {1} {2} {3} -a_srs EPSG:3035\
+ -sql "select st_intersection(f.geom, st_geomfromtext('POLYGON (({0} {1}, {2} {1}, {2} {3}, {0} {3}, {0} {1}))',3035)), \
+ f.id from NUTS2021_V22_EU_MS.NUTS2021_v22_EU_MS_3035 f where \
+ st_intersects(f.geom, st_geomfromtext('POLYGON (({0} {1}, {2} {1}, {2} {3}, {0} {3}, {0} {1}))', 3035))" \
+ PG:"host=localhost user=postgres dbname=NaturaConnect-Connectivity password=07089452" {4}""".format(xmin, ymin, xmax, ymax, outname_mask)
+
+ # Execute GDAL commands for rasterization
+ execute_gdal_command(cmd_frag, "Fragmentation Rasterization")
+ execute_gdal_command(cmd_mask, "Mask Rasterization")
+
+
+# Define the path to the input shapefile and output paths for mask and fragmentation rasters
+inshape = r"I:\biocon\Emmanuel_Oceguera\projects\2023_09_Fragmentation\data\eu_reference_grid\europe_100km.shp"
+outpath_frag = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\frag_raster"
+outpath_mask = r"S:\Emmanuel_OcegueraConchas\eu_fragmentation_forest\land_raster"
+
+# Define a command to get the total extent of the shapefile
+cmd = 'ogrinfo -al -so ' + inshape
+
+
+# Get the total extent of the shapefile
+bbox = [x for x in get_stdout(cmd) if x.startswith('Extent')][0]
+bbox = bbox.replace('Extent: ','').replace('(','').replace(')','').split(' - ')
+bbox = [float(z) for x in bbox for z in x.split(',')]
+
+x_range = bbox[2] - bbox[0]
+y_range = bbox[3] - bbox[1]
+
+xmin = bbox[0]
+ymin = bbox[1]
+
+
+
+# Create a list of sub-bounding boxes covering the entire extent
+sub_bbox_L = []
+for x in range(0, int(x_range/100000 + .5)):
+ for y in range(0, int(y_range/100000 + .5)):
+ sub_bbox_L.append([xmin + 100000 * x, ymin + 100000 * y, xmin + 100000 * (x + 1), ymin + 100000 * (y + 1) ])
+
+
+# Run the main function over all sub-bounding boxes using multiprocessing
+if __name__ == "__main__":
+ pool = Pool(10)
+ for foo in pool.imap_unordered(main, sub_bbox_L):
+ bar = 1
+
+
+
+
+#%%
+# gdal_rasterize [--help] [--help-general]
+# [-b <band>]... [-i] [-at]
+# [-oo <NAME>=<VALUE>]...
+# {[-burn <value>]... | [-a <attribute_name>] | [-3d]} [-add]
+# [-l <layername>]... [-where <expression>] [-sql <select_statement>|@<filename>]
+# [-dialect <dialect>] [-of <format>] [-a_srs <srs_def>] [-to <NAME>=<VALUE>]...
+# [-co <NAME>=<VALUE>]... [-a_nodata <value>] [-init <value>]...
+# [-te <xmin> <ymin> <xmax> <ymax>] [-tr <xres> <yres>] [-tap] [-ts <width> <height>]
+# [-ot {Byte/Int8/Int16/UInt16/UInt32/Int32/UInt64/Int64/Float32/Float64/
+# CInt16/CInt32/CFloat32/CFloat64}] [-optim {AUTO|VECTOR|RASTER}] [-q]
+# <src_datasource> <dst_filename>
+#%%
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
--
GitLab