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

[Waichler, Scott R]

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? 

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