Convert SpatialPolygonsDataFrame to projected coordinates using spTransform

Viewed 1268

Im trying to do a point pattern analysis. To do this I have to convert a SpatialPolygonsDataFrame so it contains projected coordinates instead of curved coordinates. However I keep getting the same error:

Error in as.owin.SpatialPolygons(Netherlands_flat) : Only projected coordinates may be converted to spatstat class objects

this is the data I used for a border:

download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip",destfile="ne_10m_admin_1_states_provinces.zip")
unzip("ne_10m_admin_1_states_provinces.zip",exdir="NaturalEarth")
border <- shapefile("NaturalEarth/ne_10m_admin_1_states_provinces.shp")

#extract the border of the Netherlands 
Netherlands <- border[paste(border$iso_a2)=="NL",]

Im able to plot the plot the Netherlands with the events.

#Plot 
plot(babesia$Longitude, babesia$Latitude, pch="+",cex=0.5, xlim=c(3.360782, 7.227095), ylim = c(50.723492, 53.554584))
plot(Netherlands, add = T)

Netherlands with events

But upon using the Spatstat package I keep running into this error.

I tried this code to transform the coordinates

coord_netherlands <- coordinates(Netherlands)
proj4string(Netherlands)
summary(Netherlands)

Netherlands_flat <- spTransform(coord_netherlands, CRS("+proj=longlat +datum=WGS84 +no_defs"))
Netherlands <- as.owin(Netherlands_flat)

Error in as.owin.SpatialPolygons(Netherlands_flat) : Only projected coordinates may be converted to spatstat class objects

Does anyone know how to solve this? Thank you very much in advance!

1 Answers

You are almost there. You just need to project to another coordinate system when you call spTransform. You currently ask for geographic coordinates on a spheriod model of the earth (long,lat). Instead you should ask for a flat (x,y) coordinate system. This could be utm coordinates in the appropriate zone for the Netherlands or there might well be a better alternative. Your events also need to be transformed from (long,lat) to the same coordinate system. Maybe you can look at the shapefile vignette of the spatstat package for an example. Or look under the spatstat tag on this site. I'm on my phone do I can't give detailed help. Good luck.

If your events are in a data.frame called xy you can project to UTM zone 31N like this:

xy <- data.frame(lon = 1:2, lat = 1:2)
coordinates(xy) <- ~lon+lat
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")
xy
# SpatialPoints:
#      lon lat
# [1,]   1   1
# [2,]   2   2
# Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
# +ellps=WGS84 +towgs84=0,0,0 
events.utm <- spTransform(xy, CRS("+proj=utm +zone=31N +datum=WGS84"))
events.utm
# SpatialPoints:
#           lon      lat
# [1,] 277438.3 110598.0
# [2,] 388786.7 221094.9
# Coordinate Reference System (CRS) arguments: +proj=utm +zone=31N
# +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
Related