1

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

https://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)
Community
  • 1
  • 1
LoF10
  • 1,907
  • 1
  • 23
  • 64

1 Answers1

0

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
loki
  • 9,816
  • 7
  • 56
  • 82