How to obtain DN values of Sentinel 1 Data (preferably GRD) given actual location coordinates using Python3 GDAL/rasterio/any_other_library?

Question

location of Svalbard

77.8750° N, 20.9752° E

I have sentinel-1 microwave data (in GRD FORMAT) downloaded from Copernicus. Each image is HH or HV polarized.

How do I calculate the DN value from the GeoTIFF image using Python libraries only?

Will I get multiple pixels in this process?

I could have opened the image in SNAP/QGIS and gone to the respective coordinates using mouse and noted down the DN value.

But, I have around 60 images and I'll get more. I'd like to automate this process.

I have found some relevant libraries:

I'm a CS guy and hence unfamiliar with the terminology.

Output of Gdalinfo:

Driver: GTiff/GeoTIFF
Files: s1a-ew-grd-hh-20170924t192629-20170924t192733-018521-01f35e-001.tiff
Size is 10398, 10860
Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["unnamed",6378137,298.257223560493,
            AUTHORITY["EPSG","4326"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
GCP[  0]: Id=1, Info=
          (0,0) -> (13.2631076300387,-72.5486618031426,2186.18671441078)
GCP[  1]: Id=2, Info=
          (520,0) -> (13.7597223737679,-72.4374481358884,2186.19491506647)
GCP[  2]: Id=3, Info=
          (1040,0) -> (14.2502613555536,-72.3250085612605,2186.20325196814)
GCP[  3]: Id=4, Info=
          (1560,0) -> (14.7347738839402,-72.2113664687458,2186.21172482893)
GCP[  4]: Id=5, Info=
          (2080,0) -> (15.2133063720026,-72.0965461117854,2186.11678827554)
GCP[  5]: Id=6, Info=
          (2600,0) -> (15.6859225359963,-71.9805676769945,2186.11678873375)
GCP[  6]: Id=7, Info=
          (3120,0) -> (16.1526727532145,-71.8634545990536,2186.11678911932)
GCP[  7]: Id=8, Info=
          (3640,0) -> (16.6136141284617,-71.7452287225366,2186.11678944714)
GCP[  8]: Id=9, Info=
          (4160,0) -> (17.0688052189969,-71.6259115158689,2186.1167897312)
GCP[  9]: Id=10, Info=
          (4680,0) -> (17.5183058552513,-71.5055240780671,2186.11678997707)
GCP[ 10]: Id=11, Info=
          (5200,0) -> (17.9621769868422,-71.384087142167,2186.1167901922)
GCP[ 11]: Id=12, Info=
          (5720,0) -> (18.4004805479005,-71.2616210762379,2186.11679038405)
GCP[ 12]: Id=13, Info=
          (6240,0) -> (18.8332793346905,-71.1381458834231,2186.11679055449)
GCP[ 13]: Id=14, Info=
          (6760,0) -> (19.2606368895563,-71.0136812024448,2186.11679070536)
GCP[ 14]: Id=15, Info=
          (7280,0) -> (19.6826173874283,-70.888246309615,2186.1167908432)
GCP[ 15]: Id=16, Info=
          (7800,0) -> (20.0992855238857,-70.7618601227267,2186.11679096799)
GCP[ 16]: Id=17, Info=
          (8320,0) -> (20.5107064067509,-70.6345412063144,2186.11679107975)
GCP[ 17]: Id=18, Info=
          (8840,0) -> (20.9169454561719,-70.50630777672,2186.11679118499)
GCP[ 18]: Id=19, Info=
          (9360,0) -> (21.3180683209978,-70.3771777042252,2186.11679128092)
GCP[ 19]: Id=20, Info=
          (9880,0) -> (21.7141408218795,-70.2471685082238,2186.11679136753)
GCP[ 20]: Id=21, Info=
          (10397,0) -> (22.102986825391,-70.1170548059967,2186.11679145228)
GCP[ 21]: Id=22, Info=
          (0,507) -> (12.8978938221826,-72.4019314084178,2078.27829001099)
GCP[ 22]: Id=23, Info=
          (520,507) -> (13.3923949730691,-72.2916377112545,2078.28607762326)
GCP[ 23]: Id=24, Info=
          (1040,507) -> (13.880988550501,-72.1801102029003,2078.29399506748)
.....
GCP[481]: Id=482, Info=
          (9880,10859) -> (14.4768515859086,-67.2927673124267,-2.49035656452179e-06)
GCP[482]: Id=483, Info=
          (10397,10859) -> (14.8542219755154,-67.1786646643528,-2.41026282310486e-06)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2017:09:24 20:49:42
  TIFFTAG_IMAGEDESCRIPTION=Sentinel-1A EW GRD MR L1
  TIFFTAG_SOFTWARE=Sentinel-1 IPF 002.84
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,10860.0)
Upper Right (10398.0,    0.0)
Lower Right (10398.0,10860.0)
Center      ( 5199.0, 5430.0)
Band 1 Block=10398x1 Type=UInt16, ColorInterp=Gray

I would use rasterio - https://gis.stackexchange.com/q/190423/2856 — user2856, Sep 05 '18 at 07:23
@Luke Where can I find them? Using TiffDump or any other way? — Tarun Maganti, Sep 05 '18 at 12:12
GEOGCS["WGS_84",DATUM["WGS_1984",SPHEROID["unnamed",6378137,298.2572235604902,[AUTHORITY["EPSG","4326"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]] — Tarun Maganti, Sep 05 '18 at 12:42
@Luke I have removed some of the data as it was not accepting entire output — Tarun Maganti, Sep 05 '18 at 13:10
Do you want the pixel value at that location or multiple pixel values in the image? what is the ultimate goal? Are you a Mac user? — GeoMonkey, Sep 07 '18 at 00:04
DN is not calibrated, so comparing values of multiple images is not very accurate. Consider applying calibration to the DN values first: https://sentinel.esa.int/web/sentinel/radiometric-calibration-of-level-1-products — AndyB, Sep 07 '18 at 18:40
@AndyB Yes. Thank you. I wanted to do this step by step. First pick out the DN values and then my next problem would be to calibrate and then pick out the values. — Tarun Maganti, Sep 08 '18 at 06:52

user2856 · Accepted Answer · 2018-09-06T23:36:21.910

Caveat, I have not tested this with a dataset that uses GCPs for georeferencing. If this doesn't work with GCPs, you can use gdalwarp or rio warp to project the dataset.

It also assumes your coordinates are in the same coordinate reference system as your raster dataset (in this case, they should be as you are using longitude and latitude coordinates).

The below code uses the rasterio.DatasetReader.sample method in Python to read values at the specified coordinates. You could also use the rio sample command line tool.

import rasterio

raster = 'test_1band.tif' # 1 band raster
coords = ( #Lon, Lat / X, Y order
    (149.012,-33.7372),
    (149.554,-33.7031),
    (149.606,-33.9477),
    (149.498,-34.2445),
    (149.107,-34.063),
    (149.284,-33.9046),
    (149.373,-33.7753),
    (149.215,-33.8014),
    (148.997,-33.8916),
    (149.189,-34.0239)
)

with rasterio.open(raster) as raster:
    samples = raster.sample(coords)
    for sample in samples:
        print(sample[0])  # sample is a list, with 1 element per raster band
                          # if I'd sampled a 3 band raster, sample would be a 3 element list

For a window you would use a Window/slice and an Affine transform (note completely untested):

with rasterio.open(raster) as raster:
    for coord in coords:
        # raster.transform is an Affine
        col, row = [int(round(cr)) for cr in ~raster.transform * coord]
        # 3x3 window
        win = ((row - 2, row + 2), (col - 2, col + 2))

        # read an ND numpy array
        nd_array = raster.read(window=win)

Works like a charm. I like it that it extracts multiple bands also. Is there a way to open a small kernel of pixels around it? Or should I extend it using a loop? — Tarun Maganti, Sep 06 '18 at 07:37
That it is an instance of an Affine class https://rasterio.readthedocs.io/en/latest/topics/switch.html?highlight=afine#geotransforms — user2856, Sep 08 '18 at 07:12
Thank you. I'll let you know on Monday. No way to access it now — Tarun Maganti, Sep 08 '18 at 07:19

How to obtain DN values of Sentinel 1 Data (preferably GRD) given actual location coordinates using Python3 GDAL/rasterio/any_other_library?

1 Answers1