Skip to content
Permalink
Browse files
Initial port to PROJ 6 API
  • Loading branch information
rouault committed Oct 1, 2019
1 parent 7506203 commit 680e7e0edfbd05caf1ed5812e008613b6dc484d0
Show file tree
Hide file tree
Showing 15 changed files with 896 additions and 371 deletions.
@@ -10,8 +10,26 @@ FIND_PATH(PROJ_INCLUDE_DIR proj_api.h)

FIND_LIBRARY(PROJ_LIBRARY NAMES proj proj_i)

set(PROJ_INCLUDE_DIRS ${PROJ_INCLUDE_DIR})
set(PROJ_LIBRARIES ${PROJ_LIBRARY})
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(PROJ DEFAULT_MSG PROJ_LIBRARY PROJ_INCLUDE_DIR)
mark_as_advanced(PROJ_LIBRARY PROJ_INCLUDE_DIR)


IF (PROJ_INCLUDE_DIR AND PROJ_LIBRARY)
SET(PROJ_FOUND TRUE)
ENDIF (PROJ_INCLUDE_DIR AND PROJ_LIBRARY)

IF (PROJ_FOUND)
IF (EXISTS ${PROJ_INCLUDE_DIR}/proj.h)
FILE(READ ${PROJ_INCLUDE_DIR}/proj.h proj_version)
STRING(REGEX REPLACE "^.*PROJ_VERSION_MAJOR +([0-9]+).*$" "\\1" PROJ_VERSION_MAJOR "${proj_version}")
STRING(REGEX REPLACE "^.*PROJ_VERSION_MINOR +([0-9]+).*$" "\\1" PROJ_VERSION_MINOR "${proj_version}")
STRING(REGEX REPLACE "^.*PROJ_VERSION_PATCH +([0-9]+).*$" "\\1" PROJ_VERSION_PATCH "${proj_version}")

