Opened 20 years ago
Closed 14 years ago
#647 closed defect (worksforme)
Reprojection for Lat/Long Wrong
Reported by: | warmerdam | Owned by: | warmerdam |
---|---|---|---|
Priority: | normal | Milestone: | FUTURE |
Component: | WCS Server | Version: | unspecified |
Severity: | normal | Keywords: | |
Cc: |
Description (last modified by )
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 (5)
comment:2 by , 20 years ago
Status: | new → assigned |
---|
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.
comment:3 by , 17 years ago
Description: | modified (diff) |
---|---|
Milestone: | → 5.0 release |
Priority: | high → normal |
I shall attempt to deal with this for 5.0 release.
comment:4 by , 17 years ago
Milestone: | 5.0 release → FUTURE |
---|
The prime meridian case is quite rare, so I'm deferring this to post 5.0.
comment:5 by , 14 years ago
Resolution: | → worksforme |
---|---|
Status: | assigned → closed |
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.
Note:
See TracTickets
for help on using tickets.