1

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))

enter image description here

# 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

enter image description here

enter image description here

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 2

0

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.

0

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")

enter image description here

# 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)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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