Welcome to the BEAM Forum!

We encourage you to sign in our forum and participate in the BEAM community. The forum is maintained by the BEAM project team who will most likely answer your questions within 24 hours (except during common holidays) - if not done by other community members. Collaborate, share your knowledge and learn from other users!

If you don't find what you are looking for, please also consider the following external forums:

Combination View Flat View Tree View
Threads [ Previous | Next ]
clip using .shp file
toggle
clip using .shp file
1/8/16 4:05 PM
hi, everyone
i want to clip the product with a .shp file in my java code,
but i have no idea about this, help, please!

any idea will be deeply appreciated.
thanks in advance
Attachments: LC81320352014198LGN00_sr_band3_band_1.bmp (219.8k)
Flag Flag
RE: clip using .shp file
1/11/16 3:19 PM as a reply to zhu yu.
I haven' really tried the following code but it should something like this:
 1File file = new File("path/toShapefile.shp");
 2Product product = new Product("MyProduct", "type", 500, 500);
 3       
 4FeatureCollection<SimpleFeatureType, SimpleFeature> featureCollection = FeatureUtils.loadFeatureCollectionFromShapefile(file);
 5
 6Geometry productBounds = FeatureUtils.createGeoBoundaryPolygon(product);
 7FeatureUtils.clipCollection(featureCollection, DefaultGeographicCRS.WGS84,
 8                                    productBounds, DefaultGeographicCRS.WGS84,
 9                                    null, null, new NullProgressMonitor());
10       
11ReferencedEnvelope bounds = featureCollection.getBounds();
12Rectangle2D rect = new Envelope2D(bounds);
13
14HashMap<String, Object> parameters = new HashMap<>();
15parameters.put("region", rect);
16Product subset = GPF.createProduct("Subset", parameters, product);


You need to retrieve the source product and the shapefile from somewhere else.

Hope it helps.

Marco
Flag Flag
RE: clip using .shp file
1/12/16 12:10 PM as a reply to Marco Peters.
Thank you very much, Marco!

I've tried your code, but an exception was thrown:

org.esa.beam.framework.gpf.OperatorException:
Unknown operator 'Subset'. Note that operator aliases are case sensitive.
at org.esa.beam.framework.gpf.GPF.createProduct(GPF.java:222)
......

I didn't know how to solve it, so i tried to replace the 14th-16th lines with the following code:
1
2ProductSubsetDef subsetdef=new ProductSubsetDef();
3subsetdef.setRegion(rect.getBounds());
4Product subset=product.createSubset(subsetdef, null, null);


unfortunately, another kind of exception was thrown:

java.lang.IllegalArgumentException: [value] is [134290142712] but should be in the range [0] to [4294967295]
at org.esa.beam.util.Guardian.assertWithinRange(Guardian.java:298)
at org.esa.beam.dataio.geotiff.internal.TiffValueRangeChecker.checkValue(TiffValueRangeChecker.java:48)
......

question 1:
In fact, the polygon(in shapefile) has a large range, which exceeds the boundary of this product. whether it is the reason to throw exception?
question 2:
I only want to obtain the data in the irregular polygon, not in a boundary rectangle, whether the code you post can achieve this?

best wishes!
Zhuyu
Flag Flag
RE: clip using .shp file
1/13/16 10:19 AM as a reply to zhu yu.
I think what is missing is the initialisation of the OperatorRegistry. The following ensures that the Operators are loaded and can be found.
It needs to be executed only once at the beginning of your code.
1registry = GPF.getDefaultInstance().getOperatorSpiRegistry()
2registry.loadOperatorSpis()

However, the way you have implemented it should also work.

The size of the shapefile shouldn't be a problem. During the code I've written in my first post the shapefile is clipped to the product bounds.
Probably the product it self is to big for storing it in GeoTiff. Very likely for Landsat8.
I think you have 3 options.
  1. Subsample the data.
    You can specify a subsample in X-/Y-direction (subsetdef.setSubSampling(int subSamplingX, int subSamplingY))
  2. Select specific bands.
    You can select only the bands you really need for your application. (subsetdef.addNodeNames(String[] names))
  3. Use different storage format.
    GeoTiff can only store up to 2GB of data. So you need to use an other format. You can try 'BigGeoTiff' or 'NetCDF-CF' or 'NetCDF4-CF'. You can use these names as format names.


