PROJECTION bug for geographic->geographic

[Cameron Shorter]

Stephen,

I think I've found a bug when you set up a layer projection to geographic and
the map projection to geographic as well.

What happens is the program core dumps ungracefully in the PROJ4 code.
Note that mapserv does not crash if the layer has a projection of utm, or if
the layer has no projection tags at all.

I've added a bit of debugging info.  Hope it helps.

Debugging Info:
===============

Mapserver Version: 3.3.010
Proj4 Version: PROJ.4.3.3.tar.gz
Operating System: Redhat 6.1

Key parts of the mapfile:
NAME npws.map
SIZE 550 370
EXTENT 141.745 -37.970 153.783 -27.785
SHAPEPATH /www/hosted/linuxdev/htdocs/webmap/npws/maps
#UNITS meters
PROJECTION
  geographic
END

...
LAYER
  NAME threatfa
  TYPE Point
  STATUS on
  PROJECTION
   geographic
  END
#PROJECTION
#  "proj=utm"
#END
  DATA threatfa
HEADER 'threatfa.shp_header.html'

    QUERY
        TEMPLATE threatfa.shp_query.html
    END #  QUERY

      CLASS
        Name 'Threatfa'
        COLOR 16 242 77
        OUTLINECOLOR 0 0 0
     END  # CLASS

END  # LAYER

====================

Have attached gdb to mapserv and got the following call stack:

(gdb) r

Starting program: /usr/local/src/ms_3.3.010/mapserv

Program received signal SIGSEGV, Segmentation fault.
pj_fwd (lp={lam = 2.4469395500776421, phi = -0.66270151698224689}, P=0x0)
    at pj_fwd.c:22
22      pj_fwd.c: No such file or directory.

(gdb) bt
#0  pj_fwd (lp={lam = 2.4469395500776421, phi = -0.66270151698224689}, P=0x0)
    at pj_fwd.c:22
#1  0x8069749 in msProjectPoint (in=0x0, out=0x0, point=0xbfffe428)
    at mapproject.c:24
#2  0x8069784 in msProjectRect (in=0x0, out=0x0, rect=0xbfffe45c)
    at mapproject.c:40
#3  0x8054914 in msWhichShapesProj (shp=0xbfffee04, window={
      minx = 140.19930894308945, miny = -37.969999999999999,
      maxx = 155.32869105691054, maxy = -27.785}, in=0x80bfda0,
out=0x4020e834)
    at mapsearch.c:332
#4  0x805ed88 in msDrawShapefileLayer (map=0x401e4008, layer=0x80bfd00,
    img=0x80f12a0, query_status=0x0) at maputil.c:1136
#5  0x805c7e1 in msDrawMap (map=0x401e4008) at maputil.c:401
#6  0x8050002 in main (argc=1, argv=0xbffff6c4) at mapserv.c:2160
(gdb)
(gdb) p P
$1 = (PJ *) 0x0

==================
pj_fwd.c contains:
/* general forward projection */
#ifndef lint
static const char SCCSID[]="@(#)pj_fwd.c        4.4     93/06/12        GIE
REL";
#endif
#define PJ_LIB__
#include <projects.h>
#include <errno.h>
# define EPS 1.0e-12
        XY /* forward projection entry */
pj_fwd(LP lp, PJ *P) {
        XY xy;
        double t;

        /* check for forward and latitude or longitude overange */
        if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 10.) {
                xy.x = xy.y = HUGE_VAL;
                pj_errno = -14;
        } else { /* proceed with projection */
                errno = pj_errno = 0;
                if (fabs(t) <= EPS)
                        lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
<22>            else if (P->geoc)
                        lp.phi = atan(P->rone_es * tan(lp.phi));
                lp.lam -= P->lam0;      /* compute del lp.lam */
                if (!P->over)
                        lp.lam = adjlon(lp.lam); /* adjust del longitude */
                xy = (*P->fwd)(lp, P); /* project */
                if (pj_errno || (pj_errno = errno))
                        xy.x = xy.y = HUGE_VAL;
                /* adjust for major axis and easting/northings */
                else {
                        xy.x = P->fr_meter * (P->a * xy.x + P->x0);
                        xy.y = P->fr_meter * (P->a * xy.y + P->y0);
                }
        }
        return xy;
}

Note on line 22 (where the error occurs) that P=NULL.
My guess is that there is a logic problem in mapproject.c

-- 
Cameron Shorter          Web Mapping Manager
Social Change Online
6A Nelson Street         Tel: +61 (0) 2 9557 6500
Annandale NSW 2038       Fax: +61 (0) 2 9519 8940
Sydney, Australia        http://online.socialchange.net.au