[mapserver-users] Examples using projection?

[Ed McNierney]

Mike -

You're making the (incorrect) assumption that the minima and maxima of
the output projection are related to the minima and maxima of the input
projection, and that's not necessarily true.  Blue Marble Geographics
has a unfixed bug in their GeoTransform DLL because they make the same
assumption.  I mention their bug because it's a relatively simple
illustration.

USGS topographic maps (DRGs) are raster maps in the UTM projection.
Since each map is bounded by lines of latitude and longitude, the map
image is not a rectangle in UTM space.  In most cases (all examples here
are in the Northern hemisphere) it looks a lot like a rotated rectangle
- west of the UTM zone centerline, the maps look like they're rotated
clockwise, while east of the centerline, they look like they're rotated
counterclockwise.  The minimum distortion occurs right at the centerline
(the central meridian of the UTM zone) - lines of latitude are indeed
perfectly horizontal at the central meridian, and start curving up on
either side of it.

As a result, the minimum point in UTM space for all USGS DRGs is EITHER
the southeast or the southwest corner, depending on whether the DRG is
in the western half or the eastern half of the UTM zone, respectively.
Since the 1:24,000 and 1:25,000 DRGs are 7.5 minutes each, the central
meridian of the UTM zone (an integer degree boundary) is at the right
edge or the left edge of the two DRGs that abut it.  Since no map ever
straddles the centerline, a transformation such as you describe in your
email works perfectly.  The 1:100,000 DRGs are 1 degree wide and
everything's OK.

Except.... for the 1:250,000 DRGs.  These maps are 1 degree high by 2
degrees wide, and some of them do indeed straddle the centerline of a
UTM zone.  For example, at 38 degrees North there's a DRG that runs from
98 degrees West to 100 degrees West, crossing the 99-degree-West central
meridian of UTM Zone 14.  If you look at the original DRG, the map image
curves up to the right and left of the centerline.  That is, the
northern neatline (39 degrees North) is a shallow curve that looks like
a smile; the southern neatline is almost identical.  And it's that
southern neatline (in this example) that jumps up and bites Blue Marble
(and you, too).

Since the southern neatline looks like a smile, its minimum value is in
the MIDDLE OF THE LINE!  That is, it's not the southeast corner and it's
not the southwest corner.  If you try to use Blue Marble's tool to clip
the collar from this map, the bottom will be truncated.  That's because
they presume the minimum output coordinate is the minimum of the
projected southeast and the projected southwest corners.  But it's not
either of those - both of those values are equal and they're too high.

To solve this in the general case, you have to reproject EVERYTHING or
do some pretty darn clever math.  If anyone would like me to email them
the 4.5 MB TIFF/LZW file used in this example, let me know (a picture is
indeed worth these thousand words).

	- Ed

Ed McNierney
Chief Mapmaker
TopoZone.com
ed at topozone.com
(978) 251-4242


-----Original Message-----
From: Mike Neuman [mailto:mcn at engarde.com]
Sent: Monday, July 23, 2001 3:45 PM
To: Stephen Lime; mapserver-users at lists.gis.umn.edu
Subject: RE: [mapserver-users] Examples using projection?



> How are you projecting from old extent to 
> new? Doing just the corners will not work. 

  What is the proper way to do this? My formula is something like:

  lccX1,lccY1 = project(minlat, minlon)
  lccX2,lccY2 = project(minlat, maxlon)
  lccX3,lccY3 = project(maxlat, minlon)
  lccX4,lccY4 = project(maxlat, maxlon)

  new_extent(minx) = MIN(lccX1, lccX2, lccX3, lccX4)
  new_extent(maxx) = MAX(lccX1, lccX2, lccX3, lccX4)

  new_extent(miny) = MIN(lccY1, lccY2, lccY3, lccY4)
  new_extent(maxy) = MAX(lccY1, lccY2, lccY3, lccY4)

  Thanks!

-Mike