I would like to know the process of visualization of a Sentinel 1 scene post its calibration and terrain correction in Snappy. My basic concern is to obtain the image matrix which basically gives me the pixel values and latitude & longitude associated with each pixel. This matrix is crucial to me since I will be writing my codes based on the extent of the matrix (rows & columns). I would like to have some help on this matter.
By visualization I was actually referring to how we display post loading it on Python. I have the following script to load it and perform a few basic operations like calibration and terrain correction. What I fail to understand is how the data is being stored. Being a Matlab user I knew the image is a matrix of MxN size and I was able to display and image and choose my regions of interest by merely looking at the Row-Column values. I need some clarity on this.
#####################################################################################
import snappy
from snappy import ProductIO
from snappy import HashMap
from osgeo import gdal
import matplotlib.pyplot as plt
import os, gc
from snappy import GPF
GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
HashMap = snappy.jpy.get_type('java.util.HashMap')
# Loop through all Sentinel-1 data sub folders that are located within the path
path = "C:\\Users\\Archeron Defense\\Desktop\\Python_Journey\\Opening_Sentinel_1_Data\\"
# S1A_IW_GRDH_1SDV_20180301T020654_20180301T020719_020815_023B0D_E9B3.SAFE\\
for folder in os.listdir(path):
gc.enable()
output = path + folder + "\\"
timestamp = folder.split("_")[:4] #takes the chunk from file name till the forth element. Each element is spaced by a "_"
# This enable to get S!A_IW_GRDH_1SDV from the
date = folder.split("_")[4:5] # Crops the date & time portion from the file name (ex. for file name S1A_IW_GRDH_1SDV_20180301T020654_20180301T020719_020815_023B0D_E9B3.SAFE
# it crops 20180301T020654 and puts it in a list ['20180301T020654'])
List_to_str_date = ''.join(date) # Convert the list to string to be able to extract just the date by counting characters in the string
date = List_to_str_date[0:8]
# Read Sentinel-1 data product:
# S1A_IW_GRDH_1SDV_20180301T020654_20180301T020719_020815_023B0D_E9B3.SAFE
sentinel_1 = ProductIO.readProduct(output + "\\manifest.safe")
print(sentinel_1)
# If polarization bands are available, spolit up code to process VH and VV intensity data separately.
# The first step is the calibration procedure by transforming the DN values to Sigma Nought respectively.
# Parameters to output the image in Decibels is also possible.
pols = ['VH','VV']
for p in pols:
polarization = p
### CALIBRATION
parameters = HashMap()
parameters.put('outputSigmaBand', True)
parameters.put('sourceBands', 'Intensity_' + polarization)
parameters.put('selectedPolarisations', polarization)
parameters.put('outputImageScaleInDb', False)
calib = output + '_'.join(timestamp) + '_' + date + "_calibrate_" + polarization
target_0 = GPF.createProduct("Calibration", parameters, sentinel_1)
ProductIO.writeProduct(target_0, calib, 'BEAM-DIMAP')
# Apply Range Doppler Terrain Correction to correct for layover and foreshortening effects, by using the SRTM 3 arcsecond product (90m) that is downloaded automatically.
# You could also specify an own DEM product with a higher spatial resolution from a local path:
### TERRAIN CORRECTION
calibration = ProductIO.readProduct(calib + ".dim")
parameters = HashMap()
parameters.put('demResamplingMethod', 'NEAREST_NEIGHBOUR')
parameters.put('imgResamplingMethod', 'NEAREST_NEIGHBOUR')
parameters.put('demName', 'SRTM 3Sec')
parameters.put('pixelSpacingInMeter', 10.0)
parameters.put('sourceBands', 'Sigma0_' + polarization)
terrain = output + '_'.join(timestamp) + '_' + date + "_corrected_" + polarization
target_2 = GPF.createProduct("Terrain-Correction", parameters, calibration)
#ProductIO.writeProduct(target_2, terrain, 'GeoTIFF')
ProductIO.writeProduct(target_2, terrain, 'BEAM-DIMAP')
## Image Visualization
# ds = gdal.Open('S1A_IW_GRDH_1SDV_corrected_VH.tif')
# plt.imshow(None)
#################################################################################
The fundamental issue with using Snappy is proper documentation. The preliminary stage of learning is hampered by the lack of documentation for different modules within a package. Can I get some help on where will I be able to find the detailed documentation.
There is documentation missing for snappy. That’s true.
As snappy is just a bridge into the Java API it is good to know the Java API. See SNAP API.
We have ideas on how to improve this situation but the realisation needs some more time.
Beside the guides I mentioned already above, there examples linked in the tutorial page. See External Resources section.
Another way of finding examples is to search this forum for GPF.createProduct
That’s the right way to store the data. You can also choose between different output formats. You have tried already GeoTIFF.
As a hint, you don’t need to write the data of the intermediate results. You can directly provide the result of one operator as input to the next operator.
Regarding the visualisation. There is a recent thread about it.
The actual data of band can be retrieved by:
import numpy as np
band = product.getBand(Sigma0_VV')
w = band.getRasterWidth()
h = band.getRasterHeight()
data = np.zeros(w*h, dtype=np.float32)
# this is the full image
band.readPixels(0, 0, w, h, data)
# Now you have the data in the numpy array 'data' and you can use it with e.g. matplotlib for visualisation