1

I have successfully loaded a shapefile of NYC PUMA areas into R with maptools and I want to plot 55 points on top of it that I have in another file as follows:

X   Y   pumace10    events_2008 events_2009
-73.9092456917  40.8916125162   3701    2   0
-73.8617096298  40.8899373255   3702    0   0
-73.8010284966  40.8460832277   3703    1   1

However, the points will not plot.

First I do this to plot the shapefile:

plot(nycs)

And it plots the shapefilethis is the shapefile

Then I try to plot the points on top but no matter which of the following I do it always fails:

points(nyc_data$X,nnyc_data$Y,pch=20,col="red")

or

plot(nyc_data, pch=16, col='firebrick',add=TRUE)

or

plot(nyc_data$X,nyc_data$Y,pch=20,col="red")

(that final one plots the data on a new plot that is just an X-Y scatter instead of overlaid on the shapefile)

Any ideas how to do this?

EDIT, files added (amended to working files, hopefully!):

Shapefile info: https://www.sendspace.com/file/wbqrpb Points file: https://www.sendspace.com/file/9yrrbu

the_t_test_1
  • 1,193
  • 1
  • 12
  • 28

1 Answers1

0

Anyway you could send the files you are using?

My guess is that you have one of two problems:

1) either the "points" spatial data frame or the NYC PUMA file do not have a coordinate reference system.

2) they both have reference systems, but they are different.

Likely, your problem is that the points lack a reference system.

EDIT:

The problem is that the two datasets are in different coordinates, and I'm not sure what the coordinate system for the points is. Here is my code. The problem is that the CRS for the points is clearly neither regular lat/long nor the projection used in the shapefile. If you have more info on the source of the data points maybe we can see what projection they use.

library(sp)
library(maptools)
library(rgdal)
library(rgeos)

Points=read.csv("nyc_data_sample.csv",stringsAsFactors=F)
shapefile=readOGR("shapefiles","nyu_2451_34512")



Points_Shape = SpatialPoints(Points[,c("X","Y")],proj4string=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
#Points_Shape = SpatialPoints(Points[,c("X","Y")],proj4string=CRS(proj4string(shapefile)))

Points_Shape = SpatialPointsDataFrame(Points_Shape,Points)
Mauricio Romero
  • 77
  • 9
  • I thought I had done this but now you've made me doubt it... I've uploaded the shape and points file above ^ – the_t_test_1 Mar 01 '17 at 09:24
  • Does your shapefile has a .prj , a .dbf file, etc? also, your .shp does not work. – Mauricio Romero Mar 01 '17 at 13:05
  • Yep - they are now attached also. – the_t_test_1 Mar 01 '17 at 15:19
  • not siure oif you read the second part of my comment. your .shp does not work. it doesnt open in either R or Qgis – Mauricio Romero Mar 01 '17 at 19:17
  • It works in both for me, and in R I loaded it via nycs=readShapePoly("/home/nyc/02_project/data/acs/shapefiles/nyu_2451_34512.shp",proj4string=crswgs84,verbose=TRUE) – the_t_test_1 Mar 02 '17 at 14:11
  • It doesnt work. I used the same command + tried readOGR. Likely the uploaded file does not match the one in your computer – Mauricio Romero Mar 02 '17 at 18:17
  • That's weird - I checked the uploaded files and you're right, sorry. No idea why... anyway, have updated the files above and they should work fine (I've tested). Thanks for your effort and hope this works :) – the_t_test_1 Mar 03 '17 at 10:44
  • Any advances on this? I've checked and the files work but still can't plot :/ – the_t_test_1 Mar 06 '17 at 11:56
  • I downloaded the data (it works now). The problem is that they come in different coordinates. Here is my code – Mauricio Romero Mar 06 '17 at 12:29
  • The points file are geocoded addresses through Bing (WGS84?) and the shapefile is from Baruch college CUNY (might be EPSG:2263. NAD83 / New York Long Island according to the .xml file?)... – the_t_test_1 Mar 06 '17 at 18:36
  • What code did you use to geocode the addresses? They are not in lat/long (look at the numbers, its just impossible) – Mauricio Romero Mar 07 '17 at 15:41
  • You're absolutely right - I'm genuinely terribly sorry for misleading you there... the points are actually generated from the centroids of the polygons through QGIS! Perhaps I need to change the settings by which I generate them? – the_t_test_1 Mar 07 '17 at 18:57
  • Likely, the problem is simply that you need to know what the projection used is... thats all. Currently, there are coordinates but without a reference system (aka, projection) they are meaningless – Mauricio Romero Mar 07 '17 at 20:19