[Proj] problem in gdalwrap wrt prime meridians

Sun Jun 24 13:47:04 EDT 2007

Hi,
I'm trying to reproject a geotiff from "lambert 2 etendu" to "lambert 93" (French projections) using gdalwarp from FWTools 1.3.2 under windows XP SP2.
Here is the info concerning the input file.

gdalinfo test.tif

Driver: GTiff/GeoTIFF
Files: test.tif
Size is 2000, 2000
Coordinate System is:
PROJCS["NTF (Paris) / Lambert zone II",
    GEOGCS["NTF (Paris)",
        DATUM["Nouvelle_Triangulation_Francaise_Paris",
            SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936265,
                AUTHORITY["EPSG","7011"]],
            AUTHORITY["EPSG","6807"]],
        PRIMEM["Paris",2.33722917],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4807"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",46.8],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",0.99987742],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",2200000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27572"]]
Origin = (827000.000000000000000,2224000.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DOCUMENTNAME=39-2006-0820-2221-LA2E07
  TIFFTAG_XRESOLUTION=50.799988
  TIFFTAG_YRESOLUTION=50.799988
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left  (  827000.000, 2224000.000) (  2d59'3.48"E, 46d58'37.71"N)
Lower Left  (  827000.000, 2223000.000) (  2d59'1.68"E, 46d58'5.34"N)
Upper Right (  828000.000, 2224000.000) (  2d59'50.76"E, 46d58'36.47"N)
Lower Right (  828000.000, 2223000.000) (  2d59'48.96"E, 46d58'4.11"N)
Center      (  827500.000, 2223500.000) (  2d59'26.22"E, 46d58'20.91"N)
Band 1 Block=2000x4 Type=Byte, ColorInterp=Red
Band 2 Block=2000x4 Type=Byte, ColorInterp=Green
Band 3 Block=2000x4 Type=Byte, ColorInterp=Blue

I try to reproject this image using the following command, giving the input projection system with a string being the exact copy of the EPSG definition for code 27572:

D:\>gdalwarp -s_srs "+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-
60,320,0,0,0,0 +pm=paris +units=m +no_defs" -t_srs "+init=epsg:2154" -tps -rc -tr
 0.5 0.5 -co INTERLEAVE=PIXEL test.tif test_res.tif

D:\>gdalinfo test_res.tif
Driver: GTiff/GeoTIFF
Files: test_res.tif
Size is 2040, 2040
Coordinate System is:
PROJCS["RGF93 / Lambert-93",
    GEOGCS["RGF93",
        DATUM["Reseau_Geodesique_Francais_1993",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6171"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4171"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",49],
    PARAMETER["standard_parallel_2",44],
    PARAMETER["latitude_of_origin",46.5],
    PARAMETER["central_meridian",3],
    PARAMETER["false_easting",700000],
    PARAMETER["false_northing",6600000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2154"]]
Origin = (1053910.241789964700000,6663452.274245270500000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( 1053910.242, 6663452.274) (  7d39'29.83"E, 46d58'38.29"N)
Lower Left  ( 1053910.242, 6662432.274) (  7d39'26.99"E, 46d58'5.29"N)
Upper Right ( 1054930.242, 6663452.274) (  7d40'18.05"E, 46d58'36.34"N)
Lower Right ( 1054930.242, 6662432.274) (  7d40'15.20"E, 46d58'3.34"N)
Center      ( 1054420.242, 6662942.274) (  7d39'52.52"E, 46d58'20.82"N)
Band 1 Block=2040x1 Type=Byte, ColorInterp=Red
Band 2 Block=2040x1 Type=Byte, ColorInterp=Green
Band 3 Block=2040x1 Type=Byte, ColorInterp=Blue
Creating output file that is 2040P x 2040L.
Processing input file test.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.

As one can see, the upper left corner is shifted ~ 2,3° east. The result should be near 5d19' and not 7d39'. It looks as if the prime meridian has been accounted for twice. This is incoherent with cs2cs, which gives the right result:

D:\>cs2cs -v +init=epsg:27572 +to +init=epsg:2154
# ---- From Coordinate System ----
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +init=epsg:27572 +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
# +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515
# +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs
# ---- To Coordinate System ----
#Lambert Conformal Conic
#       Conic, Sph&Ell
#       lat_1= and lat_2= or lat_0
# +init=epsg:2154 +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3
# +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
# +no_defs
827000 2224000
876410.07       6655580.63 41.71

In addition, it seems that gdalwarp "knows" an internal definition of EPSG codes giving the right results:

D:\>gdalwarp -s_srs EPSG:27572 -t_srs "+init=epsg:2154" -tps -rc
 -tr 0.5 0.5 -co INTERLEAVE=PIXEL test.tif test_res2.tif
Creating output file that is 2015P x 2015L.
Processing input file test.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.

D:\dalles_lambert93>gdalinfo test_res2.tif
Driver: GTiff/GeoTIFF
Files: test_res2.tif
Size is 2015, 2015
Coordinate System is:
PROJCS["RGF93 / Lambert-93",
    GEOGCS["RGF93",
        DATUM["Reseau_Geodesique_Francais_1993",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6171"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4171"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",49],
    PARAMETER["standard_parallel_2",44],
    PARAMETER["latitude_of_origin",46.5],
    PARAMETER["central_meridian",3],
    PARAMETER["false_easting",700000],
    PARAMETER["false_northing",6600000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2154"]]
Origin = (876401.507783832840000,6655580.627617656300000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  876401.508, 6655580.628) (  5d19'15.01"E, 46d58'37.58"N)
Lower Left  (  876401.508, 6654573.128) (  5d19'13.61"E, 46d58'4.94"N)
Upper Right (  877409.008, 6655580.628) (  5d20'2.70"E, 46d58'36.62"N)
Lower Right (  877409.008, 6654573.128) (  5d20'1.29"E, 46d58'3.98"N)
Center      (  876905.258, 6655076.878) (  5d19'38.15"E, 46d58'20.78"N)
Band 1 Block=2015x1 Type=Byte, ColorInterp=Red
Band 2 Block=2015x1 Type=Byte, ColorInterp=Green
Band 3 Block=2015x1 Type=Byte, ColorInterp=Blue

However, this is not a solution, as I would like to use my own definitions for projections (I want to use nadgrids for French projections)...

There's definitely a big problem here, don't you think so ?

Best regards,
Fabrice.

      _____________________________________________________________________________ 
Ne gardez plus qu'une seule adresse mail ! Copiez vos mails vers Yahoo! Mail 