[gdal-dev] reproject python numpy binary swath/lat/lon

[Rutger]

Hey David,


David Hoese wrote
> 
> On 6/13/12 10:00 AM, gdal-dev-request at .osgeo wrote
> I'm not sure at the moment, it's VIIRS data.  I'm playing around with 
> different existing remapping tools, I have one that works called 
> ms2gt(ll2cr/fornav) which you might be able to google for the old 
> versions, but it's slightly outdated and is single purpose.  I would 
> like to use GDAL if it can handle it, because of the various formats 
> that it supports as output.
> 

I actually already did some test with GDAL geolocation arrays and VIIRS
data. It work flawless for my target area; utm31 for the Netherlands. My
VRT's looked like this:

lat.vrt
<VRTDataset rasterXSize="6400" rasterYSize="1536">
  <VRTRasterBand dataType="Float32" band="1">
    <Metadata />
    <SimpleSource>
      <SourceFilename
relativeToVRT="0">HDF5:GITCO_npp_d20120315_t1208597_e1210238_b01974_c20120321152146856862_noaa_ops.h5://All_Data/VIIRS-IMG-GEO-TC_All/Latitude</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="6400" RasterYSize="1536"
DataType="Float32" BlockXSize="6400" BlockYSize="1536" />
      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

lon.vrt
<VRTDataset rasterXSize="6400" rasterYSize="1536">
  <VRTRasterBand dataType="Float32" band="1">
    <Metadata />
    <SimpleSource>
      <SourceFilename
relativeToVRT="0">HDF5:GITCO_npp_d20120315_t1208597_e1210238_b01974_c20120321152146856862_noaa_ops.h5://All_Data/VIIRS-IMG-GEO-TC_All/Longitude</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="6400" RasterYSize="1536"
DataType="Float32" BlockXSize="6400" BlockYSize="1536" />
      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

band5.vrt
<VRTDataset rasterXSize="6400" rasterYSize="1536">
   <Metadata domain="GEOLOCATION">
     <MDI key="X_DATASET">lon.vrt</MDI>
     <MDI key="X_BAND">1</MDI>
     <MDI key="Y_DATASET">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 dataType="UInt16" band="1">
    <Metadata />
    <SimpleSource>
      <SourceFilename
relativeToVRT="0">HDF5:SVI05_npp_d20120315_t1208597_e1210238_b01974_c20120316203244033373_noaa_ops.h5://All_Data/VIIRS-I5-SDR_All/BrightnessTemperature</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="6400" RasterYSize="1536"
DataType="UInt16" BlockXSize="6400" BlockYSize="1536" />
      <SrcRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
      <DstRect xOff="0" yOff="0" xSize="6400" ySize="1536" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

And then reprojection with gdalwarp like:
gdalwarp -s_srs "epsg:4326" -t_srs "+proj=utm +zone=31 +datum=WGS84
+units=m" \
-tr 250 250 -te 392875 5574125 808875 5940125 band5.vrt band5.tif

This works without any errors for me, and the geolocation seems accurate.
But as said, i didnt like the resulting image quality. I have attached a
simple comparison between a gdalwarp nearest neighbor, bilinear and
Pyresample gaussian reprojection. The gdalwarp reprojection seem very
patchy, and the bilinear flag (same with cubic) doesnt seem to have much
affect. Im not sure whether this is to be expected with the current
implementation or if i made any mistakes while specifying the VRT files.


David Hoese wrote
> 
> So I've played around with this for a while now and have gotten 
> somewhere.  I was able to get data into the geotiff, but with a lot of 
> guess work and small tweaking.  My t_srs is now:
> "+a=6378137 +lon_1=-95.000 +b=6356752.314 +lon_0=-123.044 +proj=lcc 
> +lat_1=25.000 +lat_0=59.844 +ellps=WGS84 +no_defs +wktext +units=m"
> 
> All of which are probably not necessary, but at this point it doesn't 
> hurt.  The origin in that projection is the upper left corner of my 
> grid.  So I now specify -te as something like:
> "-te 0.0 -5120000.0 5120000.0 0.0"
> 
> Which gets me something in the geotiff.  I've played around with 
> different extents and my best result comes from me manually changing the 
> extents to match a picture I have from another system.  The grid I'm 
> trying to match to is the first one on this page:
> http://www.nws.noaa.gov/noaaport/html/icdtb48_2.html
> 
> I have tried using pyproj and plugging in different corner points to get 
> the extent points in meters, but when I use those I get skewed images 
> (the grid should be 5120 x 5120).
> 
> So at this point, I just need to find a way to get the extents figured 
> out in a programmatic way (like 1km * 5120 grid pixels = 5120000 
> meters).  If you have any more advice that'd be great, otherwise I'll 
> just have to keep playing around and asking people around my work.  
> Thanks a lot for the help so far.
> 
> -Dave
> 

Since it worked for my case, i think its likely that the problems you are
having are related to the projection and area definition. I have no
experience with the AWIPS grid you are mentioning. What i would do in such a
case is to get a sample file with the correct projection and boundingbox,
and just read the metadata with:
gdalinfo -proj4 samplefile.tif

Optionally add some of the -norat -noct etc flags to hide the large tables.
The proj4 string can be directly used by other gdal tools. 


Regards,
Rutger

--
View this message in context: http://osgeo-org.1560.n6.nabble.com/gdal-dev-reproject-python-numpy-binary-swath-lat-lon-tp4978609p4981318.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.