[GRASS-dev] [bug #5454] (grass) ps.map scaling not correct when projection units are feet

Justin_Hickey_via_RT · January 28, 2007, 11:51pm

this bug's URL: http://intevation.de/rt/webrt?serial_num=5454
-------------------------------------------------------------------------

Subject: ps.map scaling not correct when projection units are feet

Platform: GNU/Linux/x86
grass obtained from: CVS
grass binary for platform: Compiled from Sources
GRASS Version: CVS HEAD 20070127 at 04:18 UTC

Noticed a strange anomaly in ps.map output when working in a State Plane Coordinate System with coordinates in feet. Specifying any scale in ps.map gives manifestly incorrect scaling. The scalebar produced also bears no obvious relationship to the scale of the map.

To reproduce: create a location in New Mexico State Plane Coordinates, Central Zone, NAD83 (EPSG number 2258). The result will have a PROJ_INFO file like this:
name: Transverse Mercator
proj: tmerc
datum: nad83
towgs84: 0.000,0.000,0.000
ellps: grs80
lat_0: 31
lon_0: -106.25
k: 0.999900
x_0: 500000.0001016001
y_0: 0
no_defs: defined

and a PROJ_UNITS file of:
unit: US survey foot
units: US survey foots
meters: 0.3048006096012192

Set the region to something arbitrary, for example:

g.region -p

projection: 99 (Transverse Mercator)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1485814.68884
south: 1485263.50884
west: 1525396.34529
east: 1526171.30805
nsres: 0.61310345
ewres: 0.61310345
rows: 899
cols: 1264
cells: 1136336

Now use the following ps.map file to generate postscript:

paper a0
end
scalebar f
where 1.5 1.5
length 500
height .25
segment 5
numbers 2
fontsize 12
end
mapinfo
  where 11.5 15
  fontsize 12
end
grid 100
  color black
  numbers 2
end

This should display a 100 foot grid and nothing else, and should have a
scale bar with 100 foot markings that match the grid size.

For the region specified above, ps.map will select "1:982" as the scale. Now try tinkering with the scale by adding scale commands to the ps.map file. Even selecting "1:982" will give completely different results than leaving the scale out altogether. As you change scale, the scalebar will not change size to reflect the change in the map grid.

I dug around a little in the source code, but couldn't find an obvious reason for why this is happening. I thought it might be that the G_database_units_to_meters function might be returning the wrong value, but I haven't actually tried to debug so I'm not ready to say that's it.

I do not see the same behavior in a UTM location (where units are in meters).

-------------------------------------------- Managed by Request Tracker

H_B · January 30, 2007, 2:16am

Request Tracker wrote:

this bug's URL: http://intevation.de/rt/webrt?serial_num=5454
---------------------------------------------------------------------

Subject: ps.map scaling not correct when projection units are feet

Platform: GNU/Linux/x86
grass obtained from: CVS
grass binary for platform: Compiled from Sources
GRASS Version: CVS HEAD 20070127 at 04:18 UTC

Noticed a strange anomaly in ps.map output when working in a State
Plane Coordinate System with coordinates in feet. Specifying any
scale in ps.map gives manifestly incorrect scaling. The scalebar
produced also bears no obvious relationship to the scale of the map.

To reproduce: create a location in New Mexico State Plane
Coordinates, Central Zone, NAD83 (EPSG number 2258). The result will
have a PROJ_INFO file like this: name: Transverse Mercator
proj: tmerc
datum: nad83
towgs84: 0.000,0.000,0.000
ellps: grs80
lat_0: 31
lon_0: -106.25
k: 0.999900
x_0: 500000.0001016001
y_0: 0
no_defs: defined

and a PROJ_UNITS file of:
unit: US survey foot
units: US survey foots

foots? ^^^^^^^

meters: 0.3048006096012192

Set the region to something arbitrary, for example:

> g.region -p
projection: 99 (Transverse Mercator)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1485814.68884
south: 1485263.50884
west: 1525396.34529
east: 1526171.30805
nsres: 0.61310345
ewres: 0.61310345
rows: 899
cols: 1264
cells: 1136336

Now use the following ps.map file to generate postscript:

paper a0
end
scalebar f
where 1.5 1.5
length 500
height .25
segment 5
numbers 2
fontsize 12
end
mapinfo
  where 11.5 15
  fontsize 12
end
grid 100
  color black
  numbers 2
end

note you need one more "end" to finish the command list.

I moved scalebar & info into the middle of the page:
paper a0
end
scalebar f
where 15.5 25.5
length 500
height .25
segment 5
numbers 1
fontsize 12
end
mapinfo
where 11.5 27
fontsize 12
end
grid 100
color black
numbers 2
end
end

This should display a 100 foot grid and nothing else, and should have
a scale bar with 100 foot markings that match the grid size.

yep.

For the region specified above, ps.map will select "1:982" as the
scale. Now try tinkering with the scale by adding scale commands to
the ps.map file. Even selecting "1:982" will give completely
different results than leaving the scale out altogether. As you
change scale, the scalebar will not change size to reflect the change
in the map grid.

G63> g.region -e
north-south extent: 551.180000
east-west extent: 774.962760
# looks ok (feet)

A0 paper is 33.07" x 46.77", uses a 1" margin (ps/ps.map/paper.h)

so for width, 31.07" => 774.962760' (9299.55") which is 1:299 :-/

