[GRASS5] PROJ_INFO and prime meridian

GRASS GRASS developer
1

A projection related question:

Currently I try to geocode a (precisely) scanned map from 1945,
Northern Italy. The map was made by the US Army in Lambert Conic
Orthomorphic. The map is based on a map from 1931 done by Istituto
Geografico Militare and was enhanced by aerial photo survey (I assume
also reprojected to Lambert Conic Orthomorphic).

Given on the map is:
- prime meridian Monte Mario (old!): 12d27'7.1E (today 12d27'8.4E)
- "spheroid": Bessel
- Origin: 14E, 45:54N
- False Easing: 800000 meters
- False Northing: 601000 meters
- The map provides the LCC grid and coordinates as well as
  Lat/Long corners in reference to Monte Mario Old and equator.

Using the latest CVS version, I created a location with LCC:
PROJ_INFO
name: Lambert Conformal Conic
towgs84: -104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
proj: lcc
ellps: bessel
a: 6377397.1550000003
es: 0.0066743722
f: 299.1528128000
lat_0: 45.9000000000 <- given 45:54N
lat_1: 33.0000000000 <- calculated by GRASS, given 'return'
lat_2: 45.0000000000 <- calculated by GRASS, given 'return'
lon_0: 14.0000000000 <- given 14E
x_0: 800000.0000000000
y_0: 601000.0000000000

Procedure:
- import of scanned map into XY location
- definition of above LCC location
- i.points of 4 points on the LCC grid, typing into the related
  coordinates by keyboard as read from the map
- i.rectify, 1st order into LCC
- r.proj of recent data from Gauss-Boaga (a tmerc derivate) into
  this LCC location

       i.* r.proj
    XY ======> LCC <======== Gauss-Boaga
   scan final recent GIS data

Validation tests:
- the Lat/Long corner coordinates as printed on the map
  match perfectly with d.where -l results
- but the rectified map is shifted compared to new data (see next)

Problems:
- how to define a different prime meridian (monte Mario old/new)?
- the towgs84 may be wrong (no idea what to use, I selected rome40)
- the resulting map is shifted according to recent data (orthophoto 1m):
    - around 400m to North
    - around 800m to East
- there is also a very small rotation

Ideas (comments welcome):
- Lambert Conic Orthomorphic != LCC (I hope not)
- the optional definition of the true prime meridian may solve the
  rotation problem (the map center is around 11 degree East from Monte
  Mario)
- the bessel definition of 1945 may be different from the current bessel
  in GRASS
- the US Army used a different datum for LCC (where to find out?)
- the Gauss-boaga location with recent GIS data is wrong (I hope not)

Sorry for the long mail,

Markus

2

Hello Markus
Below are a few ideas about the problem

On Wed, 4 Jun 2003, Markus Neteler wrote:

...

Using the latest CVS version, I created a location with LCC:
PROJ_INFO
name: Lambert Conformal Conic
towgs84: -104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
proj: lcc
ellps: bessel
a: 6377397.1550000003
es: 0.0066743722
f: 299.1528128000
lat_0: 45.9000000000 <- given 45:54N

lat_1: 33.0000000000 <- calculated by GRASS, given 'return'
lat_2: 45.0000000000 <- calculated by GRASS, given 'return'

I am actually not sure this is right at all; you could try deleting these
two extra lines. I think lcc should have EITHER lat_0 OR (lat_1 AND lat_2)
but NOT all 3. In fact this looks like a bug in g.setproj as 33 and 45 are
the default values for lat_1 and lat_2 (in src/libes/gis/geo_init.c) and
have not been 'calculated'. Perhaps someone more familiar with the lcc
projection could comment.

lon_0: 14.0000000000 <- given 14E
x_0: 800000.0000000000
y_0: 601000.0000000000

...

Problems:
- how to define a different prime meridian (monte Mario old/new)?

If the longitude of origin is given relative to the Monte Mario prime
meridian, then you would need to add the two together to give the
longitude of origin relative to Greenwich, i.e. 14d0 + 12d27'7.1 =
26d27'7.1"E Is this anywhere near the area covered by the map?

Or maybe the Monte Mario pm simply means the map is referenced to the
rome40 datum, in which case..

- the towgs84 may be wrong (no idea what to use, I selected rome40)

these parameters would be right but the ellipsoid would be wrong, as
rome40 uses international. Is there some other reason it should be bessel?

...

Ideas (comments welcome):
- Lambert Conic Orthomorphic != LCC (I hope not)

I think they are the same

- the Gauss-boaga location with recent GIS data is wrong (I hope not)

If it also is rome40 does it use the same towgs84 parameters?

Paul

3

