[R-sig-Geo] function to convert to/from U.S. State Plane Coordinates?

[Roger Bivand]

On Tue, 17 Sep 2013, Waichler, Scott R wrote:

> Erin,
>
> Thanks, now I see that I was trying to use spTransform() on a map object 
> instead of a Spatial one.  I followed your logic and came up with the 
> following to plot the state of Washington using state plane coordinates 
> instead of the original geographic values.  But I get an extraneous 
> triangular polygon.  I am guessing it has to do with not repeating the 
> initial x,y at the end of a polygon, or a missing row of NA.  What am I 
> missing?

You have not looked for information on converting legacy map objects to sp 
classes. Look for ?map2SpatialPolygons in the maptools package. Better 
still, find a better representation of the state boundary as a shapefile 
and use that instead.

Roger

>
> library(maps)
> library(rgdal)
> wa <- map('state', region = 'washington', plot=F, type="n")
> ind.not.na <- which(!is.na(wa$x))  # note which rows are not NA
> wa1 <- data.frame(x=wa$x,y=wa$y)  # make a Spatial object from the map object
> coordinates(wa1) <- ~x+y
> proj4string(wa1) <- CRS("+proj=longlat +datum=WGS84")
> # I got the CRS spec from http://spatialreference.org/ref/esri/102749/proj4/
> wa.spc <- spTransform(wa1, CRS("+proj=lcc +lat_1=45.83333333333334 +lat_2=47.33333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0000000002 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs"))
>
> # make another data frame with the same number of rows as the map() original
> wa2 <- data.frame(x=wa$x,y=wa$y)
> # repace the geographic values with the state plane coordinate values
> wa2$x[ind.not.na] <- wa.spc$x
> wa2$y[ind.not.na] <- wa.spc$y
> plot(wa2, type="n")
> polygon(wa2)  # gives a spurious triangular polygon--why?
>
>> Ok.  I think I see it now.
>>
>> You need to have a spatial object rather than a list.
>>
>> So you might set up
>>
>> 	pnw1 <- data.frame(x=pnw$x,y=pnw$y)
>> 	coordinates(pwn1) <- ~x+y
>> Then you'll have a spatial object.  After that, put in your projection:
>> 	proj4string(pnw1) <- CRS(xxxxx)
>> 	pnw2 <- spTransform(pnw1,CRS(yyyy))
>
> --Scott
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no