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

Tue Jun 6 15:38:20 EDT 2006

Ivan;

Yup.... A long time ago and in a galaxy far far away (6 years and at University of North Dakota) I wrote a routine to take Level 1B MODIS 250m GSD (at Nadir) data and project it (if I remember right it was into the LCC projection). Yes, I set up the polynomial (full linear) for each of the blocks of 4 geo-located corner pixels. It was fun and interesting since in the output space driven method as I showed, you don't know which input space block you'll end up in, so I had to test to find which block I was in. Also in addition to framing the whole image, I framed each scan and worked scan by scan. Another fun thing about MODIS, is that at end of scan you have complete overlap between scans (meaning a point on the earth is seen by three scans), so which of the three scans do you use, or do you average..... And at the center there is a scan gap......

Hummmm... Maybe I should have framed each block too.... Might have saved the test.... or you could combine blocks and use higher order polynomials....

The Level 1B MODIS product is still raw geometrically and since they don't give you attitude and position information, you have to use the given geo-located pixels. I would rather have created a platform/sensor/earth model to do the projecting.

I'm calling them geo-located pixels since they are not really tied to the ground (not GCPs). They are created by forward projecting the Line of Sight vector to the ground using a platform/sensor/earth model. They are only as good as the attitude and position information allows.

Higher level MODIS products are already projected, so you can use the steps I listed.

Chuck Wivell

MN DNR Grand Rapids MN

>>> "Ivan Lucena" <ILucena at clarku.edu> 6/6/2006 1:33 PM >>>
Hi Charles,

That was great! 

I didn't miss that class, but forget half of the steps anyway.

Would you be so kind to show as a canonical-text-book-method to
reproject based on a regular or irregular swath of control points?

Talking about a regular grid of control points, some NASA-EOS
HDF-4-Swath image comes with more GCP than you might want. It could be
something like one GPCS over each other 4 pixels! Should we just do the
step 1 through 13 for every each 4x4 window in a 44000x48000 image?

Thanks for your contribution.

Lucena, Ivan

-----Original Message-----
From: gdal-dev-bounces at lists.maptools.org 
[mailto:gdal-dev-bounces at lists.maptools.org] On Behalf Of Charles Wivell
Sent: Tuesday, June 06, 2006 12:42 PM
To: robert.w.rose at gmail.com 
Cc: gdal-dev at lists.maptools.org 
Subject: Re: DIY: Reprojection...Re: [Gdal-dev] raster
imagepreprocessing problem.

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 

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