Opened 20 years ago

Last modified 14 years ago

#647 closed defect

Reprojection for Lat/Long Wrong — at Initial Version

Reported by: warmerdam Owned by: sdlime
Priority: normal Milestone: FUTURE
Component: WCS Server Version: unspecified
Severity: normal Keywords:
Cc:

Description

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;
  }

Change History (0)

Note: See TracTickets for help on using tickets.