3

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?

Improve this question
5
  • gdal has built-in support for this, though my experience is limited. See here for an example: lists.osgeo.org/pipermail/gdal-dev/2012-February/031756.html I also think gdalwarp infile outfile.tif will do something.
    – mdsumner
    Commented Jun 25, 2014 at 5: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: lists.osgeo.org/pipermail/gdal-dev/2012-June/033147.html
    – Max
    Commented Jun 25, 2014 at 7:03
  • Can you point to a specific file? I am keen to check it out
    – mdsumner
    Commented Jun 25, 2014 at 8:43
  • 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
    Commented Jun 25, 2014 at 11:46
  • My email is in my profile. Send a link to the original source rather than the file if possible
    – mdsumner
    Commented Jun 25, 2014 at 13:20

1 Answer 1

Reset to default
2

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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.