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

[David Hoese]

Hi Rutger (long email below),

I figured out my projection problems.  To increase google search 
accuracy I'm going to try to include as many relevant terms as I 
searched for to describe what happened.

I was attempting to project VIIRS binary image data that I had extracted 
from the SDR hdf files.  This extraction left me with three files, the 
binary image file, the binary latitude file, the binary longitude file.  
I needed to project the data into AWIPS grids, east 211 for testing.  
The grid is 5120 x 5120 with about 1km equal area pixels.  I had 
previously been using ms2gt (a mapx based toolbox, gpd files to specify 
projection parameters) to do this, since it was the first thing 
suggested, and had a png image of what the remapped data should look 
like.  To be able to use GDAL I needed to convert the ms2gt/mapx 
projection parameters into a proj4 string.  Seemed pretty 1:1 and 
generally it was.

After spending way too much time trying to figure out why the images 
weren't looking the same, I decided to look at the ms2gt/mapx source and 
figure out how it was using the parameters to see how that differs from 
the expected behavior.  The version of ms2gt I was using was very old 
(2001), so was the mapx library.  The documentation mentions that you 
can specify the equatorial and polar radius so that's what I had been 
using in proj4.  I checked the source code of mapx...and it IGNORES it, 
it doesn't even give you an error about an unknown parameter.

So after figuring this out, I changed the proj4 parameters to say that 
+a and +b were the same.  I then made the origin of the projection the 
origin specified on the NOAA page (from previous email) and used pyproj 
to calculate the grid origin offset in meters to use for the -te (target 
extent) in the gdalwarp call.  And finally it worked, the 2 images 
matched, at least visually.  Now I just need to work on bow tie 
correction (even though "Bow ties are cool"), but some coworkers may 
know how to handle this...or gdal-dev will get another email from me.  
I'll also need to meet my original goal and do the projecting with the 
Python bindings.

Thanks a lot Rutger.

-Dave

P.S. Final gdalwarp call, may not be exactly correct for AWIPS since I 
haven't tested it yet:
gdalwarp -s_srs "EPSG:4326" -t_srs "+proj=lcc +ellps=clrk66 +a=6378137 
+b=6378137 +e=0.0818191910435 +lat_0=24.9999 +lon_0=-95 +lat_1=24.9999 
+lat_ts=25.0001 +units=m +no_defs" -te -1952976.3246 -828316.5944 
3248431.6754 4373091.4056 -of GTiff -geoloc -overwrite -r bilinear -ts 
5120 5120 image.vrt out.tif

Note: Clark1866 is what ms2gt says its default is. And the 24.9999 and 
25.0001 are artifacts from previous iterations.

On 6/14/12 7:48 AM, gdal-dev-request at lists.osgeo.org wrote:
> 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