Extracting a polygon from slanted Sentinel satellite image using R

Question

I have a Sentinel image, the coordinates of the boundary are given. I want to extract/subset a polygon from it. I am using R and don't have any prior experience in GIS.

I have downloaded some Sentinel 1 images from Copernicus https://scihub.copernicus.eu/dhus/#/home. The associated KML file has the coordinates of the image corners.

<coordinates>-118.362823,34.353683 -121.078079,34.752659 -120.786476,36.256836 -118.020294,35.860096</coordinates>

If I plot the coordinates on a map using DaftLogic Area Calculator it looks like the screenshot below:

The shape drawn above isn't parallel/perpendicular to the latitude and longitude lines, its inclined because of the satellite orbit. Also, two of the sides have a difference of ~1 km. So the shape is not in rectangular format.

I am opening the satellite image using raster package in R. The extent(georeferencing) of the image/raster doesn't have the actual geo coordinates. But, it represents the dimensions of this image.

library(raster)

s <- brick('s1a-iw-grd-vv-20161030t230720-20161030t230745-013726-016066-001.tiff')

s
#class       : RasterBrick 
#dimensions  : 16716, 25485, 426007260, 1  (nrow, ncol, ncell, nlayers)
#resolution  : 1, 1  (x, y)
#extent      : 0, 25485, 0, 16716  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#data source : s1a-iw-grd-vv-20161030t230720-20161030t230745-013726-016066-001.tiff 
#names       : s1a.iw.grd.vv.20161030t230720.20161030t230745.013726.016066.001 
#min values  : 0 
#max values  : 65535

extent(s)
# class       : Extent 
# xmin        : 0 
# xmax        : 25485 
# ymin        : 0 
# ymax        : 16716

There is a function to assign the extent, setExtent(). This would have worked if the shape/polygon was rectangular and not inclined.

What I need to do after the correct georeferencing is done, is to get a subset of the raster based on the polygon using extract(). Since I am not able to set the extent properly, I am not able to do so.

Why R? I can use any other language as well. I am using a language because I have to repeat the process with multiple images.

One more problem to deal with is, some of the images are mirrored image, this adds to the problem. Note for mirrored image https://sentinel.esa.int/web/sentinel/faq

Is this the right approach? Should I work with other package or language? If somebody has already done it, please help.

Answer 1

The raster package in R does not deal well with non-axis aligned GeoTIFFs. The way they work is that the file contains a number of Ground Control Points (GCPs). These map pixel coordinates to world coordinates.

If you have gdal installed fully you can use the gdalinfo command to see the GCPs. Here's the output from a sample raster of mine:

GCP[  0]: Id=1, Info=
          (1,1) -> (0,0,0)
GCP[  1]: Id=2, Info=
          (457,1) -> (3,3,0)
GCP[  2]: Id=3, Info=
          (457,322) -> (0,6,0)
GCP[  3]: Id=4, Info=
          (1,322) -> (-3,3,0)

I did download a file from Copernicus and it had a lot of GCPs.

If I load my test raster and plot it I don't see it in the right place:

> r = raster("./screenA5.tiff")
> plot(r)

(Note the image is a raster of one of our departmental leaflets!)

R has ignored the GCPs and plotted it in pixel coordinates.

What I can do is warp the raster using the gdalwarp command and that will use the GCPs to create a new, axis-aligned raster:

$ gdalwarp  screenA5.tiff screenA5warp.tiff
Processing input file screenA5.tiff.
0...10...20...30...40...50...60...70...80...90...100 - done.

Load this in and plot and it is rotated correctly and projected to lat-long:

> r2 = raster("./screenA5warp.tiff")
> plot(r2)

The corners have been padded out with NA missing values. You may have to adjust the parameters for gdal_warp depending on your data, and it also might take ages to warp the sort of 1Gb raster I got from the Copernicus site...