[Python] Efficiently Rasterizing a Shapefile into Slant Range with Sentinel-1 SAR Data

I'm currently working with Sentinel-1 SAR data and need to rasterize a shapefile into the slant range. *

To achieve this, I created a GeoTIFF image with two bands (latitude and longitude) using the SNAP (ESA) software. The BandMaths tool was utilized to calculate these coordinates, which are derived from the precise orbit of the satellite available for download at step.esa.int. The challenge I'm facing is that my GeoTIFF is not associated with any Coordinate Reference System (CRS) because it remains in slant range rather than being projected into ground range.

I implemented a nested loop to extract information related to specific points (i, j) in my GeoTIFF matrix using GeoPandas, but this approach is quite slow. To improve performance, I attempted to parallelize the operation with Python's multiprocessing module, but it still does not perform as efficiently as I would like.

Here is my code :

# Library importations
from osgeo import gdal
import matplotlib.pyplot as plt
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
from multiprocessing import Pool, Array

# Read raster (Slant Range)
raster_lat_lon = gdal.Open('subset_lat_lon.tif', gdal.GA_ReadOnly)
lat, lon = raster_lat_lon.GetRasterBand(1), raster_lat_lon.GetRasterBand(2)
arr_lat, arr_lon = lat.ReadAsArray(), lon.ReadAsArray()
arr_lat_line, arr_lat_col = arr_lat.shape

# Read shapefile (Ground Range)
shapefile = gpd.read_file('tbe.shp')
shapefile_yearly = shapefile[shapefile['ANNEE'] == 2020]

# Get pixel value

def get_pixel_value(lin_col):
i, j = lin_col
# Create a point with shapely
pixel_point = Point(arr_lon[i, j], arr_lat[i, j])
# Test if the point is in the shapefile
sbw = shapefile_year[shapefile_year.geometry.contains(pixel_point)]

# If a line in the shp is associated, put it in the empty array
if len(sbw) != 0:
value = sbw['Ia'].iloc[0]
empty_array[i * arr_lat_col + j] = value # Put the value into the 1D array


# Set empty array to fill + set args
empty_array = Array('d', arr_lat_line * arr_lat_col)
args = [(i,j) for i in range(arr_lat_line) for j in range(arr_lat_col)]

# Use multiprocessing to parallelise
with Pool(50) as p:
p.map(get_pixel_value, args)

# Convert the array into NumPy array to visualise it
empty_array_np = np.frombuffer(empty_array.get_obj()).reshape((arr_lat_line, arr_lat_col))

# Visualisation of the "Slant Range" shapefile rasterized
plt.imshow(arr_vide_np, cmap='gray')
plt.title('Spruce budworm 2020 rasterized, slant range')
plt.show()



I was expecting to find a more efficient method or tool to optimize the rasterization of the shapefile into slant range data. By using only "ground range" data, there are tools available on QGIS but I wasn't able to find anything that can help me rasterize it without having to use the ground range. I would appreciate any suggestions or insights on how to optimize this rasterization process further.

Continue reading...