[R-sig-Geo] Transforming from cartographic to arbitrary local coordinate system

[Dylan Beaudette]

On Monday 19 November 2007, Don MacQueen wrote:
> 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

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

-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341