DIY: Reprojection...Re: [Gdal-dev] raster image preprocessing problem.

[Robert Rose]

On Jun 6, 2006, at 6:42 AM, Charles Wivell wrote:

> Hi all;
>
> There's been a lot of mail about reprojecting from one map  
> projection to another. Writing a ditty (scientific term for quick  
> program) to reproject an image from one map projection to another  
> is pretty straight forward, once you have a projection library.
>
> The projection library I use is called GCTP (General Cartographic  
> Transformation package) it is free from the USGS.
>
> ftp://edcftp.cr.usgs.gov/pub/software/gctpc/
>
> I'll pseudo code the process below (the steps are the same no  
> matter what language you use).
>
> 1) Read in the image into a memory buffer (input space)
> 2) Get the current map projection information (Projection, datum,  
> Upper left corner projection parameters, size of pixel in X and Y)
> 3) Frame your output space
>     - Calculate the output space upper left and lower right values
>       - Lat, Lon, projection x/y, output space pixel size
> 4) Create a memory buffer for the output space
> 5) For each pixel in the output space
> 6) Find the projection x/y using the pixel size and offset of the  
> pixel from the upper left corner
> 7) Convert output space projection x/y to lat, lon using GCTP  
> (inverse projection)
> 8) Convert lat, lon to input space projection x/y using forward  
> projection routine in GCTP
> 9) Convert projection x/y into line, sample (row, column)
> 10) Resample DN values in input space near the calculated input  
> space line, sample
> 11) Place resampled DN value into the output space at current pixel  
> location (line, sample)
> 12) Get next output pixel location (back to step 5)
> 13) Write the output space buffer to a file
>
>
> I'm use to using the terms line, sample since my background is from  
> satellite imagery (thats was usually how they collected data a line  
> scanner). It is just the same as row, column or y,x (line is in the  
> y direction and sample is in the x direction)
>
> These are also the basic steps to warp an image from one space to  
> another, just that steps 7 and 8 are replaced with a mathematical  
> relationship, usually a polynomial of some kind. Or if I'm taking  
> raw satellite imagery or photography and projecting it into a map  
> projection then  steps 7 and 8 are replaced with a platform/sensor/ 
> earth mathematical model.
>
> In step 10 you can use Nearest Neighbor (just round to the nearest  
> pixel location) or for a little better looking final image use  
> bilinear resampling.
>
> Oops... I just gave away all my secrets of being a Geospatial  
> programmer.... ;-)


Hi Chuck,  I sent out that email last week about re-projecting a  
Plate Carree image of the Earth into Mercator... (I still haven't  
gotten a sufficient answer btw)..

Sure, I could use the process you describe above to re-project the  
image myself programmatically, but isn't the process you describe  
above exactly what gdalwarp is for?  Or do I have gdalwarp completely  
mistaken?  If so, that would explain why I haven't gotten it to do  
what I want! :-)

Actually I have had success using gdalwarp to convert UTM  
orthoimagery to Mercator, but that's obviously a much smaller area  
than the Earth.  Does gdalwarp have any limitations when you give it  
an image that's suppose to cover the entire Earth?  Maybe that's why  
I can't get this to work.

-robert