Search code examples
rgisshapefile

Looking up polygons in a shapefile that point belong in... why does it work with some shapefiles, not with others?

I am trying to look up census tract info by lat/lon coordinates using the census.gov shapefile located here:

This has been working quite successfully with code proposed by @lukeA for the ZCTA shapefile. See the original question and resolution here:

However, when I changed the shapefile to census tract, suddenly the lookup functionality broke, and I can't seem to figure out why. The documentation on the shapefile indicates that it is using NAD83 projection, which is the same as the ZCTA shapefile. The code doesn't return any errors, per se... but the output is all NA.

Here is the code:

library(maptools)
library(maps)
library(sp)
library(rgdal)

mn.zip.map <- readShapePoly("tl_2015_27_tract.shp")
proj4string(mn.zip.map) <- CRS("+proj=utm +zone=15 +datum=NAD83")
mn.zip.map <- spTransform(mn.zip.map, CRS("+proj=longlat"))
latlong <- data.frame(matrix(0,nrow=3,ncol=1))
latlong$lat <- c(44.730178, 44.784711, 44.943008)
latlong$long <- c(-93.235381, -93.476415, -93.303201)
coordinates(latlong) = ~long+lat
proj4string(latlong) <- CRS(proj4string(mn.zip.map))
over(latlong, mn.zip.map)

With output:

STATEFP COUNTYFP TRACTCE GEOID NAME NAMELSAD MTFCC FUNCSTAT ALAND AWATER INTPTLAT INTPTLON
1    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>
2    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>
3    <NA>     <NA>    <NA>  <NA> <NA>     <NA>  <NA>     <NA>    NA     NA     <NA>     <NA>

Again, clearly I am missing some nuance of GIS and geocoding.

Solution

  • You are doing something wrong with the projections. Use rgdal::readOGR() instead of maptools::readShapePoly(), this handles the projection information automatically:

    library("rgdal")
mn.zip.map <- readOGR(".", "tl_2015_27_tract")
proj4string(mn.zip.map)
# [1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
latlong <- data.frame(lat=c(44.730178, 44.784711, 44.943008),
                      long=c(-93.235381, -93.476415, -93.303201))
coordinates(latlong) <- ~long+lat
proj4string(latlong) <- CRS(proj4string(mn.zip.map))
over(latlong, mn.zip.map)
#   STATEFP COUNTYFP TRACTCE       GEOID   NAME            NAMELSAD MTFCC
# 1      27      037  060812 27037060812 608.12 Census Tract 608.12 G5020
# 2      27      139  080301 27139080301 803.01 Census Tract 803.01 G5020
# 3      27      053  108000 27053108000   1080   Census Tract 1080 G5020
#   FUNCSTAT    ALAND  AWATER    INTPTLAT     INTPTLON
# 1        S  3654009  299711 +44.7246681 -093.2316029
# 2        S 29771331 2237086 +44.7886836 -093.4454977
# 3        S   678381       0 +44.9444245 -093.3014752

plot(mn.zip.map, axes=TRUE)
plot(latlong, add=TRUE, col="red")

    plot