Skip to content
Permalink
Browse files
msProjectShapeLine(): take into account antimeridian when reprojectin…
…g lines from Polar Stereographic to geographic / WebMercator
  • Loading branch information
rouault committed Jan 15, 2020
1 parent 44dde80 commit 9ad9824c2f6f3b8b12e9ab58f0e7144c3e551aba
Show file tree
Hide file tree
Showing 2 changed files with 234 additions and 0 deletions.
@@ -66,6 +66,8 @@ struct reprojectionObj
projectionObj* in;
projectionObj* out;
PJ* pj;
int should_do_line_cutting;
shapeObj splitShape;
int bFreePJ;
};

@@ -366,6 +368,7 @@ reprojectionObj* msProjectCreateReprojector(projectionObj* in, projectionObj* ou
reprojectionObj* obj = (reprojectionObj*)msSmallCalloc(1, sizeof(reprojectionObj));
obj->in = in;
obj->out = out;
obj->should_do_line_cutting = -1;

/* -------------------------------------------------------------------- */
/* If the source and destination are simple and equal, then do */
@@ -448,6 +451,7 @@ void msProjectDestroyReprojector(reprojectionObj* reprojector)
return;
if( reprojector->bFreePJ )
proj_destroy(reprojector->pj);
msFreeShape(&(reprojector->splitShape));
msFree(reprojector);
}

@@ -490,6 +494,8 @@ struct reprojectionObj
{
projectionObj* in;
projectionObj* out;
int should_do_line_cutting;
shapeObj splitShape;
int no_op;
};

@@ -503,6 +509,7 @@ reprojectionObj* msProjectCreateReprojector(projectionObj* in, projectionObj* ou
obj = (reprojectionObj*)msSmallCalloc(1, sizeof(reprojectionObj));
obj->in = in;
obj->out = out;
obj->should_do_line_cutting = -1;

/* -------------------------------------------------------------------- */
/* If the source and destination are equal, then do nothing. */
@@ -548,6 +555,7 @@ void msProjectDestroyReprojector(reprojectionObj* reprojector)
{
if( !reprojector )
return;
msFreeShape(&(reprojector->splitShape));
msFree(reprojector);
}

@@ -1134,6 +1142,96 @@ static int msProjectSegment( reprojectionObj* reprojector,
return MS_SUCCESS;
}

/************************************************************************/
/* msProjectShapeShouldDoLineCutting() */
/************************************************************************/

/* Detect projecting from north polar stereographic to longlat or EPSG:3857 */
#ifdef USE_GEOS
static int msProjectShapeShouldDoLineCutting(reprojectionObj* reprojector)
{
if( reprojector->should_do_line_cutting >= 0)
return reprojector->should_do_line_cutting;

projectionObj *in = reprojector->in;
projectionObj *out = reprojector->out;
if( !(!in->gt.need_geotransform && !msProjIsGeographicCRS(in) &&
(msProjIsGeographicCRS(out) ||
(out->numargs == 1 && strcmp(out->args[0], "init=epsg:3857") == 0))) )
{
reprojector->should_do_line_cutting = MS_FALSE;
return MS_FALSE;
}

int srcIsPolar;
double extremeLongEasting;
if( msProjIsGeographicCRS(out) )
{
pointObj p;
double gt3 = out->gt.need_geotransform ? out->gt.geotransform[3] : 0.0;
double gt4 = out->gt.need_geotransform ? out->gt.geotransform[4] : 0.0;
double gt5 = out->gt.need_geotransform ? out->gt.geotransform[5] : 1.0;
p.x = 0.0;
p.y = 0.0;
srcIsPolar = msProjectPointEx(reprojector, &p) == MS_SUCCESS &&
fabs(gt3 + p.x * gt4 + p.y * gt5 - 90) < 1e-8;
extremeLongEasting = 180;
}
else
{
pointObj p1;
pointObj p2;
double gt1 = out->gt.need_geotransform ? out->gt.geotransform[1] : 1.0;
p1.x = 0.0;
p1.y = -0.1;
p2.x = 0.0;
p2.y = 0.1;
srcIsPolar = msProjectPointEx(reprojector, &p1) == MS_SUCCESS &&
msProjectPointEx(reprojector, &p2) == MS_SUCCESS &&
fabs((p1.x - p2.x) * gt1) > 20e6;
extremeLongEasting = 20037508.3427892;
}
if( !srcIsPolar )
{
reprojector->should_do_line_cutting = MS_FALSE;
return MS_FALSE;
}

pointObj p;
double invgt0 = out->gt.need_geotransform ? out->gt.invgeotransform[0] : 0.0;
double invgt1 = out->gt.need_geotransform ? out->gt.invgeotransform[1] : 1.0;
double invgt3 = out->gt.need_geotransform ? out->gt.invgeotransform[3] : 0.0;
double invgt4 = out->gt.need_geotransform ? out->gt.invgeotransform[4] : 0.0;

lineObj newLine = {0,NULL};
const double EPS = 1e-10;

p.x = invgt0 + -extremeLongEasting * (1 - EPS) * invgt1;
p.y = invgt3 + -extremeLongEasting * (1 - EPS) * invgt4;
msProjectPoint(out, in, &p);
pointObj first = p;
msAddPointToLine(&newLine, &p );

p.x = invgt0 + extremeLongEasting * (1 - EPS) * invgt1;
p.y = invgt3 + extremeLongEasting * (1 - EPS) * invgt4;
msProjectPoint(out, in, &p);
msAddPointToLine(&newLine, &p );

p.x = 0;
p.y = 0;
msAddPointToLine(&newLine, &p );

msAddPointToLine(&newLine, &first );

msInitShape(&(reprojector->splitShape));
reprojector->splitShape.type = MS_SHAPE_POLYGON;
msAddLineDirectly(&(reprojector->splitShape), &newLine);

reprojector->should_do_line_cutting = MS_TRUE;
return MS_TRUE;
}
#endif

/************************************************************************/
/* msProjectShapeLine() */
/* */
@@ -1197,7 +1295,49 @@ msProjectShapeLine(reprojectionObj* reprojector,
#undef p_y
#endif

#ifdef USE_GEOS
if( shape->type == MS_SHAPE_LINE &&
msProjectShapeShouldDoLineCutting(reprojector) )
{
shapeObj tmpShapeInputLine;
msInitShape(&tmpShapeInputLine);
tmpShapeInputLine.type = MS_SHAPE_LINE;
tmpShapeInputLine.numlines = 1;
tmpShapeInputLine.line = line;

shapeObj* diff = NULL;
if( msGEOSIntersects(&tmpShapeInputLine, &(reprojector->splitShape)) )
{
diff = msGEOSDifference(&tmpShapeInputLine, &(reprojector->splitShape));
}

tmpShapeInputLine.numlines = 0;
tmpShapeInputLine.line = NULL;
msFreeShape(&tmpShapeInputLine);

if( diff )
{
for(int j = 0; j < diff->numlines; j++ )
{
for( i=0; i < diff->line[j].numpoints; i++ ) {
msProjectPointEx(reprojector, &(diff->line[j].point[i]));
}
if( j == 0 )
{
line_out->numpoints = diff->line[j].numpoints;
memcpy(line_out->point, diff->line[0].point, sizeof(pointObj) * line_out->numpoints);
}
else
{
msAddLineDirectly(shape, &(diff->line[j]));
}
}
msFreeShape(diff);
msFree(diff);
return MS_SUCCESS;
}
}
#endif

wrap_test = out != NULL && out->proj != NULL && msProjIsGeographicCRS(out)
&& !msProjIsGeographicCRS(in);
@@ -123,9 +123,103 @@ def test_pattern():

return 'success'

###############################################################################
# Test reprojection of lines from Polar Stereographic and crossing the antimerdian

def test_reprojection_lines_from_polar_stereographic_to_geographic():

shape = mapscript.shapeObj( mapscript.MS_SHAPE_LINE )
line = mapscript.lineObj()
line.add( mapscript.pointObj( -5554682.77568025, -96957.3485051449 ) ) # -179 40
line.add( mapscript.pointObj( -5554682.77568025, 96957.3485051449 ) ) # 179 40
shape.add( line )

polar_proj = mapscript.projectionObj("+proj=stere +lat_0=90 +lat_ts=60 +lon_0=270 +datum=WGS84")
longlat_proj = mapscript.projectionObj("+proj=longlat +datum=WGS84")

if shape.project(polar_proj, longlat_proj) != 0:
pmstestlib.post_reason('shape.project() failed')
return 'fail'

part1 = shape.get(0)
part2 = shape.get(1)
if not part1 or not part2:
pmstestlib.post_reason('should get two parts')
return 'fail'

point11 = part1.get(0)
point12 = part1.get(1)
point21 = part2.get(0)
point22 = part2.get(1)
if abs(point11.x - -179.0) > 1e-7:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point12.x - -180.0) > 1e-7:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point21.x - 180.0) > 1e-7:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point22.x - 179.0) > 1e-7:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
return 'success'

