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

[David Hoese]

On 6/13/12 10:00 AM, gdal-dev-request at lists.osgeo.org wrote
> David Hoese wrote
>> >  
>> >  I'm attempting to project satellite data to an AWIPS grid.  Assuming I
>> >  have the grid defined correctly with proj4 parameters, I've been using
>> >  "+proj=latlong +datum=wgs84" as my s_srs parameters.  Does that sound
>> >  correct for satellite data?  So my command line call (still for testing)
>> >  becomes:
>> >  
> Is it publicly available satellite data? I will give it a go myself if i
> have some sample data.
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.
>
>
> David Hoese wrote
>> >  
>> >  gdalwarp -s_srs "+proj=latlong +datum=wgs84" \
>> >  -te 59.844 -123.044 14.335 -65.091 \
>> >  -ts 5120 5120 \
>> >  -t_srs "+proj=lcc +lon_1=-95.000  +lon_0=-113.133 +lat_2=25.0001
>> >  +lat_1=24.9999 +lat_0=16.369" \
>> >  -geoloc -of GTIFF \
>> >  image.vrt out.tif
>> >  
>> >  The output grid is 5120x5120 and I wasn't sure if -te was lat/lon
>> >  points, but that's what I used.  -te has lat/lon of upper left pixel of
>> >  the grid and lat/lon of the lower right pixel of the grid.  Does any of
>> >  this look correct?  Thanks for the help.
>> >  
> Geographic WGS84 (lat/lon) is an obvious choice if you're not sure, but its
> not necessarily correct. As you can see, your Lambert target projection has
> coordinates in a similar value range, so you cannot be sure just based on
> the values alone.
> The syntax for gdalwarp seems alright to me. Because you're outputting to a
> rather large area (geographically), it might be useful to test it for a
> small subset first to rule out difficulties in finding a proper
> transformation, although im not sure if it could help in this case.
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