Map project a raster having separate latitude and longitude raster bands

Question

I have a satellite image where one of the channels is a raster having the longitude values, another channel a raster containing the latitude values.

What would be the best way for map projecting such files? Is there some standard routines or tools (I use Python / gdal preferrably)?

My only thought so far is to loop through the pixels, convert the lat/long value of a given pixel and drop it into the closest matching pixel of a raster already in the projection I want to have (where it seems not that straightforward to actually find the closest matching value) -- but somehow I wonder if there is a better and quicker way?

gdal has built-in support for this, though my experience is limited. See here for an example: http://lists.osgeo.org/pipermail/gdal-dev/2012-February/031756.html I also think gdalwarp infile outfile.tif will do something. — mdsumner, Jun 25 '14 at 05:34
thanks a lot -- had not seen that gdal supports this and that the word to search for is "geolocation arrays" --- but there is very little info on how to use it. One has to create descriptive files like this, but have not gotten it to work yet: http://lists.osgeo.org/pipermail/gdal-dev/2012-June/033147.html — Max, Jun 25 '14 at 07:03
I use AMSR2 brightness files at gcom-w1.jaxa.jp and created a modified vrt file according to the link in my previous comment (I think these have some errors) ... not sure how to send private messages, but I could send you the file. gdalwarp does something, but not correctly yet ... — Max, Jun 25 '14 at 11:46
My email is in my profile. Send a link to the original source rather than the file if possible — mdsumner, Jun 25 '14 at 13:20

Max · Accepted Answer · 2014-06-26T08:53:41.660

The comments above lead me towards an answer (but if there are other possibilities, please let me know):

gdal can utilise bands containing longitude and latitude values, which is called there geolocation arrays. The answer below may not be thorough and I myself don't understand all the details:

The original image is a orbital swath in a raster (AMSR2 data) and there is a raster band containing latitude and a raster band containing longitude values:

enter image description here

This link tells me that you have to create vrt-files (as described here) describing the image raster, the lat band and the lon band, also connnecting the lat/lon band to the image

So I create one file called GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt:

<vrtdataset rasterxsize="486" rasterysize="2017">
   <metadata domain="GEOLOCATION">
     <mdi key="X_DATASET">C:\Users\max\Documents\lon.vrt</mdi>
     <mdi key="X_BAND">1</mdi>
     <mdi key="Y_DATASET">C:\Users\max\Documents\lat.vrt</mdi>
     <mdi key="Y_BAND">1</mdi>
     <mdi key="PIXEL_OFFSET">0</mdi>
     <mdi key="LINE_OFFSET">0</mdi>
     <mdi key="PIXEL_STEP">1</mdi>
     <mdi key="LINE_STEP">1</mdi>
   </metadata> 
   <vrtrasterband band="1" datatype="UInt16">
   <metadata>
    <simplesource>
      <sourcefilename relativetovrt="1">HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Brightness_Temperature_(89.0GHz-A,H)</sourcefilename>
      <sourceband>1</sourceband>
      <sourceproperties blockxsize="486" blockysize="1" datatype="UInt16" rasterxsize="417" rasterysize="2017">
      <srcrect xoff="0" xsize="486" yoff="0" ysize="2017">
      <dstrect xoff="0" xsize="486" yoff="0" ysize="2017">
    </dstrect></srcrect></sourceproperties></simplesource>
  </metadata></vrtrasterband>
</vrtdataset>

another file called lon.vrt

<vrtdataset rasterxsize="486" rasterysize="2017">
  <vrtrasterband band="1" datatype="Float32">
    <metadata>
    <simplesource>
      <sourcefilename relativetovrt="1">HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Longitude_of_Observation_Point_for_89A</sourcefilename>
      <sourceband>1</sourceband>
      <sourceproperties blockxsize="486" blockysize="1" datatype="Float32" rasterxsize="486" rasterysize="2017">
      <srcrect xoff="0" xsize="486" yoff="0" ysize="2017">
      <dstrect xoff="0" xsize="486" yoff="0" ysize="2017">
    </dstrect></srcrect></sourceproperties></simplesource>
  </metadata></vrtrasterband>
</vrtdataset>

and a file called lat.vrt

<vrtdataset rasterxsize="486" rasterysize="2017">
  <vrtrasterband band="1" datatype="Float32">
    <metadata>
    <simplesource>
      <sourcefilename relativetovrt="1">HDF5:GW1AM2_201301311114_050A_L1SGBTBR_1110110.h5://Latitude_of_Observation_Point_for_89A</sourcefilename>
      <sourceband>1</sourceband>
      <sourceproperties blockxsize="486" blockysize="1" datatype="Float32" rasterxsize="486" rasterysize="2017">
      <srcrect xoff="0" xsize="486" yoff="0" ysize="2017">
      <dstrect xoff="0" xsize="486" yoff="0" ysize="2017">
    </dstrect></srcrect></sourceproperties></simplesource>
  </metadata></vrtrasterband>
</vrtdataset>

Now I can use gdalwarp, note that you call the vrt file rather than the hdf!

gdalwarp -geoloc -t_srs EPSG:4326 C:\GW1AM2_201301311114_050A_L1SGBTBR_1110110.vrt C:\test100.tif

The resulting image is this one:

enter image description here

See also this page for this description and possible updates.

Map project a raster having separate latitude and longitude raster bands

1 Answers1

Linked