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

[Charles Wivell]

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.... ;-)


Chuck Wivell

MN DNR Grand Rapids MN








>>> "Zhonghai Wang" <zhonghaiw at gmail.com> 6/6/2006 5:34 AM >>>
Hi there,

I have a couple of GIFs with wld files, but the original Projection is not
the one I want, I tried to use the Proj4 to reproject them on-the-fly, but
it takes really long to perform the map request, therefore I tried
to preprocess the images with GDAL tools, but I've got some problems here:

1. problem with gdalwarp: trying to mosaic the images into one output file,
but only the first image has the correct image color, the other images are
rendered with strange colors, is this a common problem if using gdalwarp to
mosaic images?

2.problem with gdal_merge.py: after trying the gdalwarp command, I've tried
the gdal_merge.py command to merge all the images together with the -pct
parameter, unfortunately, the merged image is the same as the first try. is
it possible to merge images using the command with the correct color table?

BTW: I can not use the gdal_translate command, since there are gaps between
the reprojected images.

thanks for your info.

zhonghai