[GRASS-user] Correctly Converting LL to LCC

[Michael Barton]

On Dec 29, 2009, at 10:00 AM, grass-user-request at lists.osgeo.org wrote:

> Message: 8
> Date: Tue, 29 Dec 2009 15:42:08 +0000
> From: Glynn Clements <glynn at gclements.plus.com>
> Subject: Re: [GRASS-user] Correctly Converting LL to LCC
> To: Rich Shepard <rshepard at appl-ecosys.com>
> Cc: grass-users at lists.osgeo.org
> Message-ID: <19258.9040.993225.362218 at cerise.gclements.plus.com>
> Content-Type: text/plain; charset=us-ascii
> 
> 
> 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.

I don't want to point fingers at the 'other' GIS, but I and others have run into problems with ArcView and .e00 files before. If the original data were recorded incorrectly, ArcView doesn't enforce correct projections the way GRASS does and bad projection info can end up in these files. When imported into GRASS, of course, it assumes that they are correct. So it might be worth your while to look at the original data source before tearing your hair out dealing with the imported data. If there are problems with the source, you might be able to deal with by reimporting and setting your input parameters differently.

Michael