Los ejes de los shapefiles renderizados no muestran los números que son de latitud y longitud

Question

Estoy usando rgdal para representar algún shapefile. El primero es un conjunto de zonas de construcción de NYC Open Data. El segundo es el Línea de costa de NYC desde OSM .

Los archivos de la Zona parecen renderizar generalmente un eje X entre 850000 y 1150000. Y un eje Y con números entre 150000 y 250000. La línea de costa OSM renderiza un eje X entre -8400000 y -8100000 (nótese la diferencia en el orden de magnitud además de ser negativo). El eje Y del OSM está entre 4900000 y 5050000.

Estoy usando este como mi ejemplo y el autor no hace nada especial para formatear los ejes.

Mi código está aquí:

# 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)
}

¿Cómo consigo que estas formas se tracen con números de longitud y latitud para los ejes en lugar de esos números?

Answer 1

Sus shapefiles tienen las siguientes proyecciones:

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"

Reproyecte sus datos a WGS84 con spTransform (paquete 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")