R Shapefile – Transforming Longitude Latitude Coordinates

coordinateslatitude longitudershapefile

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?

Best Answer

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

Related Solutions

[GIS] locating latitude-longitude point in a piste-map with matrix transforms

If the image is in the least realistic--and it appears to be--then it represents a projection from a portion of Euclidean three-space onto a portion of two-space. We cannot expect this projection to have a nice mathematical formula, because it is likely the artist has introduced irrregular small distortions, but we can expect its formula to somewhat approximate formulas for such projections. These formulas tend to be linear or projective and they also tend to be simpler when Cartesian coordinates are used.

Therefore, what you should consider doing at the outset is this:

Convert all (lat, lon, elevation) values to geocentric (x, y, z) values.
Associated with each (x, y, z) value is a pair of image coordinates (u, v). We probably lose little by estimating u and v independently of each other, which reduces the problem to one of interpolating samples of a scalar function. The previous considerations suggest that you seek interpolating functions that locally are of the form
```
u = a x + b y + c z + d
```
(linear) or even
```
u = (a x + b y + c z + d) / (a' x + b' y + c' z + 1)
```
(projective) and similarly for v. This reduces the problem to one of performing an interpolation in three dimensions, which is sufficiently routine to give you access to a large body of methods and software.

Alternatively, you might consider triangulating the (x,y,z) points and interpolating their associated (u,v) values separately across the triangles; even a linear interpolation across each triangular face could do a pretty good job. This is a 3D analog of what ESRI's 3D Analyst does in two dimensions with TINs.

Please note, though, that unless you have the elevation of a given location, you cannot hope to plot it accurately on this map: changing the elevation will create a nearly vertical line on the map, whence your location could be plotted practically anywhere along this line. With considerable work you could solve this problem for points on the ground by first digitizing the evident discontinuities in the inverse projection--that is, representing the outlines of mountain ridges as polylines--and constructing a local inverse interpolator (respecting these breaklines) that associates to each (u,v) on the image a location (x,y,z) in the world. You could exploit this to identify which point along any vertical line has a (lat,lon) close to the one you started with.

[GIS] Convert Shape Length Area Data to Latitude Longitude

The layer is provided as a polygon so you'll need to convert the polygons to points and calculate the XYs of each point.

Enable 'on the fly' CRS projections in Project Properties
Ensure that the Selected CRS is the CRS you want to project the data to
Add vector layer
Right click layer and "Save vector layer as..."
Select the CRS you want the layer exported to
Save layer and re-add the layer to the project
Convert the polygon to points using Vector-->Geometry Tools-->Extract nodes
Ensure the points layer is added to the project
Open Attributes table for point layer
Toggle editing mode on
Open field calculator
Add two new calculated fields using the $x and $y funcitons under Geometry for each respective field
Save edits and toggle off edit mode

This should give you a point file with the vertices for each polygon and the respective calcualted X and Y values for each point in the chosen CRS.

Best Answer

Related Solutions

[GIS] locating latitude-longitude point in a piste-map with matrix transforms

[GIS] Convert Shape Length Area Data to Latitude Longitude

Related Question