[Gdal-dev] Re: World Wind

[Frank Warmerdam]

On 12/10/05, Chris Maxwell <cmaxwell at arc.nasa.gov> wrote:
> My name is Chris Maxwell, and I'm the lead software developer on the NASA
> World Wind project.  I was told you were the guy to speak to when it comes
> to Proj and GDAL.  Let me first say that our group readily uses the GDAL
> suite of utility programs for much of our data process flow, and are eager
> to try to use it inside of the actual World Wind client software (in our
> next version).  I've very impressed with the performance of the GDAL
> programs, and they easily compete with quite a few high-end commercial
> packages we've used.

Chris,

Thanks.

> I've been trying to figure out how to reproject the NASA Blue Marble from
> it's native Plate Carree projection into Polar Stereographic projection (for
> just the upper latitudes for the north and south pole).  I'm running into
> problems with this, and I'm wondering if you have some experience that you
> could pass my way.  Is this something that GDAL cannot handle, or is it more
> of a user-error?

There are some tricky aspects to warping around the polar point.

I have an old blue marble with georeferencing that looks like this:

warmerda at gdal2200[100]% gdalinfo world.tif
Driver: GTiff/GeoTIFF
Size is 2048, 1024
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000,90.000000)
Pixel Size = (0.17578125,-0.17578125)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0'0.00"W, 90d 0'0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0'0.00"W, 90d 0'0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0'0.00"E, 90d 0'0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0'0.00"E, 90d 0'0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Band 1 Block=2048x4 Type=Byte, ColorInterp=Red
Band 2 Block=2048x4 Type=Byte, ColorInterp=Green
Band 3 Block=2048x4 Type=Byte, ColorInterp=Blue

I tried reprojecting a north pole region out of this, and found I needed
the following to be successful:

gdalwarp -t_srs EPSG:32661 -te 0 0 4000000 4000000 world.tif
northpole.tif -ts 500 500 -wo SAMPLE_GRID=YES -wo SOURCE_EXTRA=50

The "-t_srs EPSG:32661" sets the projection to Universal Polar Stereographic
North (which has a false easting and northing of 2000000m).

The -te establishes a default extents, and the -ts says it should be
generated with an output size of 500 pixels by 500 lines.    I found
that gdalwarp gets confused about pixel sizes in the polar regions so
I set the extents and size explicitly.  You could also use the -tr
instead of -ts to set a target pixel resolution instead of target file
size.   The 0 0 4000000 4000000 box gives a 4000km by 4000km
region centered on the pole.

The -wo SAMPLE_GRID=YES forces gdalwarp to determine the source
region based on a sampling of locations through the output square, not
just around the edges.  If you just walk the edges of a box centered on
the pole, you won't realize you need the data that is very near the pole.

The -wo SOURCE_EXTRA adds 50 pixels to any source rectangle
computed to make up for the inperfectness of the sampling scheme.

With that I was able to get a nice result out.

  http://home.gdal.org/~warmerda/northpole.png

I have taken the liberty of cc:ing this response to gdal-dev, since
other folks also run into problems with these sorts of warps.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent