[GRASS-user] Correctly Converting LL to LCC

GRASS GRASS user
1

   I apologize in advance for cross-posting this message, but I need to
quickly resolve this issue so I can proceed with the project. I posted this
to the proj4j mail list Friday (yes, I know it was a holiday), but have seen
no traffic on that list at all, and no responses to this request. Perhaps
someone here can help me.

   There are three points whose geographic coordinates I have in long/lat. I
need to convert these to the Oregon Lambert Conformal Conic projection with
the Washington-Oregon adjustment. When I try doing this with the cs2cs tool I
get incorrect values as a result.

   This is the input file:

#longitude|latitude
122:30:32.43W|45:19:19.49N
122:30:17.92W|45:18:52.45N
122:29:34.08W|45:18:47.16N

   This is the command line I used:

cs2cs +proj=latlong +datum=NAD83 +to +proj=lcc +datum=NAD83 +ellps=GRS80 \
+lat_1=43.0 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=1312336 +y_0=0 +nadgrids=WO

The above values are taken from the metadata for the Oregon LCC:

Lambert Conformal Conic
    stdparll 43.000000
    stdparll 45.500000
    long center -120.500000
    latprj origin 41.750000
    false east 1312335.958000 (International feet)
    false north 0.000000

    North American Datum of 1983
    Geodetic Reference System 80
    semiaxis 6378137.000000
    denflat 298.257222

   And these are the results:

#Easting|Northing
1154847.51|398797.24
1155143.02|397955.05
1156093.62|397768.55

   Not only do these seem reversed, they are not at all within the bounds of
the other project maps.

   I have the middle location on a separate map that was imported to GRASS
from an ESRI file, and its coordinates are approximately 796758.58,
1305792.74. Those values are in easting, northing order (at least, that's
the way I presume the GRASS-6.4 wxPython display presents them), and differ
in magnitude from what cs2cs produced.

   What have I done incorrectly? How should I properly translate the long/lat
values to lcc for my project?

Rich

_______________________________________________
Proj4j mailing list
Proj4j@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/proj4j

2

Can you setup a lat/long in GRASS, import lat/long points, then reproject within GRASS to your LCC projection?

Mark

On Dec 28, 2009, at 10:46 AM, Rich Shepard <rshepard@appl-ecosys.com> wrote:

I apologize in advance for cross-posting this message, but I need to
quickly resolve this issue so I can proceed with the project. I posted this
to the proj4j mail list Friday (yes, I know it was a holiday), but have seen
no traffic on that list at all, and no responses to this request. Perhaps
someone here can help me.

There are three points whose geographic coordinates I have in long/lat. I
need to convert these to the Oregon Lambert Conformal Conic projection with
the Washington-Oregon adjustment. When I try doing this with the cs2cs tool I
get incorrect values as a result.

This is the input file:

#longitude|latitude
122:30:32.43W|45:19:19.49N
122:30:17.92W|45:18:52.45N
122:29:34.08W|45:18:47.16N

This is the command line I used:

cs2cs +proj=latlong +datum=NAD83 +to +proj=lcc +datum=NAD83 +ellps=GRS80 \
+lat_1=43.0 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=1312336 +y_0=0 +nadgrids=WO

The above values are taken from the metadata for the Oregon LCC:

Lambert Conformal Conic
  stdparll 43.000000
  stdparll 45.500000
  long center -120.500000
  latprj origin 41.750000
  false east 1312335.958000 (International feet)
  false north 0.000000

  North American Datum of 1983
  Geodetic Reference System 80
  semiaxis 6378137.000000
  denflat 298.257222

And these are the results:

#Easting|Northing
1154847.51|398797.24
1155143.02|397955.05
1156093.62|397768.55

Not only do these seem reversed, they are not at all within the bounds of
the other project maps.

I have the middle location on a separate map that was imported to GRASS
from an ESRI file, and its coordinates are approximately 796758.58,
1305792.74. Those values are in easting, northing order (at least, that's
the way I presume the GRASS-6.4 wxPython display presents them), and differ
in magnitude from what cs2cs produced.