MESSAGE(STATUS "Found Proj ${PROJ_VERSION_MAJOR}.${PROJ_VERSION_MINOR}")
ADD_DEFINITIONS(-DPROJ_VERSION_MAJOR=${PROJ_VERSION_MAJOR})
ELSE()
MESSAGE(STATUS "Found Proj 4.x")
ADD_DEFINITIONS(-DPROJ_VERSION_MAJOR=4)
ENDIF()
ENDIF (PROJ_FOUND)
@@ -985,8 +985,8 @@ int msDrawVectorLayer(mapObj *map, layerObj *layer, imageObj *image)
/* the layer projection. */
if( layer->connectiontype == MS_UVRASTER &&
!layer->projection.gt.need_geotransform &&
!(pj_is_latlong(map->projection.proj) &&
pj_is_latlong(layer->projection.proj)) ) {
!(msProjIsGeographicCRS(&(map->projection)) &&
msProjIsGeographicCRS(&(layer->projection))) ) {
rectObj layer_ori_extent;

if( msLayerGetExtent(layer, &layer_ori_extent) == MS_SUCCESS ) {
234 mapfile.c
@@ -1113,240 +1113,6 @@ static void writeGrid(FILE *stream, int indent, graticuleObj *pGraticule)
writeBlockEnd(stream, indent, "GRID");
}

/*
** Initialize, load and free a projectionObj structure
*/
int msInitProjection(projectionObj *p)
{
p->gt.need_geotransform = MS_FALSE;
p->numargs = 0;
p->args = NULL;
p->wellknownprojection = wkp_none;
#ifdef USE_PROJ
p->proj = NULL;
p->args = (char **)malloc(MS_MAXPROJARGS*sizeof(char *));
MS_CHECK_ALLOC(p->args, MS_MAXPROJARGS*sizeof(char *), -1);
#if PJ_VERSION >= 480
p->proj_ctx = NULL;
#endif
#endif
return(0);
}

void msFreeProjection(projectionObj *p)
{
#ifdef USE_PROJ
if(p->proj) {
pj_free(p->proj);
p->proj = NULL;
}
#if PJ_VERSION >= 480
if(p->proj_ctx) {
pj_ctx_free(p->proj_ctx);
p->proj_ctx = NULL;
}
#endif

msFreeCharArray(p->args, p->numargs);
p->args = NULL;
p->numargs = 0;
#endif
}

/*
** Handle OGC WMS/WFS AUTO projection in the format:
** "AUTO:proj_id,units_id,lon0,lat0"
*/
#ifdef USE_PROJ
static int _msProcessAutoProjection(projectionObj *p)
{
char **args;
int numargs, nProjId, nUnitsId, nZone;
double dLat0, dLon0;
const char *pszUnits = "m";
char szProjBuf[512]="";

/* WMS/WFS AUTO projection: "AUTO:proj_id,units_id,lon0,lat0" */
args = msStringSplit(p->args[0], ',', &numargs);
if (numargs != 4 ||
(strncasecmp(args[0], "AUTO:", 5) != 0 &&
strncasecmp(args[0], "AUTO2:", 6) != 0)) {
msSetError(MS_PROJERR,
"WMS/WFS AUTO/AUTO2 PROJECTION must be in the format "
"'AUTO:proj_id,units_id,lon0,lat0' or 'AUTO2:crs_id,factor,lon0,lat0'(got '%s').\n",
"_msProcessAutoProjection()", p->args[0]);
return -1;
}

if (strncasecmp(args[0], "AUTO:", 5)==0)
nProjId = atoi(args[0]+5);
else
nProjId = atoi(args[0]+6);

nUnitsId = atoi(args[1]);
dLon0 = atof(args[2]);
dLat0 = atof(args[3]);


/*There is no unit parameter for AUTO2. The 2nd parameter is
factor. Set the units to always be meter*/
if (strncasecmp(args[0], "AUTO2:", 6) == 0)
nUnitsId = 9001;

msFreeCharArray(args, numargs);

/* Handle EPSG Units. Only meters for now. */
switch(nUnitsId) {
case 9001: /* Meters */
pszUnits = "m";
break;
default:
msSetError(MS_PROJERR,
"WMS/WFS AUTO PROJECTION: EPSG Units %d not supported.\n",
"_msProcessAutoProjection()", nUnitsId);
return -1;
}

/* Build PROJ4 definition.
* This is based on the definitions found in annex E of the WMS 1.1.1
* spec and online at http://www.digitalearth.org/wmt/auto.html
* The conversion from the original WKT definitions to PROJ4 format was
* done using the MapScript setWKTProjection() function (based on OGR).
*/
switch(nProjId) {
case 42001: /** WGS 84 / Auto UTM **/
nZone = (int) floor( (dLon0 + 180.0) / 6.0 ) + 1;
sprintf( szProjBuf,
"+proj=tmerc+lat_0=0+lon_0=%.16g+k=0.999600+x_0=500000"
"+y_0=%.16g+ellps=WGS84+datum=WGS84+units=%s",
-183.0 + nZone * 6.0,
(dLat0 >= 0.0) ? 0.0 : 10000000.0,
pszUnits);
break;
case 42002: /** WGS 84 / Auto Tr. Mercator **/
sprintf( szProjBuf,
"+proj=tmerc+lat_0=0+lon_0=%.16g+k=0.999600+x_0=500000"
"+y_0=%.16g+ellps=WGS84+datum=WGS84+units=%s",
dLon0,
(dLat0 >= 0.0) ? 0.0 : 10000000.0,
pszUnits);
break;
case 42003: /** WGS 84 / Auto Orthographic **/
sprintf( szProjBuf,
"+proj=ortho+lon_0=%.16g+lat_0=%.16g+x_0=0+y_0=0"
"+ellps=WGS84+datum=WGS84+units=%s",
dLon0, dLat0, pszUnits );
break;
case 42004: /** WGS 84 / Auto Equirectangular **/
/* Note that we have to pass lon_0 as lon_ts for this one to */
/* work. Either a PROJ4 bug or a PROJ4 documentation issue. */
sprintf( szProjBuf,
"+proj=eqc+lon_ts=%.16g+lat_ts=%.16g+x_0=0+y_0=0"
"+ellps=WGS84+datum=WGS84+units=%s",
dLon0, dLat0, pszUnits);
break;
case 42005: /** WGS 84 / Auto Mollweide **/
sprintf( szProjBuf,
"+proj=moll+lon_0=%.16g+x_0=0+y_0=0+ellps=WGS84"
"+datum=WGS84+units=%s",
dLon0, pszUnits);
break;
default:
msSetError(MS_PROJERR,
"WMS/WFS AUTO PROJECTION %d not supported.\n",
"_msProcessAutoProjection()", nProjId);
return -1;
}

/* msDebug("%s = %s\n", p->args[0], szProjBuf); */

/* OK, pass the definition to pj_init() */
args = msStringSplit(szProjBuf, '+', &numargs);

msAcquireLock( TLOCK_PROJ );
if( !(p->proj = pj_init(numargs, args)) ) {
int *pj_errno_ref = pj_get_errno_ref();
msReleaseLock( TLOCK_PROJ );
msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
"msProcessProjection()", pj_strerrno(*pj_errno_ref), szProjBuf) ;
return(-1);
}

msReleaseLock( TLOCK_PROJ );

msFreeCharArray(args, numargs);

return(0);
}
#endif /* USE_PROJ */

