Calculate longitude/latitude for geographic centroid

Question

I want to do some spatial statistic analysis with the county-level crop yield data in Nebraska for the STAT class. For that I need the longitude and latitude of the geographic centroids of each county. Anybody know how to do it in R? I know it can be done in ArcGIS but I have no access to it now.

Rahm

This question should be moved to gis.stackexchange.com. – Permeance

Answer 1

You didn't give any details where you got your shapefile from, but I got one from here and you can use gCentroid from rgeos thusly:

library(rgdal)
library(sp)
library(rgeos)

nebraska <- readOGR("CountyBoundsUTM/", "CountyUTM")

gCentroid(nebraska, byid=TRUE)

## SpatialPoints:
##           x       y
## 0  721768.5 4636738
## 1  430938.8 4524651
## 2  698036.4 4566570
## 3  370970.6 4641340
## ...
## 89 623301.6 4603228
## 90 618883.0 4486931
## 91 439295.3 4582756
## 92 493680.8 4522680
## Coordinate Reference System (CRS) arguments: +proj=utm +zone=14 +datum=NAD83
## +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 

Barnabe

Hi @hrbrmstr. I've tried your solution to get the result as a data.frame but then I loose the id information. Like this df <-as.data.frame(gCentroid(gridriopoly500, byid=TRUE)) . Is it possible to keep the id info? – Maemaeander

You might want to post that as a separate question. There's a way to do it (easily) but I'd rather not try to find that CountyUTM file to show it (and the OP here didn't need it, apparently). – Barnabe

Thank you @hrbrmstr. I've just posted a new question, here – Maemaeander

Answer 2

You can also extract centroids of SpatialPolygons* objects with coordinates, though the centroids won't be returned as SpatialPoints as with rgeos::gCentroid.

For example:

library(rgdal)
download.file('http://dds.cr.usgs.gov/pub/data/nationalatlas/countyp020_nt00009.tar.gz', 
              f <- tempfile()) # ~ 4.5 Mb
untar(f, exdir=tempdir())
counties <- readOGR(tempdir(), 'countyp020')
xy <- coordinates(counties)
head(xy)

#        [,1]     [,2]
# 0 -153.3905 69.30193
# 1 -156.0582 71.33094
# 2 -155.6695 71.24763
# 3 -155.5164 71.23148
# 4 -155.1846 71.18189
# 5 -155.6126 71.00725

Note that, as pointed out by @Spacedman in the comments, the polygons should be projected to a planar coordinate system first.

Lieu

Those coordinates are computed using lat-long as planar geometry which is not the right way to go about it. Convert to a planar CRS, get centroids, turn into SpatialPoints, convert back to Lat-long. In this case the difference is pretty small but its something you should know about. – Marabout

Answer 3

Recent versions of the sf package (I think since version 1) use the S2 library from Google for spherical geometry calculations. The advantage is that centroid calculations are not simply planar. The relevant method is st_centroid(). An example for a region with a significant spatial extent:

library(rnaturalearth)
# ne_countries() returns 'sp'-type data by default
nc <- ne_countries(continent = "Asia", returnclass = "sf")
  
library(sf)
# long-lat data in WGS84
st_crs(nc)
# use st_geometry() to plot only the polygons and not the associated data
plot(st_geometry(nc), axes = T)
plot(st_centroid(st_geometry(nc)), pch = "+", col = "red", add = T)

# 'sf' integrates nicely with 'ggplot2':
library(ggplot2)
ggplot(nc) + geom_sf() +
    geom_sf(aes(geometry = st_centroid(st_geometry(nc))), colour = "red")

Gurnard

Answer 4

You can use the get_map() function from the ggplot2 package to extract the US county map data from the maps package to a dataframe. Then you can calculate the mid points of the ranges of the lat/lon columns by county (or whatever method you want to use to define geographic center).

Winner

mean of lat/lon ranges are only a rough close approximation to the centroid of the polygons. – Barnabe