What have I done incorrectly? How should I properly translate the long/lat
values to lcc for my project?

Rich

_______________________________________________
Proj4j mailing list
Proj4j@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/proj4j
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

On Mon, 28 Dec 2009, MS wrote:

Can you setup a lat/long in GRASS, import lat/long points, then reproject
within GRASS to your LCC projection?

Mark,

   I don't see why not. I'd still like to understand why the cs2cs didn't
work as expected.

Rich

4

On Mon, 28 Dec 2009, MS wrote:

Can you setup a lat/long in GRASS, import lat/long points, then reproject within GRASS to your LCC projection?

Mark, et al.:

   Issues to be resolved:

   1) In both -6.4svn and -6.5svn using the wxPython GUI to create the new
location using 'll' as the projection, the spin control widgets are not
implemented for central latitude and central longitude. They will not change
from 0.0.

   2) Also with both versions, trying to set the extents/bounds fails. I
enter the appropriate values in each text entry widget, but clicking on the
"OK" (or whatever that button is labeled) does nothing. No reaction. Must
click the "Cancel" button to close that window.

   Manually entering the central latitude and logitude in
<location>/PERMANENT/PROJ_INFO, I used the lat_0 and lon_0 used in the
PROJ_INFO file of another location. Then I started GRASS, and tried to
import an ASCII point file.

   3) Copying the source file (sites.txt) to <location>/PERMANENT, I try to
import the points. This fails. See attached screenshot.

Thanks,

Rich

(attachments)

screenshot.png

5

Rich Shepard wrote:

   There are three points whose geographic coordinates I have in long/lat. I
need to convert these to the Oregon Lambert Conformal Conic projection with
the Washington-Oregon adjustment. When I try doing this with the cs2cs tool I
get incorrect values as a result.

   This is the command line I used:

cs2cs ... +x_0=1312336

Versus:

    false east 1312335.958000 (International feet)

Versus:

  http://trac.osgeo.org/proj/wiki/GenParms#Units

  If one uses the internal ellipsoid figure table, ellps=, then because
  all a= entries in the table are in meters the user must refer to the
  other linear unit items in meters.

   What have I done incorrectly? How should I properly translate the long/lat
values to lcc for my project?

You forgot to convert between feet and metres.

Technically, PROJ is units-agnostic. All of the values which it deals
with are dimensionless quantities. HOWEVER: you must use the same
units for everything. As the named ellipsoids use metres for the
semi-major axis, x_0 must also be in metres (the false easting given
above is 400000m), and you need to use +to_meter= if you want the
resulting cartographic coordinates converted to feet.

Looking for similar projections in the EPSG table (/usr/share/proj/epsg),
I find:

# NAD83 / Oregon Lambert
<2991> +proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs <>
# NAD83 / Oregon Lambert (ft)
<2992> +proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048 +no_defs <>
# NAD83(HARN) / Oregon Lambert
<2993> +proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs <>
# NAD83(HARN) / Oregon Lambert (ft)
<2994> +proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999984 +y_0=0 +ellps=GRS80 +to_meter=0.3048 +no_defs <>

Note that all of these use +x_0=400000 (give or take 1.6 microns), and
those labelled "(ft)" have "+to_meter=0.3048".

Aside from issues related to PROJ itself, make sure that GRASS knows
that your coordinates are in feet (the PROJ_UNITS file should contain
"meters: 0.3048" (or "meters: 0.30480060960121920243" if you want US
survey feet).

Also, take care when dealing with elevation data. Modules will
typically assume that elevations are in metres unless told otherwise
(e.g. zmult=0.3048 for r.slope.aspect), but I wouldn't entirely rule
out the possibility of something somewhere assuming that elevations
use the same units as the horizontal coordinates. Check that the
results make sense.

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

6

Mark:

> Can you setup a lat/long in GRASS, import lat/long
> points, then reproject within GRASS to your LCC projection?

Rich wrote:

  Issues to be resolved:

  1) In both -6.4svn and -6.5svn using the wxPython
