headers
pcreso | 15 Apr 2012 22:19

EPSG:3994 problem

Hi,

I'm unable to get Proj to correctly return lat/long coordinates from Mercator41 data.
http://www.spatialreference.org/ref/epsg/3994/ does not list any proj4 string, Postgis proj4text or Mapserver values for this projection.

I have an ESRI binary grid with a projection definition (prj.adf) of:

Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters   
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

gdalinfo fails to interpret this correctly, giving:

gdalinfo w001001.adf
Driver: AIG/Arc/Info Binary Grid
Files: .
       ./#m41_to_ll#
       ./sta.adf
       ./w001001x.adf
       ./hdr.adf
       ./w001001.adf
       ./m41_to_ll
       ./whbr
       ./prj.adf
       ./dblbnd.adf
Size is 8817, 7901
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",100],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",-41],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]
Origin = (6199987.499999998137355,-3709987.500002769753337)
Pixel Size = (25.000000000000000,-25.000000000000000)
Corner Coordinates:
Upper Left  ( 6199987.500,-3709987.500)
Lower Left  ( 6199987.500,-3907512.500)
Upper Right ( 6420412.500,-3709987.500)
Lower Right ( 6420412.500,-3907512.500)
Center      ( 6310200.000,-3808750.000)
Band 1 Block=256x4 Type=Float32, ColorInterp=Undefined
  Min=-2570.453 Max=2917.197
  NoData Value=-3.4028234663852886e+38

With some obvious errors, like latitude of origin = 100... so just extracting the xyz values & trying to reproject with proj, I'm having problems.

I tried: g dal2xyz.py w001001.adf | head
which returns:
6200000.000 -3710000.000 76.1574
6200025.000 -3710000.000 76.1842
...

which seems to work fine, but when I pass these Mercator41 coordinates into proj I'm unable to get a parameter list which gives the lat/long values for where I know these data sit. There is an offset of many miles. This suggests either my source data has an error in the coordinates (unlikely) or I don't understand enough about proj parameters (more likely)

Based on the contents of the above prj.adf file I have tried the following command:

echo "6200000.000 -3710000.000" | proj -If "%.5f" +proj=merc +lon_0=100 +lat_ts=-41 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
173.69093       -40.48353

I expected values around 174.66  -41.41 (the data is a 25m multibeam bathymetry grid around Cook Strait, between the North & South Islands of New Zealand)


Any advice appreciated.

  Brent Wood



 
_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply |
headers
Mikael Rittri | 16 Apr 2012 11:27
Favicon

Re: EPSG:3994 problem

Hello Brent,

 

I don’t know why gdalinfo misinterprets the ESRI binary grid, or why spatialreference.org doesn’t help.

But EPSG:3994 uses Mercator variant B (see http://www.epsg.org/guides/docs/G7-2.pdf), which is a

somewhat unusual variant.

 

On the other hand, the epsg file of Proj version 4.8.0 has the same definition that you tried:

 

# WGS 84 / Mercator 41

<3994> +proj=merc +lon_0=100 +lat_ts=-41 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs  <>

 

which I believe is correct (lat_ts means latitude of true scale). Based on that, I can compute:

 

>proj +init=epsg:3994

 

100 0

0.00    0.00

 

101 41

84135.19        3767131.99

 

101 41.000009

84135.19        3767132.99

 

The first test point shows that the grid origin (zero easting, zero northing) is at

the equator at longitude 100 d E, which is what I would expect.

 

The other two test points shows that a difference of 0.000009 degrees in latitude,

near the latitude 41 N, gives a difference of about 1 meter in northing. This is also

what I would expect, since 41 N (and 41 S) are standard parallels, where the scale

should be true.  (Because 0.000009 degrees of latitude is about 1 meter.)

 

For your test point, I reproduce your results:

 

>proj -I +init=epsg:3994 -f "%.5f"

 

6200000.000 -3710000.000

173.69093       -40.48353

 

You say you didn’t expect these results:

 

> I expected values around 174.66  -41.41

 

but I don’t know if what you expected is correct or not. I have used Proj 4.8.0 that came with OSGeo4W.

 

On the other hand, in Proj 4.7.0 the epsg file had the entry

 

# WGS 84 / Mercator 41

<3994> +proj=merc +lon_0=100 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs  <>

 

which I think is very wrong.

    The EPSG definition of WGS 84 / Mercator 41 used to have the code 3752, but

this was deprecated in 2008 and replaced with 3994 (EPSG change request 2008.088).

I am not sure if the bad definition in Proj 4.7.0 had anything to do with that, though.

 

Best regards,

 

Mikael Rittri

Carmenta

Sweden

http://www.carmenta.com

 

From: proj-bounces <at> lists.maptools.org [mailto:proj-bounces <at> lists.maptools.org] On Behalf Of pcreso <at> pcreso.com
Sent: den 15 april 2012 22:19
To: proj <at> lists.maptools.org
Subject: [Proj] EPSG:3994 problem

 

Hi,

I'm unable to get Proj to correctly return lat/long coordinates from Mercator41 data.
http://www.spatialreference.org/ref/epsg/3994/ does not list any proj4 string, Postgis proj4text or Mapserver values for this projection.

I have an ESRI binary grid with a projection definition (prj.adf) of:

Projection    MERCATOR
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0
Parameters   
 100  0  0.0 /* longitude of central meridian
 -41  0  0.0 /* latitude of true scale
0.0 /* false easting (meters)
0.0 /* false northing (meters)

gdalinfo fails to interpret this correctly, giving:

gdalinfo w001001.adf
Driver: AIG/Arc/Info Binary Grid
Files: .
       ./#m41_to_ll#
       ./sta.adf
       ./w001001x.adf
       ./hdr.adf
       ./w001001.adf
       ./m41_to_ll
       ./whbr
       ./prj.adf
       ./dblbnd.adf
Size is 8817, 7901
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",100],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",-41],
    PARAMETER["false_northing",0],
    UNIT["METERS",1]]
Origin = (6199987.499999998137355,-3709987.500002769753337)
Pixel Size = (25.000000000000000,-25.000000000000000)
Corner Coordinates:
Upper Left  ( 6199987.500,-3709987.500)
Lower Left  ( 6199987.500,-3907512.500)
Upper Right ( 6420412.500,-3709987.500)
Lower Right ( 6420412.500,-3907512.500)
Center      ( 6310200.000,-3808750.000)
Band 1 Block=256x4 Type=Float32, ColorInterp=Undefined
  Min=-2570.453 Max=2917.197
  NoData Value=-3.4028234663852886e+38

With some obvious errors, like latitude of origin = 100... so just extracting the xyz values & trying to reproject with proj, I'm having problems.

I tried: g dal2xyz.py w001001.adf | head
which returns:
6200000.000 -3710000.000 76.1574
6200025.000 -3710000.000 76.1842
...

which seems to work fine, but when I pass these Mercator41 coordinates into proj I'm unable to get a parameter list which gives the lat/long values for where I know these data sit. There is an offset of many miles. This suggests either my source data has an error in the coordinates (unlikely) or I don't understand enough about proj parameters (more likely)

Based on the contents of the above prj.adf file I have tried the following command:

echo "6200000.000 -3710000.000" | proj -If "%.5f" +proj=merc +lon_0=100 +lat_ts=-41 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
173.69093       -40.48353

I expected values around 174.66  -41.41 (the data is a 25m multibeam bathymetry grid around Cook Strait, between the North & South Islands of New Zealand)


Any advice appreciated.

  Brent Wood



 

 

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply |

Gmane