1

Given a jp2 file from Sentinel 2-A data and a set of coordinates for a polygon, I need to cut out that polygon from the JP2 file so that I can then convert it to NDVI and run a colour table on it. I have managed to transform it to NDVI and colourise it, but I fail to see which method using GDAL I can using to crop out the polygon.

I would prefer to cut out the polygon before running the other processes as these would be easier on a smaller map, rather than the 100km/sq map from Sentinel in the original JP2 file, but if not possible, then cutting out the polygon from the colourised TIFF file will have to do.

Say I have a list of coordinates, such:

"geo_json": {
    "coordinates": [
      [
        [
          -80.54314497367054,
          26.56162862015057
        ],
        [
          -80.54328137634315,
          26.54771269306852
        ],
        [
          -80.53577924313589,
          26.547773732873694
        ],
        [
          -80.53568830955277,
          26.561669308646813
        ],
        [
          -80.54314497367054,
          26.56162862015057
        ]
      ]
    ],

How would I get that polygon from the jp2 file or the resulting colourised ndvi tif file?

7
  • Do you have something around "coordinates" in your json? Things like Polygon, geometry, Feature, FeatureCollection? If you do, include them as well.
    – user30184
    Jun 29 at 9:47
  • Updated question. Yes, geo_json was wrapping coordinates
    – hyprstack
    Jun 29 at 11:45
  • What is the coordinate system of the jp2 file? You can check it with gdalinfo.
    – user30184
    Jun 29 at 12:51
  • Yup - PROJCRS["WGS 84 / UTM zone 30N",
    – hyprstack
    Jun 29 at 13:11
  • So it is EPSG:32630 epsg.org/crs_32630/WGS-84-UTM-zone-30N.html. It requires one more step.
    – user30184
    Jun 29 at 13:14
2

Unfortunately the json is still not GeoJSON that GDAL could understand directly.

The GDAL utility that supports cropping with a cutline is gdalwarp https://gdal.org/programs/gdalwarp.html. If your input and output rasters are in EPSG:4326 coordinate system like your coordinates, you can crop with this method.

Convert the json coordinates into a WKT polygon and save it as a text file. Name it for example as "cutline.csv".

id,wkt
1,"POLYGON (( -80.54314497367054 26.56162862015057, -80.54328137634315 26.54771269306852, -80.53577924313589 26.547773732873694, -80.53568830955277 26.561669308646813, -80.54314497367054 26.56162862015057 ))"

Test with ogrinfo that GDAL understands the csv file correctly.

ogrinfo cutline.csv -al
INFO: Open of `cutline.csv'
      using driver `CSV' successful.

Layer name: cutline
Geometry: Unknown (any)
Feature Count: 1
Extent: (-80.543281, 26.547713) - (-80.535688, 26.561669)
Layer SRS WKT:
(unknown)
id: String (0.0)
wkt: String (0.0)
OGRFeature(cutline):1
  id (String) = 1
  wkt (String) = POLYGON((-80.54314497367054 26.56162862015057,-80.54328137634315 26.54771269306852,-80.53577924313589 26.547773732873694,-80.53568830955277 26.561669308646813,-80.54314497367054 26.56162862015057))
  POLYGON ((-80.5431449736705 26.5616286201506,-80.5432813763431 26.5477126930685,-80.5357792431359 26.5477737328737,-80.5356883095528 26.5616693086468,-80.5431449736705 26.5616286201506))

Now you can use gdalwarp

gdalwarp -of GTiff -co tiled=yes -co compress=deflate -cutline cutline.csv -crop_to_cutline input.jp2 output.tif

If your cutline is not in the right coordinate system you need either to re-project it, or store into such format that GDAL can automatically recognize the coordinate system. CSV is not such format. The usage of other parameters are documented in gdalwarp manual and in the document page of the output format that you want to use. In my example I used GeoTIFF and the options are documented in https://gdal.org/drivers/raster/gtiff.html.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.