An other option is to set an environment property.
Set landsat.reader.donotscaletopanresolution to true. Afterwards the data will not be scaled up to the resolution of the panchromatic band (Band 8) during the reading process. You can either set it as parameter when you run you program with the -D java option or set it in your code via System.setProperty("landsat.reader.donotscaletopanresolution", "true").
Flag Flag
RE: clip using .shp file
1/13/16 3:56 PM as a reply to Marco Peters.
thanks again, Marco. it is under your instructive idea i solved my problem.

1) as to regular clip, the correct code is:
 1
 2File shapefile = new File("path/filename.shp");
 3String productfile = "path/filename";
 4       
 5//read product and feature
 6Product product = ProductIO.readProduct(productfile);
 7FeatureCollection<SimpleFeatureType, SimpleFeature> featcoll_all = FeatureUtils.loadFeatureCollectionFromShapefile(shapefile);
 8       
 9//get the intersection between product bounds and feature bounds
10Geometry productBounds = FeatureUtils.createGeoBoundaryPolygon(product);
11FeatureCollection<SimpleFeatureType, SimpleFeature> featcoll_part = FeatureUtils.clipCollection(featcoll_all,
12        DefaultGeographicCRS.WGS84, productBounds, DefaultGeographicCRS.WGS84, null, null,new NullProgressMonitor());
13
14// get the bounds of the subset product
15ReferencedEnvelope bounds = featcoll_part.getBounds();
16Rectangle2D rect = new Envelope2D(bounds);
17
18//get image_x, image_y, image_width and image_height of the bounds
19double map_x = rect.getX();//map_x of the upper-left corner (projection coordinate system)
20double map_y = rect.getY() + rect.getHeight();//map_y of the upper-left corner
21DirectPosition map_coord = new DirectPosition2D(map_x, map_y);//create a map coordinate
22GeoCoding geocoding = product.getGeoCoding();
23DirectPosition image_coord =
24        geocoding.getImageToMapTransform().inverse().transform(map_coord, null);//transform map coord to image coord
25int x = (int) Math.floor( image_coord.getCoordinate()[0] );//image_x (image coordinate system)
26int y = (int) Math.floor( image_coord.getCoordinate()[1] );//image_y
27final int resolution=30;//resolution of the landsat8 image
28int w = (int) Math.rint( rect.getWidth()/resolution ) + 1;//image_width
29int h = (int) Math.rint( rect.getHeight()/resolution ) + 1;//image_height
30       
31//create product subset (method 1)
32OperatorSpiRegistry registry = GPF.getDefaultInstance().getOperatorSpiRegistry();
33registry.loadOperatorSpis();
34HashMap<String, Object> parameters = new HashMap<>();
35Rectangle rectangle = new Rectangle(x, y, w, h);
36parameters.put("region", rectangle);
37Product subset = GPF.createProduct("Subset", parameters, product);
38       
39//create product subset (method 2)
40// ProductSubsetDef subsetdef=new ProductSubsetDef();
41// subsetdef.setRegion(x, y, w, h);
42// Product subset=product.createSubset(subsetdef, null, null);

2) as to irregular clip, the following code can be append to the end of the above code:
 1
 2// create vector data node
 3VectorDataNode VecNod = new VectorDataNode("Shape", featcoll_part);
 4subset.getVectorDataGroup().add(VecNod);
 5
 6// create mask
 7Mask maske = subset.getMaskGroup().get(VecNod.getName());
 8maske.readRasterDataFully();
 9
10// get band
11Band band = subset.getBands()[0];
12band.readRasterDataFully();
13
14// clip irregularly
15int width = band.getRasterWidth();
16int height = band.getRasterHeight();
17for (int j = 0; j < height; j++) {
18    for (int i = 0; i < width; i++) {
19        if (maske.getPixelInt(i, j) == 0) {
20            band.setPixelInt(i, j, -9999);//-9999 is the value of the background
21        }
22    }
23}

i'm new in java and beam, if my code is Inefficient, please tell me, thanks in advance!
best wishes!
Zhuyu
Flag Flag