# Maintaining coordinates when converting SpatialPolygonsDataFrame to simple feature in R

I am ultimately trying to determine the area of a home range, but I need to remove areas that overlap land. To accomplish this, I have calculated and plotted my home range, and want to use one polygon (land) to clip the other (home range), leaving only the area over water.

For the clipping code to work, both the HR polygon and the land polygon need to be in the same format (sf) and projection. After converting my HR polygon to sf format, the coordinates change (by a lot) and no longer overlap my land area. I've put together an example below.

``````library(sp)
library(sf)
library(ggmap)
library(adehabitatHR)
library(rnaturalearth)
library(lwgeom)

data<-data.frame(x=c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),y=c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
data\$id<-as.factor("a")

#create a SpatialPointsDataFrame by defining coordinates
coordinates(data)<-c("x","y")

#set CRS
proj4string(data)<-CRS("+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

# run code to determine minimum convex polygon
mcp<-mcp(data, percent = 95, unin = c("m"),
unout = c("km2"))

#view resulting polygons by plotting them
plot(data, col=as.factor(data@data\$id), pch=16)
plot(mcp, col=alpha(1:7, 0.5),add=TRUE)

#now plot polygons on map
#import map
canada <- ne_countries(country="canada",scale = "large", returnclass = "sf")

ggplot(data = canada) +
geom_sf(data=canada, fill="gray80") +
coord_sf(xlim=c(-58,-47), ylim=c(46, 52)) + # map coordinates
labs(x = "Longitude", y = "Latitude", colour="") + #labels
geom_polygon(data=mcp, fill="goldenrod",alpha=0.7, aes(x=long, y=lat, group=group))
``````

``````# since these MCPs overlap land a lot, I want to clip away the parts that are on land. I will do this by clipping the MCP polygon with the map polygon.

#first, see what class of data the HR polygon and map polygon are
class(canada) # sf, so this is good.
class(mcp) # sp, so needs to be converted to sf

# coerce sp object to sf
mcp.poly <- st_as_sf(mcp,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",agr="constant",coords=c("x","y"))

class(mcp.poly) # now this is sf too

#next check that the CRS matches for both the HR polygon and land
#check utm format using:
st_crs(canada)
st_crs(mcp.poly)

#Since they differ, convert both to match
nl<-st_transform_proj(canada,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
mcp.poly<-st_transform_proj(mcp.poly,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

#plot these to check that they still make sense before clipping
``````

``````ggplot() +
geom_sf(data = nl, colour = "light gray", fill = "gray80")+
geom_sf(data = mcp.poly,aes(group=id),fill = "gray")

#Coordinates are wrong
# map plots correctly alone:
ggplot() +
geom_sf(data = nl, colour = "light gray", fill = "gray80")

# but the coordinates on the polygon have changed. Why?
ggplot() +
geom_sf(data = mcp.poly,aes(group=id),fill = "gray")

## from here, I want to use the below code to clip the area of the polygon over land, and then determine the area of the polygon that covers water in km2

# difference between world polygons and the mcp
difference <- st_difference(mcp.poly, nl)

# coerce back to sp
difference <- as(difference, 'Spatial')

#determine area
poly.area<-area(difference)
``````

## 2 Answers

This data looks like lat-long coordinates:

``````data<-data.frame(x=c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),y=c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
``````

but you set its CRS to a UTM 21 region:

``````#set CRS
proj4string(data)<-CRS("+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
``````

which is a coordinate system in metres. Assigning to `proj4string` in `sp` doesn't transform coordinates, so now you have data in UTM 21 metres.

You need to assign this data with its correct CRS string, which is probably something like lat-long in the WGS84 system, ie `"+init=epsg:4326"`.

If you want to transform `sp` data, use `spTransform`.

Also, when you convert to `sf`, don't force a CRS unless you know the CRS from the source is wrong. ie don't do this:

``````# coerce sp object to sf
mcp.poly <- st_as_sf(mcp,"+proj=utm +zone=21 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",agr="constant",coords=c("x","y"))
``````

The reason the `ggplot` works is because you overlay the `mcp` using `geom_polygon(data=mcp,...`, which isn't aware of the coordinate system and so plots using the underlying numbers, which are the right lat-long numbers. But once `geom_sf` enters the picture and thinks "these numbers are UTM 21, lets convert them to lat-long" gets it all wrong for you.

Thanks @Spacedman for you help -- I managed to accomplish what I was after. The data is indeed long-lat coordinates, but I wanted to use UTM as my understanding is that this is more accurate when dealing with distances and areas.

After setting the CRS of the raw data to long-lat, I converted it to UTM prior to running my home range code and converting the polygon to sf.

Here is the complete solution.

``````data<-data.frame(x=c(-50.3, -49.9, -50.0, -50.6, -55.3, -55.4, -55.5, -55.3, -54.9, -54.4, -51.5, -51.2, -50.8, -50.3),y=c(50.3, 48.8, 48.1, 47.4, 48.2, 47.4, 50.1, 48.1, 47.5, 50.7,50.4, 50.7, 50.5, 48.3))
data\$id<-as.factor("a")

#create a SpatialPointsDataFrame by defining coordinates
coordinates(data)<-c("x","y")

# currently there is no CRS for this data. Since it is lat/long, we will set it as such:
proj4string(data)=CRS("+proj=longlat +datum=WGS84 +no_defs")

# now it needs to be converted to UTM
data<- spTransform(data, CRS("+init=epsg:32621 +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 "))

# run code to determine minimum convex polygon
mcp<-mcp(data, percent = 95, unin = c("m"),
unout = c("km2"))

#view resulting polygons by plotting them
plot(data, col=as.factor(data@data\$id), pch=16, axes=TRUE)
plot(mcp, col=alpha(1:7, 0.5),add=TRUE)

#now plot polygons on map
#import map and check projection
canada <- ne_countries(country="canada",scale = "large", returnclass = "sf") # in rnaturalearth package
st_crs(canada)

#convert to utm
nl<-st_transform_proj(canada,"+init=epsg:32621 +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

#now plot to check
ggplot(data = nl) +
geom_sf(data=nl, fill="gray80") +
coord_sf(xlim=c(100000,1200000), ylim=c(5000000, 6100000)) + # map coordinates in utm
labs(x = "Longitude", y = "Latitude", colour="") + #labels
geom_polygon(data=mcp, fill="goldenrod",alpha=0.7, aes(x=long, y=lat, group=group))

# since these MCPs overlap land a lot, I want to clip away the parts that are on land. I will do this by clipping the MCP polygon with the map polygon.

#check what class of data the HR polygon and map polygon are (need to be sf for clipping)
class(nl) # sf, so this is good.
class(mcp) # sp, so needs to be converted to sf

# coerce sp object to sf
mcp.poly <- st_as_sf(mcp)

class(mcp.poly) # now this is sf too

#plot these to check that they still make sense before clipping
ggplot() +
geom_sf(data = nl, colour = "light gray", fill = "gray80")+
geom_sf(data = mcp.poly,aes(group=id),fill = "gray")
``````

``````# now clip HR polygon with land
difference <- st_difference(mcp.poly, st_union(nl))

#determine area (while in sf format)
st_area(difference)
# 106822885279 [m^2]
#  = 106822.89 [km^2]

# plot the result (convert back to sp)
difference <- as(difference, 'Spatial')
plot(difference, axes=TRUE)
``````