0

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 - gis.stackexchange.com/q/190423/2856 – user2856 Sep 5 '18 at 7:23
  • @Luke Where can I find them? Using TiffDump or any other way? – Tarun Maganti Sep 5 '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 5 '18 at 12:42
  • @Luke I have removed some of the data as it was not accepting entire output – Tarun Maganti Sep 5 '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? – Nathan Thomas Sep 7 '18 at 0:04
3

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 6 '18 at 7:37
  • 1
    @TarunMaganti see edit. – user2856 Sep 7 '18 at 7:17
  • What do you mean by raster transform is an Affine? – Tarun Maganti Sep 8 '18 at 6:53
  • That it is an instance of an Affine class rasterio.readthedocs.io/en/latest/topics/… – user2856 Sep 8 '18 at 7:12
  • Thank you. I'll let you know on Monday. No way to access it now – Tarun Maganti Sep 8 '18 at 7:19

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.