int msProcessProjection(projectionObj *p)
{
#ifdef USE_PROJ
assert( p->proj == NULL );

if( strcasecmp(p->args[0],"GEOGRAPHIC") == 0 ) {
msSetError(MS_PROJERR,
"PROJECTION 'GEOGRAPHIC' no longer supported.\n"
"Provide explicit definition.\n"
"ie. proj=latlong\n"
" ellps=clrk66\n",
"msProcessProjection()");
return(-1);
}

if (strcasecmp(p->args[0], "AUTO") == 0) {
p->proj = NULL;
return 0;
}

if (strncasecmp(p->args[0], "AUTO:", 5) == 0 ||
strncasecmp(p->args[0], "AUTO2:", 6) == 0) {
/* WMS/WFS AUTO projection: "AUTO:proj_id,units_id,lon0,lat0" */
/*WMS 1.3.0: AUTO2:auto_crs_id,factor,lon0,lat0*/
return _msProcessAutoProjection(p);
}
msAcquireLock( TLOCK_PROJ );
#if PJ_VERSION < 480
if( !(p->proj = pj_init(p->numargs, p->args)) ) {
#else
p->proj_ctx = pj_ctx_alloc();
if( !(p->proj=pj_init_ctx(p->proj_ctx, p->numargs, p->args)) ) {
#endif

int *pj_errno_ref = pj_get_errno_ref();
msReleaseLock( TLOCK_PROJ );
if(p->numargs>1) {
msSetError(MS_PROJERR, "proj error \"%s\" for \"%s:%s\"",
"msProcessProjection()", pj_strerrno(*pj_errno_ref), p->args[0],p->args[1]) ;
} else {
msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
"msProcessProjection()", pj_strerrno(*pj_errno_ref), p->args[0]) ;
}
return(-1);
}

msReleaseLock( TLOCK_PROJ );

#ifdef USE_PROJ_FASTPATHS
if(strcasestr(p->args[0],"epsg:4326")) {
p->wellknownprojection = wkp_lonlat;
} else if(strcasestr(p->args[0],"epsg:3857")) {
p->wellknownprojection = wkp_gmerc;
} else {
p->wellknownprojection = wkp_none;
}
#endif


return(0);
#else
msSetError(MS_PROJERR, "Projection support is not available.",
"msProcessProjection()");
return(-1);
#endif
}

static int loadProjection(projectionObj *p)
{
#ifdef USE_PROJ
@@ -191,7 +191,7 @@ int msGraticuleLayerWhichShapes(layerObj *layer, rectObj rect, int isQuery)
layer->project = msProjectionsDiffer(&(layer->projection), &(layer->map->projection));
if( layer->project &&
strstr(layer->map->projection.args[0], "epsg:3857") &&
pj_is_latlong(layer->projection.proj) )
msProjIsGeographicCRS(&(layer->projection)) )
{
if( rectMapCoordinates.minx < -20037508)
rectMapCoordinates.minx = -20037508;
@@ -1073,8 +1073,8 @@ static int _AdjustLabelPosition( layerObj *pLayer, shapeObj *pShape, msGraticule
msProjectShape( &pLayer->projection, &pLayer->map->projection, pShape );

/* Poor man detection of reprojection failure */
if( pj_is_latlong(pLayer->projection.proj) !=
pj_is_latlong(pLayer->map->projection.proj) )
if( msProjIsGeographicCRS(&(pLayer->projection)) !=
msProjIsGeographicCRS(&(pLayer->map->projection)) )
{
if( ptPoint.x == pShape->line->point[0].x &&
ptPoint.y == pShape->line->point[0].y )
@@ -1128,7 +1128,7 @@ static int _AdjustLabelPosition( layerObj *pLayer, shapeObj *pShape, msGraticule
/* Clamp coordinates into the validity area of the projection, in the */
/* particular case of EPSG:3857 (WebMercator) to longlat reprojection */
if( strstr(pLayer->map->projection.args[0], "epsg:3857") &&
pj_is_latlong(pLayer->projection.proj) )
msProjIsGeographicCRS(&(pLayer->projection)) )
{
if( !pLayer->map->projection.gt.need_geotransform &&
ePosition == posLeft && pShape->line->point[0].x < -20037508)
@@ -523,7 +523,7 @@ int KmlRenderer::checkProjection(mapObj *map)
{
projectionObj *projection= &map->projection;
#ifdef USE_PROJ
if (projection && projection->numargs > 0 && pj_is_latlong(projection->proj)) {
if (projection && projection->numargs > 0 && msProjIsGeographicCRS(projection)) {
return MS_SUCCESS;
} else {
char epsg_string[100];
@@ -2090,7 +2090,7 @@ int msOWSPrintEncodeParamList(FILE *stream, const char *name,
*/
void msOWSProjectToWGS84(projectionObj *srcproj, rectObj *ext)
{
if (srcproj->numargs > 0 && !pj_is_latlong(srcproj->proj)) {
if (srcproj->proj && !msProjIsGeographicCRS(srcproj)) {
projectionObj wgs84;
msInitProjection(&wgs84);
msLoadProjectionString(&wgs84, "+proj=longlat +ellps=WGS84 +datum=WGS84");

0 comments on commit 680e7e0

Please sign in to comment.