Subset of common area of Sentinel-1 GRD images

development python
#1

Hello everyone and thanks for all the insights I found in the forum so far.

I’m writing a python script that should extract the overlapping area of 2 (or more) Sentinel-1 GRD images and save each cropped image as a GeoTIFF.
I was able to calibrate the products and apply terrain correction. I can also save products in the .tif format (despite the lack of official documentation, snappy is quite easy and fun to work with, once you get the basics down).
The only thing I’m missing at this point is subsetting each image of the common area, which I would ideally perform before saving my images, in order to make things faster.

I’m pretty new to SAR processing, so every tip will be helpful!

#2

So, you don’t have a predefined area which you want to use for subset but you want to find the overlapping area and then subset to this area?
If so, you need to get the boundary first of all products and then create the intersection to find the common overlapping area.The area can be converted to a WKT string and used with the Subset operator with each product.
You can use this method to get the boundary:
org.esa.snap.core.util.ProductUtils#createGeoBoundaryPaths(org.esa.snap.core.datamodel.Product)

The resulting path (you probably only need the first path in the array) needs to be converted into a geometry
The following is a Java snippet which can be converted into python
GeometryFactory factory = new GeometryFactory()
Geometry boundary = JTS.shapeToGeometry(geoBoundaryPath, factory);
if (boundary instanceof LinearRing) {
boundary = factory.createPolygon((LinearRing) boundary, null);
}
return boundary;

If you have two or more geometries you can create the intersection:

intersection = geometry1.intersection(geometry2)

The intersection can be used as parameter “geoRegion” with the Subset operator.

You need to import some types in your script e.g.

jpy.get_type('com.vividsolutions.jts.geom.GeometryFactory')

or
jpy.get_type(‘org.geotools.geometry.jts.JTS’)

1 Like
#3

@marpet thank you so much!

I was able to do it by following your suggestion, except for the condition

if (boundary instanceof LinearRing)

which I didn’t know how to convert to python, but it worked anyway. Is it really necessary?
Thank you again.

Edit: correct line in quote

#4

If it works it is fine.

I now noticed that my sample code is a mixture of Java and Python anyway.

To use the GeometryFactory you need to do (hope now this is real Python):

JTS = jpy.get_type(‘org.geotools.geometry.jts.JTS’)
GeoFactory = jpy.get_type('com.vividsolutions.jts.geom.GeometryFactory')
gf = GeoFactory()
boundary = JTS.shapeToGeometry(geoBoundaryPath, gf)
1 Like
#5

No problem, I managed to do it by looking at similiar imports like the one for HashMap. By the way, I just noticed I copied the wrong line in the quote in my previous reply, but again, it works fine.
Thank you