Thursday, June 26, 2014

Reading and Map projecting raster data with geolocation arrays using gdal

The following describes how to geocode a raster band, whose the geolocation information is stored in one channel containing latitude values for each pixel, and another channel containing longitude values for each pixel. Opening the image, you see the unprocessed orbital swath.


Information on the hdf image can be retrieved with gdalinfo:

gdalinfo C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5

The metadata printed out contains the names of the subchannels, as seen here:
and information about a specific channel can be retrieved by using gdalinfo on exactly this subchannel-name from the metadata:

gdalinfo HDF5:"C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5"://Brightness_Temperature_(89.0GHz-A,V)

The tricky part is now, to transfer such raster bands into a map projection.

gdalwarp can  read separate latitude and longitude rasters, which is called geolocation arrays in gdal. the comment on this page brought me a step further although I think that the vrt-code and the commands given are not completely correct. One has to create a vrt file describing the image raster, also referring to the latitude and longitude bands. Then in addition a vrt file describing the latitude band and a vrt file describing the longitude band.

So I create a vrt-file (see here for info) named GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt:


   
     C:\Users\max\Documents\lon.vrt
     1
     C:\Users\max\Documents\lat.vrt
     1
     0
     0
     1
     1
    
   
   
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Brightness_Temperature_(89.0GHz-A,H)
      1
      
      
      
    
  




Then a file named lon.vrt :


  
    
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Longitude_of_Observation_Point_for_89A
      1
      
      
      
    
  



Finally a file named lat.vrt :

  
    
    
      HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Latitude_of_Observation_Point_for_89A
      1
      
      
      
    
  



With these files created and placed in the same directory, I can use gdalwarp, but note that you are calling the vrt file rather than the hdf file!
gdalwarp -geoloc -t_srs EPSG:4326 C:\Users\max\Documents\GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt C:\Users\max\Documents\test100.tif


The result is the following swath:


In order to get it onto the polar stereographic projection, in my case the NSIDC EPSG:3411 projection, I first subset the image:


gdal_translate -projwin -94 90 40 35 C:\Users\max\Documents\test100.tif C:\Users\max\Documents\test101.tif


which results in this image:



now I reproject into EPSG:3411:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3411 C:\Users\max\Documents\test101.tif C:\Users\max\Documents\test102.tif

resulting in this image:

11 comments:

  1. Are you doing this all in qGIS?

    ReplyDelete
  2. This above is done in pure Python and gdal directly, I use QGIS mostly for display and maps, less for processing

    ReplyDelete
  3. Hi Max - Thanks for your blog. Excellent to see a hands-on example, for me as a GDAL newbie. What would you do if gdalwarp returns with: ERROR 1: Unable to compute a GEOLOC_ARRAY based transformation between pixel/line and georeferenced coordinates
    Thanks

    ReplyDelete
  4. Hi Martin, thanks for that :-) what i usually domis just googling the error message. That mostly gives me answers on stackexchsnge such as http://gis.stackexchange.com/questions/137331/georeferencing-uneven-irregularly-gridded-rasters
    Maybe that or another answer there helps?

    ReplyDelete
  5. Hi Max, I have tried it in gdal, but it said that "HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Brightness_Temperature_(89.0GHz-A,H)" does not exist in the file system, and is not recognised as a supported dataset name. Why?

    ReplyDelete
    Replies
    1. my immediate thought that you possibly specify the wrong file path or are not using double slash // in the path?

      Delete
  6. Hello Max,

    I'm trying to implement this procedure for another data set, SMAP L1C_S0_HiRes, with different geolocation arrays and projection but am hitting a wall!

    Would simply updating all the above .vrt files be enough?

    I am getting a few "ERROR 1: Too many points (441 out of 441) failed to transform" errors and data values from certain data fields turn out to be zero.




    ReplyDelete
    Replies
    1. it's been a while I have been doing this, but your vrt file should contain the correct information describing your dataset ( http://www.gdal.org/gdal_vrttut.html ) Since gdalwarp is trying to do something, some info given to it seems to be incorrect(?)

      Sometimes I got your error message when trying to reproject into an "impossible" projection for my dataset, like by mistake trying to reproject an Arctic location into an Antarctic projection.

      Delete
    2. Thanks for the note. It looked like the source and target projections weren't defined correctly.

      A follow-up question! Which do you think is a better option: (a) make multiple .vrt files for different data files, perform gdalwarp the merge them together, or (b) make a single .vrt file linking all the data sets, the gdalwarp? Each file is a single swath, covering part of the globe.

      Thanks!

      Delete
  7. Hi Max,
    Thanks for the tutorial. I would like to automate the the process using python. Is there a way to perform it?
    I trying to search for ways to call osgeo shell commands from python script or alternative gdal api in python.

    Will be grateful if u suggest solution.

    ReplyDelete
    Replies
    1. I guess that would be my post just following right after that one above: http://geoinformaticstutorial.blogspot.no/2014/07/reading-amsr-2-data-into-geotiff-raster.html

      I see the referred github link has vanished, I can try to fix that ... but the post may indicate what to do (for AMSR at least)

      Delete