RFC 5 - Horizon clipping during reprojection
Frank Warmerdam
warmerdam at POBOX.COM
Fri Oct 21 13:15:48 PDT 2005
On 10/21/05, Bruce Raup <brauplists at gmail.com> wrote:
> Hi Frank,
>
> RFC-5 was just today called to my attention. I think it's great that this
> is being (has been) addressed. Two closely related issues (I think) are these:
>
> A. How are polygon vertices joined? 1) By straight lines in the
> projection, or 2) as great-circle arcs? Can this be specified in the
> mapfile? If not, it would be great to have a LAYER-level parameter,
> something like JOINVERTICES=greatcircle or JOINVERTICES=straightline.
> The GMT program "psxy" offers such a choice (the -A option).
Bruce,
Currently a straight line in the output projection is created between the
point at which the polygon boundary went "over the edge", and where it
came back. This does result in anomolies as we can see a bit in the
sample images I provided.
The reason something smarter (like great circles or some other way
to track the horizon) wasn't done was limited time. The client decided
they didn't want to fund that degree of fidelity unless they find the
current approach unsatisfactory.
> B. When making a selection in the mapserver viewport (via jbox, for
> example), only the corner coordinates of that box get converted back
> to the coordinate system of the data, which means that what is
> selected can be very different from what you see. I wrote about this
> before. See http://spot.colorado.edu/~braup/MapServer/ for an example
> of polar weirdness. Has this issue made it to the bug list? I can do
> that if you guide me on how.
I am aware of this problem, but there has been no effort (other than
a bit of thinking) to resolve it. I'm not sure if there is a bug report on
this issue or not. You can find info at the following url on how to submit
bugs.
http://ms.gis.umn.edu/development/bugs/
The case of reprojecting bounding boxes is handled specially
via the msProjectRect() (in mapproject.c) and this function already
has lots of logic about sampling along edges, and through the rect
if some points fail to reproject. But it does not have any logic
to deal with special polar issues. If we could force it to use the
internal sampling for such cases that would help though not
completely resolve the issue.
I have included the code in case you would like to skim it and
consider possible fixes.
There is a wiki area for discussion of "world mapping problems" where
I liked to accumulate such issues, but with the wiki locked down now
there isn't an easy place to catalog such things.
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
ack. This does result in anomolies as we can see a bit in the
sample images I provided.
The reason something smarter (like great circles or some other way
to track the horizon) wasn't done was limited time. The client decided
they didn't want to fund that degree of fidelity unless they find the
current approach unsatisfactory.
> B. When making a selection in the mapserver viewport (via jbox, for
> example), only the corner coordinates of that box get converted back
> to the coordinate system of the data, which means that what is
> selected can be very different from what you see. I wrote about this
> before. See http://spot.colorado.edu/~braup/MapServer/ for an example
> of polar weirdness. Has this issue made it to the bug list? I can do
> that if you guide me on how.
I am aware of this problem, but there has been no effort (other than
a bit of thinking) to resolve it. I'm not sure if there is a bug report on
this issue or not. You can find info at the following url on how to submit
bugs.
http://ms.gis.umn.edu/development/bugs/
The case of reprojecting bounding boxes is handled specially
via the msProjectRect() (in mapproject.c) and this function already
has lots of logic about sampling along edges, and through the rect
if some points fail to reproject. But it does not have any logic
to deal with special polar issues. If we could force it to use the
internal sampling for such cases that would help though not
completely resolve the issue.
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
/************************************************************************/
/* msProjectGrowRect() */
/************************************************************************/
#ifdef USE_PROJ
static void msProjectGrowRect(projectionObj *in, projectionObj *out,
rectObj *prj_rect, int *rect_initialized,
pointObj *prj_point, int *failure )
{
if( msProjectPoint(in, out, prj_point) == MS_SUCCESS )
{
if( *rect_initialized )
{
prj_rect->miny = MS_MIN(prj_rect->miny, prj_point->y);
prj_rect->maxy = MS_MAX(prj_rect->maxy, prj_point->y);
prj_rect->minx = MS_MIN(prj_rect->minx, prj_point->x);
prj_rect->maxx = MS_MAX(prj_rect->maxx, prj_point->x);
}
else
{
prj_rect->minx = prj_rect->maxx = prj_point->x;
prj_rect->miny = prj_rect->maxy = prj_point->y;
*rect_initialized = MS_TRUE;
}
}
else
(*failure)++;
}
#endif /* def USE_PROJ */
/************************************************************************/
/* msProjectRect() */
/************************************************************************/
#define NUMBER_OF_SAMPLE_POINTS 100
int msProjectRect(projectionObj *in, projectionObj *out, rectObj *rect)
{
#ifdef USE_PROJ
pointObj prj_point;
rectObj prj_rect;
int rect_initialized = MS_FALSE, failure=0;
double dx, dy;
double x, y;
dx = (rect->maxx - rect->minx)/NUMBER_OF_SAMPLE_POINTS;
dy = (rect->maxy - rect->miny)/NUMBER_OF_SAMPLE_POINTS;
/* first ensure the top left corner is processed, even if the rect
turns out to be degenerate. */
prj_point.x = rect->minx;
prj_point.y = rect->miny;
#ifdef USE_POINT_Z_M
prj_point.z = 0.0;
prj_point.m = 0.0;
#endif /* USE_POINT_Z_M */
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&failure);
/* sample along top and bottom */
if(dx > 0) {
for(x=rect->minx; x<=rect->maxx; x+=dx) {
prj_point.x = x;
prj_point.y = rect->miny;
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&failure);
prj_point.x = x;
prj_point.y = rect->maxy;
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&failure);
}
}
/* sample along left and right */
if(dy > 0) {
for(y=rect->miny; y<=rect->maxy; y+=dy) {
prj_point.y = y;
prj_point.x = rect->minx;
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&failure);
prj_point.x = rect->maxx;
prj_point.y = y;
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&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.
*/
if( failure > 0 )
{
failure = 0;
for(x=rect->minx + dx; x<=rect->maxx; x+=dx) {
for(y=rect->miny + dy; y<=rect->maxy; y+=dy) {
prj_point.x = x;
prj_point.y = y;
msProjectGrowRect(in,out,&prj_rect,&rect_initialized,&prj_point,
&failure);
}
}
if( !rect_initialized )
{
if( out == NULL || out->proj == NULL
|| pj_is_latlong(in->proj) )
{
prj_rect.minx = -180;
prj_rect.maxx = 180;
prj_rect.miny = -90;
prj_rect.maxy = 90;
}
else
{
prj_rect.minx = -22000000;
prj_rect.maxx = 22000000;
prj_rect.miny = -11000000;
prj_rect.maxy = 11000000;
}
msDebug( "msProjectRect(): all points failed to reproject,
trying to fall back to using world bounds ... hope this helps.\n" );
}
else
{
msDebug( "msProjectRect(): some points failed to reproject,
doing internal sampling.\n" );
}
}
rect->minx = prj_rect.minx;
rect->miny = prj_rect.miny;
rect->maxx = prj_rect.maxx;
rect->maxy = prj_rect.maxy;
if( !rect_initialized )
return MS_FAILURE;
else
return(MS_SUCCESS);
#else
msSetError(MS_PROJERR, "Projection support is not available.",
"msProjectRect()");
return(MS_FAILURE);
#endif
}
More information about the MapServer-users
mailing list