On Wed, Jun 04, 2003 at 08:16:54PM +0100, Paul Kelly wrote:

Hello Markus
Below are a few ideas about the problem

On Wed, 4 Jun 2003, Markus Neteler wrote:

...
> Using the latest CVS version, I created a location with LCC:
> PROJ_INFO
> name: Lambert Conformal Conic
> towgs84: -104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
> proj: lcc
> ellps: bessel
> a: 6377397.1550000003
> es: 0.0066743722
> f: 299.1528128000
> lat_0: 45.9000000000 <- given 45:54N

> lat_1: 33.0000000000 <- calculated by GRASS, given 'return'
> lat_2: 45.0000000000 <- calculated by GRASS, given 'return'

I am actually not sure this is right at all; you could try deleting these
two extra lines. I think lcc should have EITHER lat_0 OR (lat_1 AND lat_2)
but NOT all 3. In fact this looks like a bug in g.setproj as 33 and 45 are
the default values for lat_1 and lat_2 (in src/libes/gis/geo_init.c) and
have not been 'calculated'. Perhaps someone more familiar with the lcc
projection could comment.

In fact it is not obvious (yet) how to define LCC 1SP or LCC 2SP.
So I left it empty.

New test:
"old" PROJ_INFO:

g.region -l
long: 11.31206 lat: 46.51243 (north/west corner)
long: 11.51722 lat: 46.51646 (north/east corner)
long: 11.52172 lat: 46.40260 (south/east corner)
long: 11.31693 lat: 46.39858 (south/west corner)
rows: 127
cols: 158
Center Longitude: 11:25:01.148843E [11.41699]
Center latitude: 46:27:27.066733N [46.45752]

Deleting
lat_1: 33.0000000000
lat_2: 45.0000000000

g.region -l
Unable to initialise PROJ.4 with the following parameter list:
proj=lcc lat_0=45.9000000000 lon_0=14.0000000000 x_0=800000.0000000000
y_0=601000.0000000000 a=6377397.155000 rf=299.152813
towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
The error message was 'conic lat_1 = -lat_2'
ERROR: Can't get projection key values of current location

http://www.remotesensing.org/geotiff/proj_list/
-> LCC 1SP
+proj=lcc +lat_1=Latitude of natural origin
             +lon_0=Longitude of natural origin
             +k_0=Scale factor at natural origin
             +x_0=False Origin Easting
             +y_0=False Origin Northing

