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

[Charles Wivell]

Hi Robert;

I really couldn't tell you what GDALWarp can do, or is for. I'm a control freak.... I write the stuff myself. I only use GDAL to read and write data from/to the really icky file formats.

But warp to me says it uses a polynomial to warp from one projection to another (but maybe I'm wrong), where as the correct (more accurate way) is to use a projection library, like I outlined.

I think if it is doing a warp it might have problems (what is the degree of the polynomial used in the warp...?) and like was pointed out by another replier, since the Mercator projection is not defined at the poles, a polynomial might go off into lala land (technical term) quicker than using the projection library....?

But I'm guessing about GDALWarp....


Chuck Wivell

MN DNR Grand Rapids MN




>>> Robert Rose <robert.w.rose at gmail.com> 6/6/2006 11:00 AM >>>

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


_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org 
http://lists.maptools.org/mailman/listinfo/gdal-dev