3

I am using rgdal to render some shapefile. The first is a set of building zones from NYC Open Data. The second is the NYC coast line from OSM.

The Zone files seems to generally render an X axis between 850000 and 1150000. And a Y axis with numbers from 150000 and 250000. The OSM coastline renders an X-axis between -8400000 and -8100000 (note the difference in order of magnitude as well as being negative). The OSMs Y axis is between 4900000 and 5050000.

I'm using this as my example and I don't the author doing anything special to format the axes.

My code is here:

# render a plot based of shapes in nycgiszoningfeatures_201311_shp/
# shapes downloaded from https://nycopendata.socrata.com/api/geospatial/kdig-pewd?method=export&format=Shapefile
require(rgeos)
require(maptools)
require(rgdal)

# OSM Coastline of NYC Data from http://osm-extracted-metros.s3.amazonaws.com/new-york.coastline.zip
mapBase <- readOGR('new-york-city.coastline', 'new-york')

shapeDir <- 'nycgiszoningfeatures_201311_shp'
shapeFiles <- dir(shapeDir, '\\.shp$')
shapeLayers <- sub('\\.shp$', '', shapeFiles)
# to see the files sprintf("%s/%s", shapeDir, shapeFiles)
# TODO: I'm told for loops are evil in R. Find a better way
for (i in 1:length(shapeLayers) ) {
    map <- readOGR(shapeDir, shapeLayers[i])
    #plot(mapBase, col="khaki2",border="grey", axes=TRUE, xlim=c(-8300000,-8200000), ylim=c(4900000, 5050000))
    #plot(map, add=TRUE, col="gray",border="blue")
    plot(map, col="gray",border="blue", axes=TRUE)
}

How do I get these shapes to plot with longitude and latitude numbers for the axes instead of those numbers.

3

Your shapefiles have the following projections:

R> proj4string(map)
[1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"
R> proj4string(mapBase)
[1] "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Reproject your data to WGS84 with spTransform (package sp):

crs <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

i <- 3
map <- readOGR(shapeDir, shapeLayers[i])
mapBase <- readOGR('new-york.shp', 'new-york')

plot(spTransform(mapBase, crs), axes=TRUE)
plot(spTransform(map, crs), add=TRUE, border="red")

plot

| improve this answer | |
  • aka the NYC data is using the NAD 1983 State Plane New York Long Island coordinate reference system. – mkennedy Jan 8 '14 at 16:52
  • @mkennedy wow, so basically before computers, there were geographically local coordinate systems that assumed the earth was flat to make the math easier? That's amazing to think about today. – Justin Dearing Jan 8 '14 at 20:22
  • @rcs that worked but the renderings don't line up. – Justin Dearing Jan 9 '14 at 3:22

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.