So this should work, but... (I didn't try yet with 'proj' nor cs2cs).

> lon_0: 14.0000000000 <- given 14E
> x_0: 800000.0000000000
> y_0: 601000.0000000000
>
...
> Problems:
> - how to define a different prime meridian (monte Mario old/new)?

If the longitude of origin is given relative to the Monte Mario prime
meridian, then you would need to add the two together to give the
longitude of origin relative to Greenwich, i.e. 14d0 + 12d27'7.1 =
26d27'7.1"E Is this anywhere near the area covered by the map?

No, the area is at 11dsomething' E (see above).

According to:
http://www.remotesensing.org/proj/gen_parms.html
"pm - Prime Meridian

A prime meridian may be declared indicating the offset between the prime
meridian of the declared coordinate system and that of greenwich. A prime
meridian is clared using the "pm" parameter, and may be assigned a symbolic
name, or the longitude of the alternative prime meridian relative to
greenwich.

Currently prime meridian declarations are only utilized by the
pj_transform() API call, not the pj_inv() and pj_fwd() calls. Consequently
the user utility cs2cs does honour prime meridians but the proj user utility
ignores them. "

Monte Mario prime meridian is hardcoded in PROJ4:

grep rome *.c
pj_datums.c: "rome", "12d27'8.4\"E",

(where I would need 12d27'7.1). Maybe proj4 needs a custom prime
meridian (couldn't find if it is already there).

Or maybe the Monte Mario pm simply means the map is referenced to the
rome40 datum, in which case..

> - the towgs84 may be wrong (no idea what to use, I selected rome40)

these parameters would be right but the ellipsoid would be wrong, as
rome40 uses international. Is there some other reason it should be bessel?

Right, rome40 uses international, but the US Army defined Bessel.
So it might be a wrong datum. Probably I have to ask on the PROJ4
list as well (tomorrow).

...
> Ideas (comments welcome):
> - Lambert Conic Orthomorphic != LCC (I hope not)

I think they are the same

> - the Gauss-boaga location with recent GIS data is wrong (I hope not)

If it also is rome40 does it use the same towgs84 parameters?

In fact the reference location in Gauss-Boaga uses rome40 datum
with international ellipsoid. So, yes, the same towgs84 which could
indicate a problem.

Thanks for your comments,

Markus

4

Markus Neteler wrote:
> A projection related question:

> - prime meridian Monte Mario (old!): 12d27'7.1E (today 12d27'8.4E)
> - "spheroid": Bessel
> - Origin: 14E, 45:54N

Origin Meridian is not the prime Meridian... It's unusual

> - False Easing: 800000 meters
> - False Northing: 601000 meters

Strange False northing too (did the surveyors find a nice cave full
of good wine ? ;-))

[...]
> lat_0: 45.9000000000 <- given 45:54N
> lat_1: 33.0000000000 <- calculated by GRASS, given 'return'
> lat_2: 45.0000000000 <- calculated by GRASS, given 'return'

Hmm... The problem may be here : I don't understand the values of lat_1
and lat_2 here : The conic projection may be tangent => Lat_0, LAt_1 and Lat_2
should be the same, or the conic projection is secant => Lat_0 is between
Lat_1 qnd Lat_2... BUT some peoples use a tangent projection, but with
a scale (at the tangent parallel) lesser than 1 (0.998...) which gives
results very closes of a secant projection (that's true for the french
LAmbert projection

> Ideas (comments welcome):
> - Lambert Conic Orthomorphic != LCC (I hope not)
I think it's synonyms

> - the optional definition of the true prime meridian may solve the
> rotation problem (the map center is around 11 degree East from Monte
> Mario)

uses 12°27'7.1" in r.proj ? (but how use 14° for false easting ???...
In the same way look at the values for Lat_1 and Lat_2 (use 45.9 ?)

> - the bessel definition of 1945 may be different from the current bessel
> in GRASS

seems unprobable, the name would be different

> - the US Army used a different datum for LCC (where to find out?)

But I doubt they do else than put new détails on an old map, thus
respecting the original projection...

> - the Gauss-boaga location with recent GIS data is wrong (I hope not)

Maybe, but difficult to verify, execpt if you have the opportunity to
go with a gps receiver on a easely identifiable point on both old
and new data and look what the GPS will give (a simple Garmin or Magellan
for outdoor will be enough to give you the name of the faulty set, since
the precision of this kind of GPS is about 10 m, comparing to your 400/800m)

Hoping this will help you

Good luck,

--
Michel Wurtz - Auzeville-Tolosane

5

(cc grass5)

On Wed, Jun 04, 2003 at 08:54:45PM +0100, Paul Kelly wrote:

On Wed, 4 Jun 2003, Markus Neteler wrote:

> > > lat_0: 45.9000000000 <- given 45:54N
> >
> > > lat_1: 33.0000000000 <- calculated by GRASS, given 'return'
> > > lat_2: 45.0000000000 <- calculated by GRASS, given 'return'
> >
...
> Deleting
> lat_1: 33.0000000000
> lat_2: 45.0000000000
>
> g.region -l
> Unable to initialise PROJ.4 with the following parameter list:
> proj=lcc lat_0=45.9000000000 lon_0=14.0000000000 x_0=800000.0000000000
> y_0=601000.0000000000 a=6377397.155000 rf=299.152813
> towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
> The error message was 'conic lat_1 = -lat_2'
> ERROR: Can't get projection key values of current location
>

Mmmm obvious suggestion then:
delete lat_0 and lat_2 from PROJ_INFO; add
lat_1: 45.9

Is this any better?

Yes. Now the shift is reduces to 70m eastwards and 170m northwards.

The reference Orthophoto I have verified with a vector street map from
a different data provider, they match.

I have made a screenshot of a zoomed area:
http://mpa.itc.it/markus/tmp/lcc_1945_ortho2000.png
(1MB)
http://mpa.itc.it/markus/tmp/lcc_1945_ortho2000_small.png
(500kb)

So maybe now only the rome40 datum for the LCC location
is wrong (and the prime meridian slightly different).

Also the first time you tried it (the original PROJ_INFO) what did r.proj
print out as the parameters being used? Did it have lat_1 or lat_0 and
lat_2?

Reverting the PROJ_INFO to the original version, it prints:
r.proj in=268_09ne.rect out=test loc=bzLCC map=PERMANENT

Input Projection Parameters: +proj=lcc +lat_0=45.9000000000
+lat_1=33.0000000000 +lat_2=45.0000000000 +lon_0=14.0000000000
+x_0=800000.0000000000 +y_0=601000.0000000000 +a=6377397.155000
+rf=299.152813 +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

Output Projection Parameters: +proj=tmerc +lat_0=0.0000000000 +lon_0=9.0
+k_0=0.9996000000 +x_0=500000.0 +y_0=-5000000.0 +a=6378388.000000
+rf=297.000000 +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

Markus