###############################################################################
# Test reprojection of lines from Polar Stereographic and crossing the antimerdian

def test_reprojection_lines_from_polar_stereographic_to_webmercator():

shape = mapscript.shapeObj( mapscript.MS_SHAPE_LINE )
line = mapscript.lineObj()
line.add( mapscript.pointObj( -5554682.77568025, -96957.3485051449 ) ) # -179 40
line.add( mapscript.pointObj( -5554682.77568025, 96957.3485051449 ) ) # 179 40
shape.add( line )

polar_proj = mapscript.projectionObj("+proj=stere +lat_0=90 +lat_ts=60 +lon_0=270 +datum=WGS84")
longlat_proj = mapscript.projectionObj("init=epsg:3857")

if shape.project(polar_proj, longlat_proj) != 0:
pmstestlib.post_reason('shape.project() failed')
return 'fail'

part1 = shape.get(0)
part2 = shape.get(1)
if not part1 or not part2:
pmstestlib.post_reason('should get two parts')
return 'fail'

point11 = part1.get(0)
point12 = part1.get(1)
point21 = part2.get(0)
point22 = part2.get(1)
if abs(point11.x - -19926188.85) > 1e-2:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point12.x - -20037508.34) > 1e-2:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point21.x - 20037508.34) > 1e-2:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
if abs(point22.x - 19926188.85) > 1e-2:
print(point11.x, point12.x, point21.x, point22.x)
pmstestlib.post_reason('did not get expected coordinates')
return 'fail'
return 'success'

test_list = [
bug_673,
test_pattern,
test_reprojection_lines_from_polar_stereographic_to_geographic,
test_reprojection_lines_from_polar_stereographic_to_webmercator,
None ]

if __name__ == '__main__':

0 comments on commit 9ad9824

Please sign in to comment.