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:

Monday, June 16, 2014

Converting Coordinates between map projections using Python

Python can be used very nicely to quickly convert coordinates from one projection into another using the pyproj package, as usual available here. This example does it from command line, but you can built it into your script, of course, as well. The easiest way is to find the EPSG code for your projections, for example on this page. You then initialize the projections you want, in this example unprojected lat/long into EPSG:3575:

>>> EPSG4326 = pyproj.Proj("+init=EPSG:4326")
>>> EPSG3575 = pyproj.Proj("+init=EPSG:3575")


Alternatively you can use the Proj.4 definition, so the following two lines are doing the same, you can use either of them
EPSG3995 = pyproj.Proj("+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ")
EPSG3995 = pyproj.Proj("+init=EPSG:3995")


Now you can define the input coordinates and transform them, in this example you transform values from EPSG:4326 to EPSG:3575:

>>> lat = 79
>>> long = -5
>>> x, y = pyproj.transform(EPSG4326, EPSG3575, long, lat)
>>> x,y
(-317466.6103301885, -1184801.519458934)

Official documentation can be found here