1

For some research I am working on, I have a dataframe called charters of charter schools in New Jersey. Simplified it looks like this:

      school.name          school.id        district.id     lat     lon
oceanside charter school     008229           3400011     39.3635  -74.4350 
discovery charter school     228881           3400020     40.7343  -74.1745

I also have a SpatialPolygonsDataFrame of school districts in New Jersey which I got through the Tigris package: nj.school.districts <- tigris::school_districts("NJ"). For regular, non-charter schools, their district.id matches their geographic district.id, as one might expect. However, for charter schools a unique district.id is assigned that doesn't reflect their geographic location at all.

My goal is to assign a geographic.district.id to each school in charters that reflects the district it is actually located in. The nj.school.districts spdf looks like this:

enter image description here

So there are 342 districts, each with a unique GEOID. It is this GEOID that should be assigned to each charter school that falls within the district. I've mapped the charters over the school districts so could do it manually that way, but in reality it would be too time consuming.

I've been able to extract the coordinates of a particular polygon, and to confirm whether or not any schools are in it through this code:

nj.polygon1.coords <- as.data.frame(nj.school.districts@polygons[[1]]@Polygons[[1]]@coords) 
%>% select(lat = V2, lon = V1)

xp = as.vector(nj.polygon1.coords$lon)
yp = as.vector(nj.polygon1.coords$lat)

in.poly <- pracma::inpolygon(charters$lon, charters$lat, xp, yp)

but I don't know how to do it at scale or how to then assign the new GEOID accordingly.

Improve this question
3
  • What do you mean by doing it "at scale"?
    – Spacedman
    Sep 17, 2019 at 15:18
  • By "at scale" I mean that I want ultimately to this process not just for ~100 charter schools in NJ but actually every charter school in the US. So doing it manually is unreasonable.
    – Eli Groves
    Sep 17, 2019 at 15:44
  • The scaling limit will probably be the size and complexity of the polygons. I reckon I can do 100 schools in NJ in much less time than it takes to download the boundary data using tigris (CA is being very slow to download now...). If your school data is already broken down by state then looping over state should be fairly efficient - creating a single whole-US boundary file might get large.
    – Spacedman
    Sep 17, 2019 at 16:27

1 Answer 1

Reset to default
2

Let's reproduce your sample data:

charters=data.frame(
  school.name=c("oceanside","discovery"),  
  school.id=c("008229","228881"), 
  district.id=c("3400011","340020"), 
  lat=c(39.3635, 40.7343),
  lon=c(-74.4350,-74.1745))

then make it a spatial points data frame:

coordinates(charters)=~lon+lat

assume the lat-long in charters is the same sort of lat-long as in the NJ boundaries:

proj4string(charters) = proj4string(nj.school.districts)

Now if you do over of points on polygons you get back the polygons that each point is in:

over(charters, nj.school.districts)

so extract the field in the polygons and bung it in the points:

charters$GEOID = over(charters, nj.school.districts)$GEOID

giving:

> as.data.frame(charters)
  school.name school.id district.id     lat      lon   GEOID
1   oceanside    008229     3400011 39.3635 -74.4350 3400960
2   discovery    228881      340020 40.7343 -74.1745 3411340
Improve this answer
3
  • Thank you so much. In retrospect I was quite close to this solution but was lacking the proj4string() assignment. For clarity, this sets the 'type' of lat-long coordinates to be the same? Would it be accurate to say that it puts them on the same scale?
    – Eli Groves
    Sep 17, 2019 at 16:46
  • No, its saying that are the "same coordinate system". There are multiple latitude-longitude standards that differ because of updated measurements of the earth's shape, or moving an origin slightly. If possible you should get the PROJ string that describes the coordinate system of the data from the source of the data. If not stated its usually WGS84 coordinates, as used by GPS, with a PROJ string of "+init=epsg:4326".
    – Spacedman
    Sep 17, 2019 at 18:17
  • I see, that makes sense. Thank you for the clarification. Also, for those who may be unfamiliar, the over() function used here is sp::over()
    – Eli Groves
    Sep 18, 2019 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.