I have a shapefile of polygons that I got from here, and a CSV file with longitude, latitude data (i.e. lon= -71.4; lat=42.5).

The shapefile has the following CRS information: "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

I can plot the shapefile and the location CSV file separately, but I have not been able to do it jointly.

Any idea what is going wrong in the code below?

shp <- readOGR('zipcodes_nt', 
               'ZIPCODES_NT_POLY', stringsAsFactors = FALSE)
coords = c(-71.42043,42.50242)
sp <- SpatialPoints(coords)
proj4string(sp) <- CRS("+init=epsg:4326") # Google use WGS84.
sp84 <- spTransform(sp, CRS("+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +
                                                 y_0=750000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
plot(sp84, col="red", add=TRUE) # do not add anything!
  • This code fails with an error at sp <- SpatialPoints(coords) - do you get that error? – Spacedman Mar 25 '19 at 11:06

You have an error in your code. SpatialPoints cannot be constructed from an atomic vector of coordinates. The following works fine:


              f <- tempfile(), mode='wb')
unzip(f, exdir=tempdir())
p <- readOGR(file.path(tempdir(), 'ZIPCODES_NT_POLY.shp'))

xy <- SpatialPoints(cbind(-71.42043,42.50242), CRS('+init=epsg:4326'))
xy2 <- spTransform(xy, CRS(proj4string(p)))

points(xy2, col='red', pch=20) 

enter image description here

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.