# R shapefile transform longitude latitude coordinates

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?

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)
 "+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.

With `sf` (more modern R spatial toolkit):

``````library(sf)

## 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")

##              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)

## OGR data source with driver: ESRI Shapefile
## 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)
##  "+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"))

• I'm not sure how this explains "why this happens" to the OP. And I wouldn't refer to the `sp` packages as "tooling". – Spacedman Sep 17 '18 at 13:17