Reprojection for Lat/Long Wrong

[warmerdam]

I have a lat/long layer with the following definition:

LAYER
  NAME GTOPO
  TYPE raster
  STATUS ON
  DATA /usr2/data/jpeg2000/gtopo.jp2
  DUMP TRUE
#  OFFSITE 255 255 255
  PROJECTION
    "init=epsg:4326"
  END
  METADATA 
    "wcs_title" "GTOPO"
    "wcs_label" "GTOPO - label"
    "wcs_extent" "-100 -10 -60 40"
    "wcs_resolution" "0.0083333 0.0083333"
  END
END

When I get the capabilities from the server, it returns the wrong latlong
region:

<ContentMetadata>
  <CoverageOfferingBrief>
    <name>GTOPO</name>
    <label>GTOPO - label</label>
    <lonLatEnvelope srsName="WGS84(DD)">
      <gml:pos>-178.682 -572.958</gml:pos>
      <gml:pos>178.479 2291.83</gml:pos>
    </lonLatEnvelope>
  </CoverageOfferingBrief>
</ContentMetadata>

The problem seems to be that the following code at the end of 
msWCSGetCoverageMetadata() calls msProjectRect() with a latlong projection
as input and NULL as output but down in msProjectPoint() it doesn't know
how to reproject to/from a latlong coordinate system. 

    msInitProjection(&proj); // or bad things happen

    sprintf(projstring, "init=epsg:%.20s", cm->srs+5);
    if (msLoadProjectionString(&proj, projstring) != 0) return MS_FAILURE;
    msProjectRect(&proj, NULL, &(cm->llextent));

The mode in msProjectPoint() where the input or output coordinate system
can be NULL to mean lat/long does not seem to work properly if the other
coordinate system is also lat/long.  In particular the following logic 
fails to apply the degree to radian conversion on the input but it is
applied on the output.

  else
  {
      p.u = point->x;
      p.v = point->y;

      if(in==NULL || in->proj==NULL) { /* input coordinates are lat/lon */
          p.u *= DEG_TO_RAD; /* convert to radians */
          p.v *= DEG_TO_RAD;  
          p = pj_fwd(p, out->proj);
      } else {
          if(out==NULL || out->proj==NULL) { /* output coordinates are lat/lon */
              p = pj_inv(p, in->proj);
              p.u *= RAD_TO_DEG; /* convert to decimal degrees */
              p.v *= RAD_TO_DEG;
          } else { /* need to go from one projection to another */
              p = pj_inv(p, in->proj);
              p = pj_fwd(p, out->proj);
          }
      }

      point->x = p.u;
      point->y = p.v;
  }

[fwarmerdam]

As a first step, I have corrected msProjectPoint() to have a pj_is_latlong()
check on the non-NULL projection object:

      /* nothing to do if the other coordinate system is also lat/long */
      if( in == NULL && out != NULL && pj_is_latlong(out->proj) )
          return MS_SUCCESS;
      if( out == NULL && in != NULL && pj_is_latlong(in->proj) )
          return MS_SUCCESS;

This corrects the error and I now get back proper results from the capabilities
document. 

[fwarmerdam]

As a further step, I have modified msGetWCSCoverageMetadata() to do a
pj_is_latlong() test and not even bother reprojecting if it finds it is
true. 

  // Already in latlong .. use directly.
  if( layer->projection.proj != NULL && pj_is_latlong(layer->projection.proj))
  {
      /* no change */
  }
  else ...

Note, that we should really be creating a real WGS84 coordinate system 
definition and reprojecting to it.  This will ensure that any transformations
of prime meridian or datum will be applied.  The datum shifts aren't too
important since the datum differences are small, and this coverage info
is mainly used for gross location on the earth.  But currently any data
stored in a geographic coordinate system with a non-greenwich prime meridian
will be grossly misplaced. 

I am leaving this open with the intent of revisiting this later. 

[warmerdam]

I shall attempt to deal with this for 5.0 release.

[warmerdam]

The prime meridian case is quite rare, so I'm deferring this to post 5.0.

[warmerdam]

I have attempted this now with the map:

MAP

NAME TEST
STATUS ON
SIZE 800 800
EXTENT -100 -10 -60 40
IMAGECOLOR 255 255 0

IMAGETYPE jpeg

WEB
  METADATA
   # OWS stuff for server
   "ows_updatesequence"   "2007-10-30T14:23:38Z"
   "ows_title"            "First Test Service"
   "ows_fees"             "NONE"
   "ows_accessconstraints" "NONE"
   "ows_abstract"         "Test Abstract"
   "ows_keywordlist"      "keyword,list"
   "ows_service_onlineresource" "http://198.202.74.215/cgi-bin/wcs_demo"
   "ows_contactorganization" "OSGeo"
   "ows_contactperson"    "Frank Warmerdam"
   "ows_contactposition"  "Software Developer"
   "ows_contactvoicetelephone" "(613) 754-2041"
   "ows_contactfacsimiletelephone" "(613) 754-2041x343"
   "ows_address" "3594 Foymount Rd"
   "ows_city" "Eganville"
   "ows_stateorprovince" "Ontario"
   "ows_postcode" "K0J 1T0"
   "ows_country" "Canada"
   "ows_contactelectronicmailaddress" "warmerdam@pobox.com"
   "ows_hoursofservice" "0800h - 1600h EST"
   "ows_contactinstructions" "during hours of service"
   "ows_role" "staff"

   # OGC:WCS
   "wcs_label"    "Test Label"
   "wcs_description" "Test description"
   "wcs_onlineresource"    "http://devgeo.cciw.ca/cgi-bin/mapserv/ecows"
   "wcs_metadatalink_href" "http://devgeo.cciw.ca/index.html"
  END
END

LAYER
  NAME GTOPO
  TYPE raster
  STATUS ON
  DATA /usr2/data/jpeg2000/gtopo.jp2
  DUMP TRUE
#  OFFSITE 255 255 255
  PROJECTION
    "init=epsg:4326"
  END
  METADATA 
    "wcs_title" "GTOPO"
    "wcs_label" "GTOPO - label"
    "wcs_extent" "-100 -10 -60 40"
    "wcs_resolution" "0.0083333 0.0083333"
  END
END

END # of map file

and the request:

../mapserv QUERY_STRING='&map=test.map&SERVICE=WCS&VERSION=1.0.0&REQUEST=GetCapabilities'

and the relavent part of the result:

<ContentMetadata>
  <CoverageOfferingBrief>
    <name>GTOPO</name>
    <label>GTOPO - label</label>
    <lonLatEnvelope srsName="urn:ogc:def:crs:OGC:1.3:CRS84">
      <gml:pos>-100 -10</gml:pos>
      <gml:pos>-60 40</gml:pos>
    </lonLatEnvelope>
  </CoverageOfferingBrief>
</ContentMetadata>

looks fine. I conclude that either the problem has been resolved in the intervening years, or there was some particular required to reproduce the problem that I failed to note.