Wednesday, February 10, 2010

GDAL 1.7.0 is dead, long live GDAL 1.7.1

Due to a serious bug in GDAL 1.7.0 resulting in all Erdas Imagine files produced by the release being unreadable to any non-GDAL applications, including Erdas Imagine, the project has decided to retract the 1.7.0 release. In it's place today we have issued a GDAL/OGR 1.7.1 release. I sincerely hope that use of 1.7.0 in the wild will not end up producing unreadable Erdas Imagine files that haunt us for years.

The introduction of the bug was due to a very subtle issue, and it didn't affect GDAL's ability to read the files. So our extensive regression tests for this file format did not detect any problems. But the lesson would appear to be that we still need mechanisms to do cross testing with other packages as often as possible. We are looking into ways of more vigorous and organized cross testing at major releases.

Tuesday, February 9, 2010

World Mapping

Over much of the last decade I have fought with world mapping issues in MapServer. This falls into three broad areas - problems with requests crossing the dateline, problems with requests including the poles, and problems with requests extending beyond the domain of validity of a projection. Typically these are not a big issue when working on local mapping projects, but do become issues when doing world level views, or when exploring areas far from the domain of validity of some of the data.

Several years ago I made a modest effort to collect some hints and tips, and knowledge about the issues into a topic in the old MapServer wiki. With the planned decommissioning of the old mapserver web site at UMN (and in light of recent work in this area), I have migrated the old wiki topic to the MapServer Trac Wiki as OldWorldMappingIssues, and started an updated wiki topic on the issues.

The problems I have been investigating in recent months on behalf of a client, are related to providing WMS service in web mercator projection with MapServer, but cascading to another WMS providing data in EPSG:4326 (geographic). The initial problem report was with WMS requests in the Alaska/Siberia area that were oddly pixelated. The problem was that the request region from -21000000 to -17000000 easting was getting reprojected to geographic as -179 to 179 in longitude. This was occurring because value further east than the dateline (roughly -20000000 meters) were getting returned as positive longitudes since PROJ.4 returns longitude values in the range -180 to 180. This transformed BBOX was sent to the cascaded WMS, but with a width and height that reflected the approximate size of the original request. The result was that the map pulled from the cascaded WMS had very poor horizontal resolution (since 360 degrees of longitude were being reflected in relatively few pixels instead of the roughly 10 degrees that should have been fetched).

The solution I developed for this as part of ticket #3179 was to utilize MapServer's special "dateline line/polygon wrapping logic" in msProjectRect(). This logic basically makes it so that small polygons crossing the dateline will be kept local even if the longitudes go outside the -180 to 180 default range. The special dateline wrapping logic was implemented several years ago but was only used when reprojecting shapes. This worked well for this specific problem, but I avoided incorporating it into MapServer 5.6 since I was concerned it might be quite a disruptive change to drop in near release time.

Last week my client contacted me with a related problem. The default view of the world is roughly -21000000,-11000000,21000000,11000000 in web mercator. This is slightly beyond the real extent of the world, so by default PROJ.4 maps -2100000 easting to around +179 longitude even though I might prefer to have thought of it as -181 longitude. The simplest transformation of the web mercator BBOX to lat/long gives an area of roughly -179 to 179 longitude, and -90 to 90 latitude. This is a reasonably close approximation, but loses a bit of data at the edges. This is an issue that has come up in several forms in the past and is sometimes is helped by using of the LOAD_WHOLE_IMAGE processing option.

But my client was running into a different problem. They were seeing WMS requests with a BBOX of -179.79,-89.00,323.96,85.07. The funny item being the longitude of 323. It turns out that this is due to relatively poor handling of world scale polygons in the dateline wrapping logic. It rewraps some points and then realizes that the polygon is really not a small local polygon and it leaves some points unaltered. This can result, as here, in polygons more than 360 degrees wide. This causes some sort of flipping anomaly which I did not investigate.

After some analysis and hair pulling I came to the conclusion that it was essentially impossible to make the polygon dateline wrapping robust in the face of mercator to geographic transformations that had already wrapped things into the -180 to 180 domain. Even if I backed out the approach entirely, I'd still end up with substantial over-requesting in some cases. For this particular problem I did institute a bit of logic in msProjectRect() that said a geographic rect that is wider than 360 degrees should be remapped as minx=-180 and maxx=180. But I was able to construct other requests that would still behave badly.

What I *really* wanted was for PROJ.4 to reproject mercator values of less than -20000000 easting as less than -180 longitude, and more than 20000000 easting as more than 180 longitude. In desperation I finally delved into PROJ.4 to better understand how this wrapping occurs and if their might be some way of avoiding it. I was flabbergasted to discover PROJ.4 was doing the -180/180 wrapping in pj_inv() and that there was a +over switch that could be used to turn this wrapping behavior off!

For my client, I am now advising they use +over in their web mercator projection definition, and this helps avoid many dateline related reprojection problems. I've also added some documentation on this in the PROJ.4 General Parameters documentation. (BTW, a shout out to Hamish who has improved this doc quite a bit!)

"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137
+units=m +over +no_defs"

I'm planning to incorporate my reproject rectangle as a polygon logic in MapServer trunk and I hope people will keep it in mind as a potential source of other problems. But I think wider knowledge of, and use of +over may be the big win of this whole effort! I also hope people will contribute to the World Mapping Issues topic in the MapServer wiki.