SubsetOp with .shp

development python
1

Hello,

i want to subset a Sentinel 2 image with a ESRI shapefile with snappy (with many .shp). Is this possible with the SubsetOp or is a WKT string necessary for that. is it in general possible to create a subset with a shp. via snappy? Which Operators can be used tu fullfill this task.

2

Currently a shapefile is not supported, but we have it on the agenda (SNAP-369)

Maybe the links provided in this post can help:

3

I just want to point out that geoPackage is slowly becoming the standard geospatial data format in some GIS softwares, instead of shapefiles.

4

How can I convert a .shp to a wkt if possible with snappy.

5

You can try the following approach for only a few shapefiles:

or use an online conversion tool like: SHP to WKT Converter Online - MyGeodata Cloud

For snappy there is no good API yet. It is better to convert it upfront.

6

i am looking at this code and dont get my head arround this part:

if len(sys.argv) != 3:
print(“usage: %s " % sys.argv[0])
print(” %s ./TEST.N1 “POLYGON((15.786082 45.30223, 11.798364 46.118263, 10.878688 43.61961, 14.722727”
“42.85818, 15.786082 45.30223))”" % sys.argv[0])
sys.exit(1)

file = sys.argv[1]
wkt = sys.argv[2]

Why is this part necessary. Were can i insert my wkt file?

7

This is only used to check that 3 arguments are specified with the script call. If not a usage example is printed.

The script does not read from a WKT file it uses a provided WKT string.
But you can simply read the string from a text file with the usual Python API.

10

Sorry for the monolog :slight_smile:

The following code is running know:
import sys
import snappy
from snappy import ProductIO, WKTReader

SubsetOp = snappy.jpy.get_type(‘org.esa.snap.core.gpf.common.SubsetOp’)

wkt = “MULTIPOLYGON (((417085.99999912607 5852152.0000229785, 416986.74999912473 5852142.0000229785, 416916.9999991239 5852141.000022978, 417065.7499991256 5852260.500022979, 416931.24999912427 5852351.00002298, 416765.49999912234 5852225.500022979, 416746.99999912194 5852237.0000229785, 416887.7499991235 5852350.000022979, 416812.9999991228 5852395.000022979, 416795.4999991227 5852395.50002298, 416772.9999991223 5852389.5000229785, 416752.2499991221 5852377.500022979, 416722.2499991217 5852352.500022979, 416673.2499991213 5852412.50002298, 416666.49999912106 5852420.5000229785, 416775.2499991223 5852513.00002298, 416784.74999912234 5852506.500022981, 416948.4999991245 5852396.00002298, 416954.7499991245 5852392.00002298, 417103.99999912584 5852291.000022981, 417270.24999912805 5852179.0000229785, 417274.99999912817 5852173.000022979, 417225.49999912747 5852167.000022979, 417191.52739912714 5852163.347022979, 417085.99999912607 5852152.0000229785)))”
geom = WKTReader().read (wkt)
print(“Reading…”)
product = ProductIO.readProduct(‘C:\Users\Savi\image2.dim’)

op = SubsetOp()
op.setSourceProduct(product)
op.setGeoRegion(geom)

sub_product = op.getTargetProduct()

print(“Writing…”)
ProductIO.writeProduct(sub_product, ‘C:\Users\Savi\subset’, “BEAM-DIMAP”)

print(“Done.”)

But i recieve a product without bands amd following error message:

C:\Python27\python.exe C:/Users/Niels/.PyCharmCE2018.1/config/scratches/SubsetOp.py
INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters
Reading…
INFO: org.hsqldb.persist.Logger: dataFileCache open start
WARNING: org.esa.snap.core.gpf.common.SubsetOp: No intersection with source product boundary image2
Writing…
Done.

Process finished with exit code 0

11

Your region does not intersect with the product. That’s why it is empty

WARNING: org.esa.snap.core.gpf.common.SubsetOp: No intersection with source product boundary image2

You need to specify lat/lon coordinates. UTM-coordinates do not work.

12

To subset with .shp you can convert shapefile to wkt using osgeo: here is an example i used to cut and it worked. i do hope it is useful for you.

def GetSubset(productst):

# Convert shapefile to Wkt
drv = ogr.GetDriverByName("ESRI Shapefile")

inSHP = drv.Open(shp, 1)
shpLayer = inSHP.GetLayer()
shpLayer.ResetReading() 

for element in xrange(shpLayer.GetFeatureCount()):
    geometry = shpLayer.GetFeature(element)
    wkt = geometry.GetGeometryRef().ExportToWkt()

subset_list = []
subset_folder = os.path.join(home, 'Subset')
CreateFolder(subset_folder)

WKTReader = snappy.jpy.get_type('com.vividsolutions.jts.io.WKTReader')
geom = WKTReader().read(wkt)
for product in products:
    subset_name    = product .getName() 
    print subset_name
    subset_path     = os.path.join(subset_folder, subset_name)
    if not os.path.isfile(subset_path):
        subset_hashmap = HashMap()
        subset_hashmap.put('copyMetadata', True)
        subset_hashmap.put('geoRegion', geom)
        subset_hashmap.put('outputImageScaleInDb', False)
        subset = GPF.createProduct('Subset', subset_hashmap, product )
        subset_out = subset_folder + subset.getName()
        subset_list.append(subset)
        ProductIO.writeProduct(subset, subset_out, 'BEAM-DIMAP')
print 'SUBSET PRODUCT..............................DONE'
return(subset_list)
2 Likes
13

did you setup qgis.core and pyqgis for youre python interpreter were you run snappy?
If yes how you did that? Would be really nice to get to know youre solution

14

i did not use QGIS, you just need install gdal, osgeo and then import them

from osgeo import gdal,ogr, osr
gdal.UseExceptions()
ogr.UseExceptions()
osr.UseExceptions()
os.environ['GDAL_DATA'] = 'C:\Python27\Lib\site-packages\osgeo\data\gdal'
16

Hello again,

i tried out youre code. I keep geting the sam error like before:
5860069.5,439760.75 5860069.0,439741.75 5860068.5,439668.5 5860067.0,439659.75 5860088.5,439600.0 5860232.0,439553.0 5860344.0,439504.5 5860459.5,439481.5 5860514.5,439434.75 5860625.9999,439388.25 5860739.0,439366.75 5860791.5,439284.25 5860992.0,439277.0 5861009.0,439243.328836564 5861088.24611345,439374.477662009 5861110.11888024,439659.968037906 5861157.73245506,439769.5 5861176.0,439923.25 5860804.5,439965.26061362 5860701.69668902,439996.78477457 5860624.55456131,440034.240594783 5860532.89719186,440083.162376539 5860413.18170309,440144.904511087 5860262.09379192,440158.862363297 5860227.93781805,440157.75 5860227.5,440158.5396 5860225.244)))
WARNING: org.esa.snap.core.gpf.common.SubsetOp: No intersection with source product boundary after_reprojected_reprojected