GUI to create the new
location using 'll' as the projection, the spin control
widgets are not
implemented for central latitude and central longitude.
They will not change from 0.0.

They are locked at 0,0 because you typically want 0 lat, 0 lon to be
somewhere off the west coast of Africa. don't worry about it.

They are shown because some other projections with locked terms are more
interesting..

  2) Also with both versions, trying to set the
extents/bounds fails. I
enter the appropriate values in each text entry widget, but
clicking on the
"OK" (or whatever that button is labeled) does nothing. No
reaction. Must click the "Cancel" button to close that window.

I think that will be the wxGUI location wizard still needing some work.
there's a trac ticket for it. The work-around is to enter the bounds
as decimal degrees not ddd:mm:ss.sssH. from the old text based "wizard"
you can enter in either format.

  Manually entering the central latitude and logitude in
<location>/PERMANENT/PROJ_INFO, I used the lat_0 and
lon_0 used in the PROJ_INFO file of another location.

don't do that.

Then I started GRASS, and tried to import an ASCII point file.

  3) Copying the source file (sites.txt) to
<location>/PERMANENT, I try to
import the points. This fails. See attached screenshot.

can you cut and paste the whole command line? it gets cut off.

what does `head sites.txt` look like?

Hamish

7

Rich wrote:

  I apologize in advance for
cross-posting this message, but I need to
quickly resolve this issue so I can proceed with the
project. I posted this
to the proj4j mail list Friday (yes, I know it was a
holiday), but have seen
no traffic on that list at all, and no responses to this
request.

I don't think it made it.
  http://lists.maptools.org/pipermail/proj/2009-December/date.html

ah,

_______________________________________________
Proj4j mailing list
Proj4j@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/proj4j

that's not the one listed on the PROJ.4 homepage, that's the mailing list
for a java port of proj.4.
http://lists.osgeo.org/pipermail/proj4j/2009-November/000000.html
http://trac.osgeo.org/proj4j/
https://trac.osgeo.org/proj4j/ticket/13

  This is the command line I used:

cs2cs +proj=latlong

just as a matter of style I always like to use +proj=longlat there, it
reminds me that the input order should be "x y" not "lat long".

Those values are in easting, northing order (at least, that's
the way I presume the GRASS-6.4 wxPython display presents them),

yup, easting,northing.

Hamish

8

On Mon, 28 Dec 2009, Hamish wrote:

They are locked at 0,0 because you typically want 0 lat, 0 lon to be
somewhere off the west coast of Africa. don't worry about it.

Hamish,

   Ah, yes! This is in Lat/Long, not a projection. Fixed it.

I think that will be the wxGUI location wizard still needing some work.
there's a trac ticket for it. The work-around is to enter the bounds as
decimal degrees not ddd:mm:ss.sssH. from the old text based "wizard" you
can enter in either format.

   OK.

can you cut and paste the whole command line? it gets cut off.

   Surprisingly, I cannot copy from the GUI window. The only missing text is
the full name of the output file: features

what does `head sites.txt` look like?

   It looks like the entire file:

#longitude|latitude|name
122:30:32.43W|45:19:19.49N|Feature 1 Name
122:30:17.92W|45:18:52.45N|Feature 2 Name
122:29:34.08W|45:18:47.16N|Feature 3 Name

Thanks,

Rich

9

On Mon, 28 Dec 2009, Glynn Clements wrote:

cs2cs ... +x_0=1312336

Versus:

    false east 1312335.958000 (International feet)

Versus:
  http://trac.osgeo.org/proj/wiki/GenParms#Units
You forgot to convert between feet and metres.

Glynn,

   Thank you. This is a point of confusion for me because the metadata from
the source shows the false easting in International feet, yet it states that
units should be in meters. I've never understood why.

