[Mapserver-dev] three reprojection bugs
Dan Stahlke
dstahlke at gi.alaska.edu
Tue Jun 8 22:31:58 EDT 2004
--=-2+S175RXjQo4aqeBe52s
Content-Type: text/plain
Content-Transfer-Encoding: 7bit
There are three bugs that show up when using an orthographic projection
(other projections are affected too) in WMS-client mode with a source
in geographic projection:
1. If the projection is valid along the borders, but contains the
north pole, the source bounding box will not be correct. In
these cases the whole grid needs to be sampled in order to
determine the source bounds. Perhaps this should be a
processing parameter option?
2. If there are invalid points on a particular scanline, the whole
scanline fails. The problem is that pj_transform returns a
failure and the whole scanline is discarded by mapserver. It
seems like pj_transform should transform valid points and
set invalid ones to HUGE_VAL, but this is not the case. Maybe
this is a bug in proj4... the manpage doesn't really say what
the behavior should be in this case.
3. If the destination image contains the -180/180 border, there
will be a stripe of invalid data there. The pixels that
are near 180 degrees longitude will correspond to points
in the source image that are up to 0.5 pixels past the
right-hand edge of the image. It seems to me that somewhere
in GDALGetGeoTransform there is a 0.5 that is inappropriately
added. A quick fix is to allow a one-pixel margin of error
in msSimpleRasterResampler.
Attached to this message is an example map file and perl script
that will demonstrate the problems. The zoomed-in image will
have a blank area near the north pole and the zoomed-out image
will only draw the center third: the top and bottom will be
blank. There is also a quick-hack patch that shows a workaround
for these problems.
As a side note, it would be nice to have an option that causes
a higher-res source image to be downloaded in order to avoid
the oversampling effect that sometimes happens (the demo script
included with this message produces a noticably pixelated image).
- Dan Stahlke
--=-2+S175RXjQo4aqeBe52s
Content-Disposition: attachment; filename=cascade.map
Content-Type: text/plain; name=cascade.map; charset=utf-8
Content-Transfer-Encoding: 7bit
MAP
NAME CASCADE
STATUS ON
DEBUG ON
IMAGECOLOR 255 255 255
IMAGETYPE PNG24
WEB
LOG "/home/dstahlke/mapserver/cascade/tmp/log.txt"
IMAGEPATH "/home/dstahlke/mapserver/cascade/tmp/"
END
LAYER
NAME GLCF
TYPE RASTER
STATUS ON
DEBUG ON
# these don't take effect until after the source is downloaded...
PROCESSING "LOAD_WHOLE_IMAGE=YES"
PROCESSING "LOAD_FULL_RES_IMAGE=YES"
CONNECTIONTYPE WMS
CONNECTION "http://onearth.jpl.nasa.gov/gis-wms.cgi?"
METADATA
wms_srs "EPSG:4326"
wms_name "global_mosaic"
wms_server_version "1.1.1"
wms_format "image/jpeg"
END
END
END
--=-2+S175RXjQo4aqeBe52s
Content-Disposition: attachment; filename=test.pl
Content-Type: application/x-perl; name=test.pl
Content-Transfer-Encoding: 7bit
#!/usr/bin/perl -w
use warnings;
use strict;
use mapscript;
gen_map('zoom_in.png', 4000000);
gen_map('zoom_out.png', 6000000);
sub gen_map {
my($out_fn, $range) = @_;
unlink $out_fn if -e $out_fn;
my $map = new mapscript::mapObj('cascade.map') or die('Unable to open mapfile.');
#$map->setProjection('proj=qsc,face=0,a=4000000');
$map->setProjection('proj=ortho,lon_0=-147.84944,lat_0=64.85944');
my $rect= new mapscript::rectObj();
$rect->{minx} = $rect->{miny} = -$range;
$rect->{maxx} = $rect->{maxy} = $range;
$map->{extent} = $rect;
$map->{width} = 800;
$map->{height} = 800;
$map->{debug} = 1;
my $img = $map->draw() or die('Unable to draw map');
$img->save($out_fn);
-e $out_fn or die 'Unable to save map';
}
--=-2+S175RXjQo4aqeBe52s
Content-Disposition: attachment; filename=quickfix.diff
Content-Type: text/x-patch; name=quickfix.diff; charset=utf-8
Content-Transfer-Encoding: 7bit
Common subdirectories: ms420-orig/fonts and ms420-mod/fonts
diff -u ms420-orig/mapproject.c ms420-mod/mapproject.c
--- ms420-orig/mapproject.c 2004-03-15 20:10:00.000000000 -0900
+++ ms420-mod/mapproject.c 2004-06-08 17:15:55.334247168 -0800
@@ -135,6 +135,11 @@
}
}
+// for some projections, sampling around the edges isn't
+// sufficient - for example, orthographic with north pole
+// at the center
+failure++;
+
/*
** If there have been any failures around the edges, then we had better
** try and fill in the interior to get a close bounds.
@@ -151,10 +156,31 @@
}
}
+// even sampling a grid doesn't get all the way to the
+// north pole... so expand by 10%
+ if( rect_initialized )
+ {
+ dx = prj_rect.maxx - prj_rect.minx;
+ dy = prj_rect.maxy - prj_rect.miny;
+ prj_rect.minx -= dx * .1;
+ prj_rect.maxx += dx * .1;
+ prj_rect.miny -= dy * .1;
+ prj_rect.maxy += dy * .1;
+
+ if( out == NULL || out->proj == NULL
+ || pj_is_latlong(out->proj) )
+ {
+ if(prj_rect.minx < -180) prj_rect.minx = -180;
+ if(prj_rect.maxx > 180) prj_rect.maxx = 180;
+ if(prj_rect.miny < -90) prj_rect.miny = -90;
+ if(prj_rect.maxy > 90) prj_rect.maxy = 90;
+ }
+ }
+
if( !rect_initialized )
{
if( out == NULL || out->proj == NULL
- || pj_is_latlong(in->proj) )
+ || pj_is_latlong(out->proj) ) // MOD - bug?
{
prj_rect.minx = -180;
prj_rect.maxx = 180;
diff -u ms420-orig/mapresample.c ms420-mod/mapresample.c
--- ms420-orig/mapresample.c 2004-05-21 09:50:17.000000000 -0800
+++ ms420-mod/mapresample.c 2004-06-08 18:14:08.273239728 -0800
@@ -296,6 +296,9 @@
nSrcX = (int) x[nDstX];
nSrcY = (int) y[nDstX];
+ // workaround for date-line crossover
+ if(nSrcX == nSrcXSize) nSrcX = nSrcXSize-1;
+
/*
* We test the original floating point values to
* avoid errors related to asymmetric rounding around zero.
@@ -552,19 +555,19 @@
if( psPTInfo->psDstProj != NULL
&& psPTInfo->psSrcProj != NULL )
{
- double *z;
-
- z = (double *) calloc(sizeof(double),nPoints);
- if( pj_transform( psPTInfo->psDstProj, psPTInfo->psSrcProj,
- nPoints, 1, x, y, z) != 0 )
- {
- free( z );
- for( i = 0; i < nPoints; i++ )
+// transform points one at a time - fixes problem where scan
+// lines are not drawn when they contain invalid points,
+// such as orthographic projection where viewpoint is zoomed
+// out to the point where the globe doesn't cover the whole screen
+ double z;
+ for( i = 0; i < nPoints; i++ ) {
+ z = 0;
+ if( pj_transform( psPTInfo->psDstProj, psPTInfo->psSrcProj,
+ 1, 1, x+i, y+i, &z) != 0 )
+ {
panSuccess[i] = 0;
-
- return MS_FALSE;
- }
- free( z );
+ }
+ }
}
for( i = 0; i < nPoints; i++ )
Common subdirectories: ms420-orig/mapscript and ms420-mod/mapscript
Common subdirectories: ms420-orig/symbols and ms420-mod/symbols
Common subdirectories: ms420-orig/tests and ms420-mod/tests
--=-2+S175RXjQo4aqeBe52s--
More information about the mapserver-dev
mailing list