Process finished with exit code 0

I put a print(wkt) in the code to see the wkt. what is in front of the error is a part of my polygon. I think that is lat/lon. I reprojected the product to the same reference system that my .shp has. Any Idea whats the mistake?

17

It is still the same problem:

Your coordinates are in meter not in geographical degrees. Change to Lat/Lon then it should work.

18

I searched the net but dont found a solution. How would you do that? Is there a operator for snappy?

19

You need to change your coordinates system from meter to degrees using QGIS. (convert to WGS84)

1 Like
20

I finally solved my problem with the land sea mask operator:

import vector to b4 after

file1 = (‘C:/Users/Niels/Extractby_mask/wald3.shp’)

separateShapes = False

GPF.getDefaultInstance().getOperatorSpiRegistry().loadOperatorSpis()
HashMap = jpy.get_type(‘java.util.HashMap’)
parameters = HashMap()
parameters.put(‘vectorFile’, file1)
parameters.put(‘separateShapes’, separateShapes)

import_vector_before = GPF.createProduct(‘Import-Vector’, parameters, B4_before)

#subset b4_after with impoted vector

parameters.put(‘geometry’, ‘wald3’)
parameters.put(‘invertGeometry’, False)
parameters.put(‘byPass’, False)

subset_b4_after = GPF.createProduct(‘Land-Sea-Mask’, parameters, import_vector_after)

maybe that is usefull for someone.