Technically, PROJ is units-agnostic. All of the values which it deals with
are dimensionless quantities. HOWEVER: you must use the same units for
everything. As the named ellipsoids use metres for the semi-major axis,
x_0 must also be in metres (the false easting given above is 400000m), and
you need to use +to_meter= if you want the resulting cartographic
coordinates converted to feet.

   No, I want them in meters so everything matches.

Looking for similar projections in the EPSG table (/usr/share/proj/epsg),
I find:

# NAD83 / Oregon Lambert
<2991> +proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5
+x_0=400000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs <>

   The above is what I do want.

Thank you,

Rich

10

On Mon, 28 Dec 2009, Hamish wrote:

that's not the one listed on the PROJ.4 homepage, that's the mailing list
for a java port of proj.4.

   Oops! Guess I saw the wrong one on the list of lists. Wondered about that
final 'j'. Thanks.

cs2cs +proj=latlong

just as a matter of style I always like to use +proj=longlat there, it
reminds me that the input order should be "x y" not "lat long".

   Good point.

Rich

11

On Mon, 28 Dec 2009, Glynn Clements wrote:

cs2cs ... +x_0=1312336

Versus:

    false east 1312335.958000 (International feet)

Versus:
  If one uses the internal ellipsoid figure table, ellps=, then because
  all a= entries in the table are in meters the user must refer to the
  other linear unit items in meters.

   I still have errors in that command line. Changing to +x_0=400000
produces:

[rshepard@salmo /usr4/keypoints]$ cs2cs +proj=longlat +datum=NAD83 +to +proj=lcc +datum=NAD83 +ellps=GRS80 \

+lat_1=43.0 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 \
+nadgrids=WO <<EOF
122d30'32.43"W 45d19'19.49"N
122d30'17.92"W 45d18'52.45"N
122d29'34.08"W 45d18'47.16"N
EOF

242511.51 398797.24 -0.00
242807.02 397955.05 -0.00
243757.62 397768.55 -0.00

and these are still off. I must still have parameters incorrect.

Rich

12

Hamish:

> can you cut and paste the whole command line? it gets cut off.

Rich wrote:

  Surprisingly, I cannot copy from the GUI window.

try right-click + copy from the context menu, or a copy button.

then ^V to paste instead of the middle click.

I'm not sure if that's entirely by design or just the way it's working so
far.

the module gui box should have a separate copy button for the command.

re. dbmi error:
did you mess with db.connect at all or are the DB settings original?
what does "db.connect -p" say?

Hamish

13

On Mon, 28 Dec 2009, Hamish wrote:

  Surprisingly, I cannot copy from the GUI window.

try right-click + copy from the context menu, or a copy button.
then ^V to paste instead of the middle click.

Hamish,

   Next time. Apparently, correcting the lat_0 and lon_0 cleared the error.

re. dbmi error: did you mess with db.connect at all or are the DB settings
original? what does "db.connect -p" say?

   I've ignored the db aspect for now. One step at a time.

Thanks,

Rich

14

On Mon, 28 Dec 2009, Rich Shepard wrote:

+lat_1=43.0 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 \
+nadgrids=WO <<EOF
122d30'32.43"W 45d19'19.49"N
122d30'17.92"W 45d18'52.45"N
122d29'34.08"W 45d18'47.16"N
EOF

242511.51 398797.24 -0.00
242807.02 397955.05 -0.00
243757.62 397768.55 -0.00

   Please excuse any of my frustration that I inadvertently express here, but
I just cannot get these three points correctly into GRASS so they display
with the other maps in the project area.

   I used the GUI to import the long/lat values. Those appear OK (but I've
not checked the decimal degrees with the original DMS input). However, when
I run 'v.proj input=features location=sitesLL mapset=PERMANENT
output=keysites --o' the same incorrect positions as above are produced.
Therefore, when I try to display that map with the basin boundary and
streams the error message is that those coordinates are outside the bounds
of the current region.

   The middle point should have values of 796758.58 (instead of 242807.02)
