Opened 21 years ago

Closed 16 years ago

#367 closed defect (fixed)

gdalwarp problem: output is shifted by 900 kilometers (wrong handling of +pm value)

Reported by: Markus Neteler Owned by: warmerdam
Priority: high Milestone:
Component: GDAL_Raster Version: unspecified
Severity: normal Keywords:
Cc: paul-grass@…, Mateusz Łoskot, martinoty

Description (last modified by Mateusz Łoskot)

At time I struggle to reproject a LANDSAT TM 5 subscene with
gdalwarp (from UTM to Gauss-Boaga). The resulting map is
shifted around 920km to West and also slightly distorted.

The input map 'lsat1_utm.tif' was created with gdal_translate from
the original LANDSAT TM 5 scene (from Geotiff to Geotiff). The
smaller subscene is located in Nothern Italy and looks like this:

gdalinfo lsat1_utm.tif
Driver: GTiff/GeoTIFF
Size is 3816, 3461
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (636861.000000,5152686.000000)
Pixel Size = (28.500000,-28.500000)
Corner Coordinates:
Upper Left  (  636861.000, 5152686.000) ( 10d47'2.74"E, 46d30'49.78"N)
Lower Left  (  636861.000, 5054047.500) ( 10d45'20.60"E, 45d37'35.33"N)
Upper Right (  745617.000, 5152686.000) ( 12d12'0.22"E, 46d28'58.64"N)
Lower Right (  745617.000, 5054047.500) ( 12d 8'57.27"E, 45d35'47.57"N)
Center      (  691239.000, 5103366.750) ( 11d28'20.39"E, 46d 3'25.67"N)
Band 1 Block=3816x2 Type=Byte, ColorInterp=Gray


The Lat/Long corners look o.k. But then, after gdalwarp:

#26591: Monte Mario (Rome) / Italy zone 1
gdalwarp -t_srs '+init=epsg:26591
+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68' \
  lsat1_utm.tif lsat1_gb.tif

The datum comes from the official site:
http://crs.ifag.de/

gdalinfo lsat1_gb2.tif
Driver: GTiff/GeoTIFF
Size is 4311, 4015
Coordinate System is:
PROJCS["Monte Mario (Rome) / Italy zone 1",
    GEOGCS["Monte Mario (Rome)",
        DATUM["Monte_Mario_Rome",
            SPHEROID["International 1924",6378388,297.000000000005,
                AUTHORITY["EPSG","7022"]],
            AUTHORITY["EPSG","6806"]],
        PRIMEM["Rome",12.45233333333333],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4806"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3.45233333333333],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26591"]]
Origin = (666266.468468,5206815.783256)
Pixel Size = (28.697167,-28.697167)
Corner Coordinates:
Upper Left  (  666266.468, 5206815.783) ( 14d19'12.01"W, 46d29'41.38"N)
Lower Left  (  666266.468, 5091596.659) ( 14d 7'12.46"W, 45d28'34.07"N)
Upper Right (  789979.954, 5206815.783) ( 12d43'50.79"W, 46d38'10.47"N)
Lower Right (  789979.954, 5091596.659) ( 12d33'32.53"W, 45d36'45.44"N)
Center      (  728123.211, 5149206.221) ( 13d25'57.64"W, 46d 3'27.81"N)
Band 1 Block=4311x1 Type=Byte, ColorInterp=Gray

Looking at LatLong here, it is heavily misplaced (especially shifted to
West).

To me "Central Meridian" in the Gauss-Boaga file looks wrong.

Looking into
/usr/local/share/proj/epsg
I am not sure if
+lon_0=21.45233333333333
is right. Shouldn't it be
+lon_0=9.0
?

The LANDSAT UTM subscene is in
gdal.velocet.ca
pub/incoming/lsat1_utm.tif.bz2

Any hint is welcome to reproject this LANDSAT subscene.

Best regards
 Markus Neteler

Change History (11)

comment:1 by warmerdam, 21 years ago

Markus, 

I appologise for losing track of this.  The problem relates to the use of the
rome prime meridian, but I got distracted and never came back to dig all the
way through this.  I will try to look into it today.


comment:2 by neteler@…, 21 years ago

Thanks, Frank.

Related to this: it were great if there is a possibility to
enter also a "user defined prime meridian" for historic maps etc.
Maybe I should post in in the PROJ4 bugzilla.

Markus

comment:3 by warmerdam, 21 years ago

Markus, 

You can use a user defined prime meridian by just giving a longitude to the
+pm keyword.  Eg. +pm=-17 to use a prime meridian at 17 degrees west of
greenwich.

I looked through this and found a number of problems with prime meridian handling.
According to the document I posted at the following URL, the longitude values
for projection parameters in WKT should be relative to Greenwich, not the
prime meridian of the GEOGCS.  

  http://bugzilla.remotesensing.org/process_bug.cgi

However, OGRSpatialReference class seems to make the opposite assumption, and
methods like SetTM() seem to be trying to convert longitudes from a greenwich
relative value to a prime meridian relative value.  That is to say "normalized"
form according the SetNormProjParm() is relative to greenwich, but the internal
storage form in WKT is relative to the GEOGCS prime meridian. 

The code to initialize from EPSG gives a longitude that is relative to
greenwich, but that appears to be accomplished by taking EPSG longitudes that
are relative to greenwich, offseting them relative to the prime meridian in the
ogr_from_epsg.cpp code, and then having the SetNormProjParm() method undo that
shift again. 

It seems to me that the longitude values in PROJ.4 format need to be relative 
to the prime meridian.  But the OGRSpatialReference::importFromProj4() method
does uses SetNormProjParm() to place PROJ.4 longitude values into WKT.  This
would appear to apply a translation in exactly the *wrong* direction.  I think
that OSR_GDV() in ogr_srs_proj4.cpp should be offseting longitude values to
make them relative to greenwich before passing them to methods like SetTM(). 

The OGRSpatialReference::exportToProj4() method also does just the wrong thing,
using GetNormProjParm() to fetch longitude values from the OGRSpatialReference.
Thus the values get set in PROJ.4 format with longitudes relative to greenwich
instead of relative to the prime meridian as would be required in PROJ.4.

What I need to do is:
 1) Verify whether WKT projection PARAMETER longitudes are supposed to be
    relative to the prime meridian or greenwich once and for all. 

 2) Correct the SetTM() and other projection specific set methods to take the
    arguments given (which should be relative to greenwich) and transform them 
    (or not) as appopriate according to the first item. 

 3) correct EPSG reading code to *not* apply any offset from greenwich to 
    longitudes before setting them with stuff like SetTM().  The EPSG longitudes
    are relative to greenwich, not the prime meridian in effect (ok ... please
    verify). 

 4) Ensure that the PROJ.4 translation (to and from) applies greenwich offsets
    after getting, or before setting since PROJ.4 longitudes need to be relative
    to their prime meridian (please verify). 

 5) Review to ensure we haven't hosted prime meridian handling for FME, GeoTIFF,
    MapInfo or any other general purpose SRS translators. 


comment:4 by neteler@…, 21 years ago

Frank,

[the PROJ related URL is not valid - your previous message]

I tried the prime meridian switch with no avail:

Projection from lat/long to GaussBoaga (Rome40/GB/Zone1):
a) current prime meridian Monte Mario

cs2cs +proj=latlong +datum=WGS84 +to +init=epsg:26591
+towgs84=-85.88,-28.85,+49.67,-1.003,-2.383,-1.808,-27.82 +pm="12d27'8.4E"
11d10'01E 46d02'54N
1667672.33      5101723.06 -0.21


b) old prime meridian Monte Mario
cs2cs +proj=latlong +datum=WGS84 +to +init=epsg:26591
+towgs84=-85.88,-28.85,+49.67,-1.003,-2.383,-1.808,-27.82 +pm="12d27'7.1E"
11d10'01E 46d02'54N
1667672.33      5101723.06 -0.21

I feel that "+pm=..." is ignored (CVS version of PROJ4).

Markus

comment:5 by paul-grass@…, 17 years ago

(In reply to comment #3)
Hello Frank,
I wonder if you ever got anywhere near the bottom of this. I was just about to report a new bug, but had a quick browse through bugzilla and you same to have come across the issue before. Included below is the bug I was going to report.

Paul

Inconsistent handling of non-Greenwich Prime Meridians between OSRImportFromProj4() and OSRExportToProj4().

OSRImportFromProj4() puts longitude of origin relative to Greenwich into the WKT string, rather than preserving longitude relative to the Prime Meridian specified by +pm=xxxx. This would actually be OK if OSRExportToProj4() also works the same way, but it doesn't. It puts the longitude unchanged into the PROJ.4 string and also adds the +pm value, so PROJ.4 will treat the longitude as relevant to the PM.

Either way would be acceptable, but it's the fact that the two work inconsistently that's the problem - converting from one format (WKT or PROJ.4) to the other and back again shifts the value of longitude of origin by the difference between the PM and the Greenwich PM!

I'm not sure what's the best way to solve this. ISTR some discussion in the past that the longitiude values in a WKT file were *always* relative to Greenwich, regardless of the value of the PRIMEM key, but in that case if OSRExportToProj4() is going to include the +pm in the PROJ string it should also adjust the longitude.

I *think* my preferred solution would be that the when converting to PROJ.4 format, the prime meridian information is simply stripped out and the longitude values adjusted accordingly (if necessary, depending on how it's "supposed" to be specified in the WKT!). I mean who really wants longitudes referenced to a non-Greenwich meridian anyway? ;) IMHO it just leads to endless confusion!

Test C program
==============

--START------------------------------------------------------------
#include <stdio.h>
#include <ogr_srs_api.h>

const char *original = "+proj=krovak +a=6377397.155 +rf=299.1528128 "
  "+lat_0=49.5 +lon_0=42.5 +alpha=30.288129722222 +k=0.9999 +x_0=0 +y_0=0 "
  "+pm=ferro +no_defs";

int main()
{
    char *translated, *wkt;
    OGRSpatialReferenceH hSRS;

    fprintf(stderr, "Original PROJ.4 string:\n%s\n", original);

    hSRS = OSRNewSpatialReference(NULL);
    OSRImportFromProj4( hSRS, original );

    OSRExportToPrettyWkt( hSRS, &wkt, 0 );

    fprintf(stderr, "Intermediate WKT (note longitude of origin already "
                    "changed):\n%s\n", wkt);

    OSRExportToProj4( hSRS, &translated );

    fprintf(stderr, "Translated PROJ.4 string:\n%s\n", translated);

    return 0;
}
--END--------------------------------------------------------------

Output
======

Original PROJ.4 string:
+proj=krovak +a=6377397.155 +rf=299.1528128 +lat_0=49.5 +lon_0=42.5 +alpha=30.28
8129722222 +k=0.9999 +x_0=0 +y_0=0 +pm=ferro +no_defs
Intermediate WKT (note longitude of origin already changed):
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6377397.155,299.1528128]],
        PRIMEM["ferro",-17.666666666668],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.833333333332],
    PARAMETER["azimuth",30.288129722222],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Translated PROJ.4 string:
+proj=krovak +lat_0=49.5 +lon_0=24.833333333332 +alpha=30.288129722222 +k=0.9999
 +x_0=0 +y_0=0 +ellps=bessel +pm=-17.666666666668 +units=m +no_defs

The +lon_0 values are different!

comment:6 by martinoty, 17 years ago

Summary: gdalwarp problem: output is shifted by 900 kilometersgdalwarp problem: output is shifted by 900 kilometers (wrong handling of +pm value)

I agree with your interpretation of the bug, however I think the problem is more related to ImportFromProj4() than to ExportToProj4(). Indeed, in ImportFromProj4(), dfFromGreenwich is added to lon_0, and at the same time the prime meridian value is not modified: this is inconsistent. One should either not add the prime meridian value and keep the +pm information, or define the prime meridian value as "Greenwhich" once the value is added to lon_0.

A workaround could be: In ogr/ogr_srs_proj4.cpp#L731 : After :

/* -------------------------------------------------------------------- */

/* Set the ellipsoid information. */

/* -------------------------------------------------------------------- */

Add:

dfFromGreenwich = 0.0;

pszPM = "Greenwich";

I think the severity of this bug is rather high, since every projection using +pm is handled wrongly.

comment:7 by Mateusz Łoskot, 16 years ago

Cc: Mateusz Łoskot added
Description: modified (diff)

comment:8 by Markus Neteler, 16 years ago

Reporter: changed from neteler@… to Markus Neteler

comment:9 by Even Rouault, 16 years ago

Cc: martinoty added

I think this bug (the part related to lon_0 and PM handling) is fixed since GDAL 1.4.4, at r12535 (for bug #1940)

comment:10 by martinoty, 16 years ago

You're right, this bug is indeed fixed since GDAL 1.4.4. I'm positive you can close this bug.

comment:11 by Even Rouault, 16 years ago

Resolution: fixed
Status: assignedclosed
Note: See TracTickets for help on using tickets.