Tutorial: Spatial subset using Shapefile in snappy

Show Room
1

Snappy use a WKT file for spatial subset. I found everywhere but there was no option for read shapfile in snappy. Here I am sharing my experience how to convert shapefile to WKT file and subset image in snappy using polygon. Hope it will be helpful for you.

import snappy
from snappy import ProductIO
from snappy import GPF
from snappy import HashMap
from snappy import ProductUtils
from snappy import WKTReader
import gc
import os
from snappy import jpy
import geopandas as gpd


# read image and do all your preprocessing with image.
 img=ProductIO.readProduct(r'image with path')

 #Now we have to read the shapfile.

 shp=gpd.read_file(r'shapefile with path')     # here you have to read your polygon shapefile
 
shp     # output for my case.
 Out[111]: 
    Id                                           geometry
 0   0  POLYGON ((88.13664 22.50243, 88.12930 22.56486...

geom=str(shp['geometry'][0])     # get the geometry of polygon shapefile as string. 

# Here we create thw WKT file from the Shapefile.

geom = WKTReader().read(geom)

# now we can subset the image with the bounded area of shapefile.

par = HashMap()
par.put('copyMetadata', True)
par.put('geoRegion', geom)
im_sub = GPF.createProduct('Subset', par, image) #Spatial subset.

Here is the output of image Width and Height for my case:

Width:7730 Height:5388 # main corrected image
Width:1131 Height:1335 # after subset

This is how you can perform a spatial subset in snappy module using polygon shapefile.

@SubhadipDatta
Satellite Image Analyst

4 Likes
2

Thank you for sharing. I moved the topic to the show room.

3

Does the polygon need to be at a specific format for the geometry to be coordinates? If so what should it be and how do I convert it to the format WKTReader recognizes?

4

The reader reads the WKT as specified on WIkipedia.
Well-known text representation of geometry - Wikipedia

But, in SNAP Desktop, you can also import the shapefile into a product which covers the whole shapefile geometry and then select the imported geometry and export it as WKT.
https://step.esa.int/main/wp-content/help/?version=10.0.0&helpid=geometryWKT

5

There must be something more to it? I have a polygon with WKT geometry that appears to be using a different lat/lon format.
Example: “POLYGON ((263543.42676330917 1632506.1222286283, …”
and the WKT documentation does not seem tied to a specific coordinate system?
am I missing or misunderstanding something?

6

Ah, yes. Your coordinates look like UTM coords. WKT doesn’t care about it, but SNAP requires WGS84 Lat/Lon coordinates.

7

Great thanks for the clarification.

So to add to this tutorial…
Before using WKTReader on the geopandas dataframe, one would need to reproject the shpfile in geopandas to the proper coordinate system like so:

target_crs = "EPSG:4326"
shp = gpd.read_file(shp_file_path) #read shp file
shp = shp.to_crs(target_crs) #convert to target CRS
geom = str(shp['geometry'][0]) # take polygon text of shp file

The next question would then be how to get the target product’s CRS before applying the subset operator…