and 130579274 (instead of 397955.05).

   It must be something quite simple that I keep missing, but I cannot find
where I've gone wrong.

Rich

15

Rich wrote:

  The middle point should have values of 796758.58
(instead of 242807.02)
and 130579274 (instead of 397955.05).

it's a feet to meters thing.

242852.015184 / 0.3048
ans = 796758.58

check for correct units with "g.proj -p" from both locations.

or use v.transform to scale things and quietly walk away..

Hamish

16

On Mon, 28 Dec 2009, Hamish wrote:

it's a feet to meters thing.

242852.015184 / 0.3048
ans = 796758.58

   That's clear.

   For future reference, how can I avoid this? I see no unit option in v.proj
or cs2cs that allows me to specify meters when the long/lat degrees are
projected.

check for correct units with "g.proj -p" from both locations.

   For the source (location sitesLL) the units are degrees, for the target
(location Oregon) the units are degrees.

or use v.transform to scale things and quietly walk away..

   Think I did this incorrectly. I started GRASS in the Oregon location,
mapset project_name. Killed the GUI and ran 'v.transform -s input=keysites
output=projpts xscale=0.3048 yscale=0.3048 --o', then re-started the GUI,
loaded the output map and saw that the coordinates of a displayed point had
not changed. Ergo, I must not have properly expressed the scaling factor.

Rich

17

Rich Shepard wrote:

>> +lat_1=43.0 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 \
>> +nadgrids=WO <<EOF
>> 122d30'32.43"W 45d19'19.49"N
>> 122d30'17.92"W 45d18'52.45"N
>> 122d29'34.08"W 45d18'47.16"N
>> EOF
> 242511.51 398797.24 -0.00
> 242807.02 397955.05 -0.00
> 243757.62 397768.55 -0.00

   Please excuse any of my frustration that I inadvertently express here, but
I just cannot get these three points correctly into GRASS so they display
with the other maps in the project area.

   I used the GUI to import the long/lat values. Those appear OK (but I've
not checked the decimal degrees with the original DMS input). However, when
I run 'v.proj input=features location=sitesLL mapset=PERMANENT
output=keysites --o' the same incorrect positions as above are produced.
Therefore, when I try to display that map with the basin boundary and
streams the error message is that those coordinates are outside the bounds
of the current region.

   The middle point should have values of 796758.58 (instead of 242807.02)
and 130579274 (instead of 397955.05).

  242807.02 / 796758.58 = 0.3047435272049408
  397955.05 / 1305792.74 = 0.3047612670905185

These look suspiciously close to 0.3048 (foot/metre conversion).

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

18

On Tue, 29 Dec 2009, Glynn Clements wrote:

  242807.02 / 796758.58 = 0.3047435272049408
  397955.05 / 1305792.74 = 0.3047612670905185

These look suspiciously close to 0.3048 (foot/metre conversion).

Glynn,

   You and Hamish point out that the problem is with conversion of the units:
feet vs. meters. What I'm trying to understand is why there is this problem
and what to do about it.

   The imported Long/Lat location is PROJ_UNITS of degrees, which is what I
would expect. The location with the other maps has PROJ_UNITS of meters.
Therefore, when I run v.proj to convert the ll map coordinates in degrees to
lcc coordinates in meters, why do I end up with feet instead? I see no
length unit option in either v.proj or cs2cs. My assumption is user error,
and I'd like to learn what error I made so I don't repeat it in the future.

   Hamish suggested I scale the units using v.transform. I tried this with
both xscale and yscale at 0.3048 as well as 1.0/0.3048 but neither resulted
in units within the appropriate range.

   It may be that all my maps from the state GIS data repository are
incorrect because soils data from a different agency also cannot be
overlaid. Since I used the ArcView .prj and Arc/Info .e00 files for the
projection and boundary information I need to learn what I've done
incorrectly so I can re-import all these data and be able to move ahead.

   If someone will work with me off the list to get straightened out I would
