LoF10 LoF10 - 28 days ago 7
R Question

How to recode xy points into a new unit in a polygon spatial join?

So I'm looking to map xy points from a data frame and then recode them with a marker related to a new polygon so I can link them to that unit; from my understanding this requires a spatial join. I've been following these links:

R convert zipcode or lat/long to county

http://gis.stackexchange.com/questions/137621/join-spatial-point-data-to-polygons-in-r

I am missing something though regarding the conversion process of the dataframe and I'm not sure where I am making the error. I'm pasting in quick and dirty reproducible code from the open data API as an example. The shape file I am using is the community districts NYC file found here:

http://www1.nyc.gov/site/planning/data-maps/open-data/districts-download-metadata.page

code is below:

In this case as a reproducible sample lets take these points from open data and map them to the community districts file for NYC. Appreciate the help guys!

#libraries--------------------------

library(ggplot2)
library(ggmap)
library(sp)
library(jsonlite)
library(RJSONIO)

#call api data--------------------------

df = fromJSON("https://data.cityofnewyork.us/resource/24t3-xqyv.json?$query= SELECT Lat, Long_")
df = df[1:10000]
df <- data.frame(t(matrix(unlist(df), nrow=length(unlist(df[1])))))
names(df)[names(df) == 'X2'] = 'x'
names(df)[names(df) == 'X1'] = 'y'
df = df[, c("x", "y")]
df$x = as.numeric(as.character(df$x))
df$y = as.numeric(as.character(df$y))
df$x = round(df$x, 4)
df$y = round(df$y, 4)
df$x[df$x < -74.2] = NA
df$y[df$y < 40.5] = NA
df = na.omit(df)

#map data----------------------------




cd = readOGR("nycd.shp", layer = "nycd")
cd = spTransform(cd, CRS("+proj=longlat +datum=WGS84"))
cd_f = fortify(cd)


#map data
nyc = ggplot() +
geom_polygon(aes(x=long,
y=lat, group=group), fill='grey',
size=.2,color='black', data=cd_f, alpha=1)


nyc + geom_point(aes(x = x, y = y), data = df, size = 1)


#spatial join---------------------------


# Convert pointsDF to a SpatialPoints object
pointsSP <- SpatialPoints(df,
proj4string=CRS("+proj=longlat +datum=wgs84"))

Error in CRS("+proj=longlat +datum=wgs84") :
unknown elliptical parameter name

# Use 'over' to get _indices_ of the Polygons object containing each point
indices <- over(pointsSP, cd)

Answer

You have to provide a correct CRS. I found it on spatialreference.org (WGS 1984)

# Convert pointsDF to a SpatialPoints object 
pointsSP <- SpatialPoints(df, 
                          proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

After that, you have to reproject the points into the CRS of cd:

pointsSP <- spTransform(pointsSP, proj4string(cd))

# Use 'over' to get _indices_ of the Polygons object  containing each point 
indices <- over(pointsSP, cd)
head(indices)
#   BoroCD Shape_Leng Shape_Area
# 1    407  139168.79  328355685
# 2    407  139168.79  328355685
# 3    407  139168.79  328355685
# 4    407  139168.79  328355685
# 5    407  139168.79  328355685
# 6    105   35287.62   43796717
Comments