Search code examples
rgisrgdalr-maptools

Add lat/long or URM coordinates to shapefiles in R

I am new to GIS in R, and I'm attempting to add either lat/long or UTM coordinates to a shapefile. I downloaded the boundaries of the city of Chicago (City_Boundaries.shp)from here: http://www.cityofchicago.org/city/en/depts/doit/supp_info/gis_data.html

I loaded the maptools and rgeos libraries:

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

I brought the data into R & attempted to add a UTM code for zone 16T:

city<- readShapeSpatial("C:/Users/Luke.Gerdes/Documents/GIS Files/Chicago/City_Boundary.shp")
proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")

However, the resulting data doesn't make sense to me. When I take a look at the "coords" slots within "city", the coordinate values are, for example, X=1092925 and Y=1944820. I used an outside tool (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html) to find GPS coordinates. The results were: lat=17.511, long=-81.42. This is somewhere between Jamaica and the coast of Honduras. It seems to me that the shapefile coordinates exist in their own universe; as the name shapefile perhaps implies, the coordinates provide an accurate depiction of the city's shape, but these coordinates do not automatically map onto the globe.

My ultimate goal is to determine if a number of geo-tagged events took place within Chicago. I am comfortable converting those points, which are listed by lat/long, into UTM, as per below:

  SP <- SpatialPoints(cbind(-87.63044, 41.79625), proj4string=CRS("+proj=longlat +datum=WGS84"))
  SP <- spTransform(SP, CRS("+proj=utm +north +zone=16T +datum=WGS84"))

I'm also willing to try to work with the events data in their original format if that makes more sense. Ultimately, I need help to turn the Chicago city boundaries into a polygon that I can check against my events data. Something like this:

  gContains(city, SP)

How do I get my shapefiles and events data into a common (and accurate) frame of reference so that I can check whether points are in the city? Any help you can render will be greatly appreciated.

Solution

  • The general advice is to read shapefiles with rgdal::readOGR, as in

    library(rgdal)
city = readOGR(".", "City_Boundary")

    That will not only give city its proper CRS:

    proj4string(city)
[1] "+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

    but also warn you if you want to change it without reprojecting:

    proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")
Warning message:
In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
  A new CRS was assigned to an object with an existing CRS:
+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform in package rgdal

    meaning that after doing this, city is wrong. You reproject city by using spTransform, as in

    city = spTransform(city, CRS("+proj=utm +north +zone=16T + datum=WGS84"))

    Your rgeos::gContains call might fail because the two CRS strings for utm differ; remove the space in + datum=WGS84 in the former and you'll be fine (TRUE).