Skip to content

Transforming from cartographic to arbitrary local coordinate system

3 messages · Don MacQueen, Dylan Beaudette, Roger Bivand

Original
Don MacQueen
# 
I have virtually no experience with cartographic coordinate systems, 
and transforming them, so I would like to ask for some help.

Someone has supplied me with a shapefile. When I open it in Qgis, 
Qgis tells me its "spatial reference system" is:

   +proj=lcc +lat_1=35.4666666666667 +lat_2=34.0333333333333 
+lat_0=33.5 +lon_0=-118 +x_0=1999999.9371016 +y_0=499999.9843516001 
+ellps=GRS80 +to_meter=0.3048006 +no_defs

I can read it into R using
    tmp <- readShapeLines('filename',
                       proj4string=CRS(pstring))

where pstring is that thing that Qgis gave me.
plot(tmp) then gives a picture that makes sense.

I would like to give the thing a convenient (to me) local coordinate 
system. I have quite a few such objects (all with the same proj.4 
spatial reference string (or coordinate specification?). Some are 
lines objects, some are polygon objects, and the polygon objects have 
nested polygons ("holes" or "islands").

Pretend for the moment that the lines of my object represent borders 
of a rectangular shaped plot of land . The rectangle is tilted 
relative to north/south.

What I now would like to do is rotate and shift the rectangle so that 
its borders appear vertical and horizontal (my local axes are 
parallel to its borders), and so that it has a local origin near one 
of the corners. I guess this isn't a true cartographic projection (?).

Is there a way I can do this with spTransform(), by supplying an 
appropriate CRS argument?


Thanks
-Don
Dylan Beaudette
#
On Monday 19 November 2007, Don MacQueen wrote:
Hi Don, 

This looks like a US State-Plane coordinate system: Lambert Conformal Conic, 
with linear units defined in terms of feet. Try "re-projecting" the data to a 
local UTM zone first, or alternatively to achieve the rotation you are 
talking about, an Albers Equal Area projection. You will find the GDAL 
tool 'ogr2ogr' useful in this context. Here are some examples:

# convert shapefile whith defined projection to UTM Zone 10:
ogr2ogr -t_srs "+proj=utm +zone=10 +datum=NAD83" utm.shp original.shp 

Cheers,

Dylan
Roger Bivand
#
On Mon, 19 Nov 2007, Don MacQueen wrote:

            

      
      
      
    

        
The new elide methods in maptools would get you some of the way there, but 
only do 90 degree rotation, not arbitrary rotation. In general, sp classes 
expect north to be upwards, and no provision for other constructions is 
made. If you can find an appropriate PROJ.4 description for what you want, 
spTransform will do it, but the rotation is going to be the problem.

The easiest way in is almost certainly to customise elide() and adding 
trigenometry to do the arbitrary rotation.

Roger