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

Waichler, Scott R Scott.Waichler at pnnl.gov
Tue Sep 17 18:52:41 CEST 2013

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