9

I'm trying to find some centroids of a shapefile for Danish municipalities and afterwards find the driving time between them. I use R's readShapePoly function from maptools combined with the gCentroidfunction from rgeos, and everything works. However, I get spatialpoints such as

SpatialPoints:
     x       y
1 571860.7 6225016
Coordinate Reference System (CRS) arguments: NA 

Which is clearly not something I can use in Google to grab travel times. I'm looking for a way to convert these numbers to longitude-latitude, but have no idea how.

When I read the data using readOGR from the rgdal library I get the same coordinates but it tells me the following about what I assume is the projection (but the coordinates are the same)

Slot "proj4string":
CRS arguments:
+proj=utm +zone=32 +ellps=intl +units=m +no_defs

Reproducible example: I've put the data for the example here: https://github.com/sebastianbarfort/shapefiles

This should reproduce the problem:

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

map = readShapePoly("~/Downloads/shapefiles-master/kommuner1983.shp")
centroid = gCentroid(map)
centroid
14

Use spTransform to transform the coordinates to WGS84:

library("rgdal")
library("rgeos")

map <- readOGR(".", "kommuner1983")
map_wgs84 <- spTransform(map, CRS("+proj=longlat +datum=WGS84"))
plot(map_wgs84, axes=TRUE)

plot

gCentroid(map_wgs84)
# SpatialPoints:
#       x     y
# 1 10.05 55.96
# Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
# +ellps=WGS84 +towgs84=0,0,0 

rgdal::readOGR is capable to read the projection information automatically. maptools functions neither read nor write projection information, leaving it up to you to manage these details manually.

  • This was just what I was looking for. – sBarfort Dec 3 '14 at 16:35
  • @sBarfort If the answer solves your question then you might want to accept it. That way others will know that it is the correct answer. – R.K. Dec 4 '14 at 1:03
  • You do not have the same control as gCentroid, but SpatialPolygons hold the centroid in the object. You can retrieve them as a matrix using: coordinates(map_wgs84) – Jeffrey Evans Dec 4 '14 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.