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, 2019 at 11:06

1 Answer 1


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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