R shapefile transform longitude latitude coordinates

Question

I want to find out the coordinates (latitude and longitude) of Texas precincts. I download the shapefile from here.

I am using R to read the file and extract the coordinates:

dsn.tx <- "~/.../Texas.shp"        
shp <- readOGR(dsn=dsn.tx, stringsAsFactor = F)
coords <- coordinates(shp)

But the resulting coordinates are in an unusual format:

    X         Y 

 1411246  864457.9
 1266383 1231047.3
 1213326  905663.0
 1307353 1213580.7
 1168558  814300.4

I am not sure why this happens. In other (not sure if related) posts some suggested it may be related to the proj4.

How do I transform the resultant coordiantes into "standard" longitude latitude coordinates?

Answer 1

With sf (more modern R spatial toolkit):

library(sf)

tx <- st_read("~/Downloads/Texas_Shapefile/Texas_VTD.shp", stringsAsFactors=FALSE)
## Reading layer `Texas_VTD' from data source `/Users/bob/Downloads/Texas_Shapefile/Texas_VTD.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8400 features and 21 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 372991.6 ymin: 412835.3 xmax: 1618133 ymax: 1594958
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=34.91666666666666 +lat_2=27.41666666666667 +lat_0=31.16666666666667 +lon_0=-100 +x_0=1000000 +y_0=1000000 +datum=NAD83 +units=m +no_defs

tx_ll <- st_transform(tx, "+proj=longlat +ellps=WGS84 +datum=WGS84")

head(st_coordinates(tx_ll))
##              X        Y L1 L2
## [1,] -97.72317 30.40926  1  1
## [2,] -97.72264 30.40920  1  1
## [3,] -97.72255 30.40920  1  1
## [4,] -97.72202 30.40927  1  1
## [5,] -97.72187 30.40931  1  1
## [6,] -97.72156 30.40941  1  1

With legacy tooling:

library(sp)
library(rgdal)

tx <- readOGR("~/Downloads/Texas_Shapefile/Texas_VTD.shp", stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/bob/Downloads/Texas_Shapefile/Texas_VTD.shp", layer: "Texas_VTD"
## with 8400 features
## It has 21 fields
## Integer64 fields read as strings:  CNTY COLOR VTDKEY CNTYKEY Gov_D_02 Gov_R_02 Pres_D_04 Pres_R_04 Gov_D_06 Gov_R_06 Pres_D_08 Pres_R_08 Gov_D_10 Gov_R_10 vap 

proj4string(tx)
## [1] "+proj=lcc +lat_1=34.91666666666666 +lat_2=27.41666666666667 +lat_0=31.16666666666667 +lon_0=-100 +x_0=1000000 +y_0=1000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

tx_ll <- spTransform(tx, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

head(coordinates(tx_ll))
##        [,1]     [,2]
## 0 -97.71355 30.40027
## 1 -95.73430 29.86986
## 2 -97.13773 33.22311
## 3 -97.77786 30.29460
## 4 -96.70334 33.05514
## 5 -98.25908 29.47624

Answer 2

Geospatial data is often stored using coordinates that are not lat-long, because these are spherical coordinates and geometry on a sphere is complex. If you are working in a small area and don't care about the inescapable distortion of pretending the earth is flat, then projected coordinates are used. You can see the projection description in "Proj.4" format with:

> proj4string(shp)
[1] "+proj=lcc +lat_1=34.91666666666666 +lat_2=27.41666666666667 +lat_0=31.16666666666667 +lon_0=-100 +x_0=1000000 +y_0=1000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

which is stored in the .prj file of the set of files that make up a shapefile.

The important one here is the "+proj=lcc" part, which tells us this is a "Lambert Conformal Conic" projection - take a cone, put it above the earth, project the sphere onto the cone, and unwrap the cone.

To convert from one projection to another you need the Proj.4 code for the projection you are targeting. lat-long as reported by GPS units is also known as WGS84, and is the most common lat-long system (there are others based on different shaped earths and different definitions of 0 degrees). It has a Proj.4 string of "+init=epsg:4326", which is a shortcut for several parameters like your shp object has. You can transform sp objects using spTransform:

 > shpLL = spTransform(shp, "+init=epsg:4326")

You may also want to look at the sf package which provides faster functions for reading and transforming shapefiles.