1m / 1' = 3.28

1:299 * 3.28= 1:982 !

So it's a units problem, it's assuming meters when it is really feet.

I also see the explicit scale=1:982 making something different than
the automatic "1:982" when units=feet.
explicit scale=1:982 on the command line makes the graph box correct
but the scalebar is wrong. If 1:982 is figured automatically, both
scalebar and map box are wrong.

I dug around a little in the source code, but couldn't find an obvious
reason for why this is happening. I thought it might be that the
G_database_units_to_meters function might be returning the wrong
value, but I haven't actually tried to debug so I'm not ready to say
that's it.

I do not see the same behavior in a UTM location (where units are in
meters).

checks out for me too. ps.map works perfectly when units=meters.

more soon.

Hamish

H_B · January 30, 2007, 6:59am

this bug's URL: http://intevation.de/rt/webrt?serial_num=5454
---------------------------------------------------------------------

Subject: ps.map scaling not correct when projection units are feet

fixed in 6.3-CVS and 6.2 release branch. Also fixes bug #3096.

outstanding issue, due to :

> > and a PROJ_UNITS file of:
> > unit: US survey foot
> > units: US survey foots

Hamish:

> foots? ^^^^^^^

Tom:

That's what it picked on its own. I think the function that chooses
the plural just tacks an "s" on the end except under very specific
conditions.

lib/gis/proj2.c, proj3.c have some pluralizing magic, but in this case
that's the way it's stored in the PROJ_UNITS file created from EPSG code
#2258.

G63> g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
proj : tmerc
datum : nad83
ellps : grs80
lat_0 : 31
lon_0 : -106.25
k : 0.999900
x_0 : 500000.0001016001
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : US survey foot
units : US survey foots
meters : 0.3048006096012192

the only place I see "US survey foot" in the code is one line in
lib/proj/unit_of_measure.csv.

Create new locn from EPSG code is in general/g.proj/output.c,
but I don't really follow it.

The string gets passed to the ps.map mapinfo output, so it matters.
Work-around: edit the PERMANENT/PROJ_UNITS file by hand.

Hamish

ps- should we undelete lib/init/make_location_epsg.sh.in and move it into
tools/ as a useful utility script?

Glynn_Clements1 · January 30, 2007, 11:52am

Hamish wrote:

> this bug's URL: http://intevation.de/rt/webrt?serial_num=5454
> ---------------------------------------------------------------------
>
> Subject: ps.map scaling not correct when projection units are feet

fixed in 6.3-CVS and 6.2 release branch. Also fixes bug #3096.

outstanding issue, due to :

> > > and a PROJ_UNITS file of:
> > > unit: US survey foot
> > > units: US survey foots
Hamish:
> > foots? ^^^^^^^

I know that "US English" tends to prefer simply tacking an "s" on the
end to using irregular plurals (e.g. "datums" rather than "data"), but
that's taking it a bit far

Tom:
> That's what it picked on its own. I think the function that chooses
> the plural just tacks an "s" on the end except under very specific
> conditions.

lib/gis/proj2.c, proj3.c have some pluralizing magic, but in this case
that's the way it's stored in the PROJ_UNITS file created from EPSG code
#2258.

G63> g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
proj : tmerc
datum : nad83
ellps : grs80
lat_0 : 31
lon_0 : -106.25
k : 0.999900
x_0 : 500000.0001016001
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : US survey foot
units : US survey foots
meters : 0.3048006096012192