greatly appreciate it. I have a 2-page .pdf showing the Oregon data library
formats and projection information. Since these data are in International
feet, I need to understand why the PROJ_UNITS are in meters. With someone's
gracious help I can get all this done today and return to being a happy
camper.

Rich

19

Rich Shepard wrote:

> 242807.02 / 796758.58 = 0.3047435272049408
> 397955.05 / 1305792.74 = 0.3047612670905185
>
> These look suspiciously close to 0.3048 (foot/metre conversion).

   You and Hamish point out that the problem is with conversion of the units:
feet vs. meters. What I'm trying to understand is why there is this problem
and what to do about it.

   The imported Long/Lat location is PROJ_UNITS of degrees, which is what I
would expect. The location with the other maps has PROJ_UNITS of meters.
Therefore, when I run v.proj to convert the ll map coordinates in degrees to
lcc coordinates in meters, why do I end up with feet instead?

Do you end up with feet? Or do you end up with metres but you're
comparing against a reference given in feet? Your previous comment
suggests the latter:

   The middle point should have values of 796758.58 (instead of 242807.02)
and 130579274 (instead of 397955.05).

Your numbers seem roughly correct for metres (45d19'N is ~3.567 degrees
north of lat_0, which equates to ~397km), while the "desired" numbers
are roughly correct for feet (397000 / 0.3048 ~= 1,302,493 feet).

I see no length unit option in either v.proj or cs2cs.

v.proj uses the information in the PROJ_INFO and PROJ_UNITS files.
Whether a location uses metres or feet is a property of the location.
If you use metres when you create the location, all coordinates will be
given in (or assumed to be in) metres.

cs2cs and proj are units-agnostic. If you give the semi-major axis
(+a=... option) in metres, you get coordinates in metres (or, for
cartographic->lat/lon, the input is treated as metres). If you give
the semi-major axis in furlongs, you get coordinates in furlongs.

All of the built-in ellipsoids which can be selected by the +ellps=...
option use metres. If you want proj or cs2cs to provide or receive
coordinates in units other than those used when specifying the
semi-major axis, you need to use the +to_meter=... option to specify
the scale factor, e.g. +to_meter=0.3048 if you want units of 0.3048
metres (i.e. feet).

   Hamish suggested I scale the units using v.transform. I tried this with
both xscale and yscale at 0.3048 as well as 1.0/0.3048 but neither resulted
in units within the appropriate range.

What are you comparing against? How far off are the numbers?

   It may be that all my maps from the state GIS data repository are
incorrect because soils data from a different agency also cannot be
overlaid. Since I used the ArcView .prj and Arc/Info .e00 files for the
projection and boundary information I need to learn what I've done
incorrectly so I can re-import all these data and be able to move ahead.

   If someone will work with me off the list to get straightened out I would
greatly appreciate it. I have a 2-page .pdf showing the Oregon data library
formats and projection information. Since these data are in International
feet, I need to understand why the PROJ_UNITS are in meters.

This is determined when you create the location.

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

20

On Tue, 29 Dec 2009, Glynn Clements wrote:

Do you end up with feet? Or do you end up with metres but you're comparing
against a reference given in feet? Your previous comment suggests the
latter:

Glynn,

   Probably the latter.

All of the built-in ellipsoids which can be selected by the +ellps=...
option use metres. If you want proj or cs2cs to provide or receive
coordinates in units other than those used when specifying the semi-major
axis, you need to use the +to_meter=... option to specify the scale
factor, e.g. +to_meter=0.3048 if you want units of 0.3048 metres (i.e.
feet).

   This must be where I went off. The metadata tells me that the state uses
international feet so I need to have that translated to meters, or specify
in the location definition that the units are international feet.

This is determined when you create the location.

   Your advice, please. Should I use the +to_meter conversion or specify all
distance units in international feet?

Thanks,

Rich