Processing S1 with SNAP python

Hello everyone!
I want to create a composite image in RGB from the processed Sentinel-1 using snappy.
I do calibration, filtering, and dB conversion. The result is a product with Brands = [Sigma 0_HH_db, Sigma 0_HV_db].
Please tell me how I can make an RGB image from this, similar to manual processing in SNAP: Open RGB Image Window --> Dual Pol Difference Sigma0 db HH + HV ?
Thank you very much in advance!

No idea about Snappy since it is not really reliable… I’m personnally doing this thanks to the rasterio and numpy python library: reading VV and VH bands, applying bands operations, and writing a 3band raster.

Would like to share my function, but it is not yet strong enough, sorry…

1 Like

Thank you, @AriJeannin !
I’ll try to use rasterio for this task, but I am afraid that this option will be very slow due to the large size of the images…

Have a look at Snapista or at pyroSAR.

They are python wrappers around gpt that may help you… you can try to join things via the SNAP operators or you may simply pre-process data with those and end-up with GeoTIFFs that you can further process with other Python frameworks.

1 Like

@cristianoLopes , thank you very much!
I hope it really can help me!
I’ll leave our solution here when it’s ready)))

@AriJeannin , thank you very much again!
I tried different options for creating RGB-images, including the one suggested by @marpet with snappy ImageManager.getInstance().createColoredBandImage:
Save image from Sigma0_HH_db band with snappy

And the best and most convenient result was obtained with using rasterio. It works really very fast using my functions with numba - GPU accelerator! Raster data is read very quickly and I can use my own functions to implement various image processing.

After processing I normalize each of the channels using the OpenCV library, create a numpy array and write each of the channels to it. The process of creating an image is much faster than using snappy.

Thank all of you so much for your help!

snappy is slow… hence the appearance of the different gpt wrappers.

if I understand this correctly - you processed to sigma0 with SNAP then got the outputs into rasterio for merging and opencv for cleanup?
Can you share how much faster this is in comparison with a pure SNAP graph (assuming this is even possible)?

@CDOEspiao, Glad to hear that!

Thanks for GPU accelerator tips. If you have a well running function for generating RGB, feel free to share!

I’m currently working on generating RGB thanks to EOBrowser scripts. Here some inspiration:


1 Like

@cristianoLopes , yes, that’s right. After processing the SAR image with snappy, I save it in Geotiff format (I will need this for future use):

ProductIO.writeProduct(output, dB_file, 'GeoTIFF')

And then I open it with rasterio and work as with a numpy array:

dataset = rasterio.open(file_path)
band1, band2 = dataset.read(1), dataset.read(2)

And then you can work with these arrays using the OpenCV, scikit-image, PIL, and so on libraries.
For example:

width, height = dataset.width, dataset.height
image = np.zeros((height, width), np.float32)
cv2.normalize(band1, image, 0, 255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

Unfortunately, I can’t say for sure whether it’s faster than other methods, because I haven’t done enough experiments yet… I hope that I will have some new results soon)

@AriJeannin , Wow! This is really impressive!
Thank you for the links, these are very useful sources, and it’s really cool!

To save RGB images, I use something like the following:

import rasterio
import cv2
import numpy as np
from PIL import Image

dataset = rasterio.open(file_path)
width, height = dataset.width, dataset.height

band1 = dataset.read(1)
band2 = dataset.read(2)
band3 = abs(band1 - band2)

red_norm_pixels = np.zeros((height, width),np.float32)
cv2.normalize(band1, red_norm_pixels, 0, 255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

green_norm_pixels = np.zeros((height, width),np.float32)
cv2.normalize(band2, green_norm_pixels, 0, 255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

blue_norm_pixels = np.zeros((height, width),np.float32)
cv2.normalize(band3, blue_norm_pixels, 0, 255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

new_image = np.zeros((height, width, 3))
new_image[:,:,0] = red_norm_pixels
new_image[:,:,1] = green_norm_pixels
new_image[:,:,2] = blue_norm_pixels

snimok = Image.fromarray(np.uint8(new_image), mode='RGB')
snimok.save("rgb.png")

And before saving the image, you can do any operations to process this image, depending on your needs.You can also do any other transformations as band 3, again, depending on what you want to highlight in the image.

1 Like

Great question. However, Seems the right solution is using numpy.array.
The main difference, I think, is depending on the way of processing Images.
The type of Image, for example, if want to use 6 bands layer for CNNs, I think the only way is GTiff.

I use gdal to write geocode information in each image, Maybe The speed is not that important cause I don`t have much experience too, haha.

It goes like this, contact me If you need more information

Geotiff_driver = gdal.GetDriverByName(‘GTiff’)
out_img = Geotiff_driver.Create(file_path, XSize, YSize, bands_number, gdal.GDT_Byte)
Geotiff_driver.SetProjection(xxxxxxxxxxxxx)
Geotiff_driver.SetGeoTransform(xxxxxxxxxxxx)


calculation

Out_img = Geotiff_driver.GetRasterBand(1)
Out_img.WriteArray(array1)

1 Like

@Ambrun , thank you very much for your idea!
I will try this option and let you know what results I will get!