the only place I see "US survey foot" in the code is one line in
lib/proj/unit_of_measure.csv.

Create new locn from EPSG code is in general/g.proj/output.c,
but I don't really follow it.

AFAICT, the problem is with GPJ_osr_to_grass(); lib/proj/convert.c:

592 dfToMeters = OSRGetLinearUnits( hSRS, &pszUnitsName );
593
594 /* Workaround for the most obvious case when unit name is unknown */
595 if( (strcasecmp(pszUnitsName, "unknown") == 0) && (dfToMeters == 1.) )
596 G_asprintf( &pszUnitsName, "meter" );
597
598 G_set_key_value( "unit", pszUnitsName, *projunits );
599 sprintf( szFormatBuf, "%ss", pszUnitsName );
600 G_set_key_value( "units", szFormatBuf, *projunits );

AFAICT, the OGR API doesn't provide a mechanism to get the plural
name, so the code will need to handle it as a special case, e.g.:

  if (G_strcasecmp(pszUnitsName, "foot") == 0)
    strcpy( szFormatBuf, "feet" );
  else
    sprintf( szFormatBuf, "%ss", pszUnitsName );

[Do we also need to handle "inchs"? ]

Alternatively, it could use general/g.setproj/proj-units.table, which
has the singular and plural forms for all of the likely units (and
some rather unlikely ones, e.g. decinanometers).

--
Glynn Clements <glynn@gclements.plus.com>

Paul_Kelly · January 30, 2007, 4:18pm

On Tue, 30 Jan 2007, Glynn Clements wrote:

Hamish wrote:

Create new locn from EPSG code is in general/g.proj/output.c,
but I don't really follow it.

AFAICT, the problem is with GPJ_osr_to_grass(); lib/proj/convert.c:

592 dfToMeters = OSRGetLinearUnits( hSRS, &pszUnitsName );
593
594 /* Workaround for the most obvious case when unit name is unknown */
595 if( (strcasecmp(pszUnitsName, "unknown") == 0) && (dfToMeters == 1.) )
596 G_asprintf( &pszUnitsName, "meter" );
597
598 G_set_key_value( "unit", pszUnitsName, *projunits );
599 sprintf( szFormatBuf, "%ss", pszUnitsName );
600 G_set_key_value( "units", szFormatBuf, *projunits );

AFAICT, the OGR API doesn't provide a mechanism to get the plural
name, so the code will need to handle it as a special case, e.g.:

Yes that's exactly the reason - the WKT format used for import only contains the singular, and when I originally wrote that function I couldn't see a programmatically convenient way of accounting for irregular plurals. But...

  if (G_strcasecmp(pszUnitsName, "foot") == 0)
    strcpy( szFormatBuf, "feet" );
  else
    sprintf( szFormatBuf, "%ss", pszUnitsName );

I thought of something like this at the time, but that wouldn't have handled this particular case, where the full string was "US survey foot". I have also seen variations on this, e.g. "U.S. Survey Foot". But I had an idea which I think should be moderately robust - just do the comparison on the last 4 characters in the string.

It seems to have work and have committed it to CVS - also handles inches now

Out of interest, here's the two sets of PROJ_INFO/PROJ_UNITS files for a central Mexico State Plane location created using (a) EPSG code and g.proj and (b) selecting the stp projection and following the prompts in g.setproj:

(a) Created by g.setproj (State Plane)
-PROJ_INFO-------------------------------------------------
name : State Plane
datum : nad83
towgs84 : 0.000,0.000,0.000
proj : tmerc
a : 0.6378137e+07
es : 0.66943800229e-02
x_0 : 0.5e+06
y_0 : 0
k : 0.9999e+00
lon_0 : 106d15'w
lat_0 : 31dn
-PROJ_UNITS------------------------------------------------
unit : USfoot
units : USfeet
meters : 0.30480060960121920243

(b) Created by g.proj using EPSG code
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
proj : tmerc
datum : nad83
towgs84 : 0.000,0.000,0.000
ellps : grs80
lat_0 : 31
lon_0 : -106.25
k : 0.999900
x_0 : 500000.0001016001
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : US survey foot
units : US survey feet
meters : 0.3048006096012192

For (a) the "projection" key value in the WIND region file, raster headers etc. will be 2 (i.e. State Plane) while for (b) it will be 99 (i.e. Other). I think Maciek was asking recently about cases where this could be different yet the two projections equivalent, or something along those lines.

Paul