Skip to content

Instantly share code, notes, and snippets.

@naomitress
Last active August 26, 2022 07:54
Show Gist options
  • Select an option

  • Save naomitress/ea9d80c7683a57606d720e8dcc5d8e25 to your computer and use it in GitHub Desktop.

Select an option

Save naomitress/ea9d80c7683a57606d720e8dcc5d8e25 to your computer and use it in GitHub Desktop.
OTN receivers from geoserver layer
#making a CSV of OTN receivers, filtering for date, and locations
library(tidyverse)
library(lubridate)
#Set project bounds ####
proj_start <- ymd("20190101")
proj_end <- ymd("20200101")
proj_long_upp <- -40.00
proj_long_low <- -70.00
proj_lat_upp <- 60.00
proj_lat_low <- 40.00
#import it once so we can do different things
geoserver_receivers <-readr::read_csv('https://members.oceantrack.org/geoserver/otn/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=otn:stations_receivers&outputFormat=csv', guess_max = 13579)
#use all the bounds to filter our data set
otn_stations <- geoserver_receivers %>%
filter(!is.na(deploy_date)) %>% #remove if deploy_date is null
filter((deploy_date > proj_start & deploy_date < proj_end)|
(recovery_date < proj_end & recovery_date > proj_start) | #keep regardless of recover_date is null, since it might still be in the water
(deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(18,'months')) | #vr2w / default
(grepl('VR3', model) & deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(4,'years'))| # duration related to model in order to capture
(grepl('VR4', model) & deploy_date < proj_end & is.na(recovery_date) & deploy_date > proj_start - duration(6,'years')) # records which may not have been fully reported
)%>%
filter(stn_lat >= proj_lat_low & stn_lat <= proj_lat_upp & #filter for spatial bounds
stn_long >= proj_long_low & stn_long <= proj_long_upp)
#maybe you want to export it...
#write_csv(
# otn_stations,'otn_stations_geoserver_filtered.csv' )
@bowlbyh
Copy link

bowlbyh commented Feb 17, 2022

Hi Naomi,
This is so easy to use it is almost scary. I was wondering about incorporating a couple of things, if they were interesting to you. I really like the filter above by dates. I was wondering if we could make the filter by positions a bit more complex by bounding it using the Canadian EEZ. I have a function that creates the Canadian/US EEZ and accounts for the keyhole for St .Pierre and Miqulon, and puts the 200m bathymetric contour on a map. Forgive me that I just copy and pasted below. This function relies heavily on the 'sf' library. Once the EEZ is defined, then we can just intersect (command sf::st_intersect(eez.CAN.clip, otn_stations) to specifically return only stations in Canadian waters. I think that these will only be in the marine environment, because they would be bounded by the land raster.

create.eez<-function()
{
#Get the Canadian EEZ filter and plotting sf objects
tem <- tempfile()
download.file("https://raw.githubusercontent.com/Mar-scal/GIS_layers/master/EEZ/EEZ.zip",tem)
tem2 <- tempfile()
unzip(zipfile = tem,exdir = tem2)
eez.CAN <- sf::st_read(paste0(tem2,"/EEZ.shp"))
rm(tem,tem2)
#Need to take one version for the intersection filter and one for pretty maps#
#Pretty maps:
eez.can.plot <- eez.CAN
eez.can.plot <- sf::st_transform(eez.CAN, 2960)
sf::st_crs(eez.can.plot)
#EEZ intersection filter:
eez.CAN <- concaveman::concaveman(eez.CAN)
eez.CAN <- sf::st_transform(eez.CAN, 2960)

#USA EEZ *Note: The site which provides the shape file used here has the download behind a sign-up and can't be pulled
#straight via download.file. You can find the zip for this at https://www.marineregions.org/downloads.php#eez
#under "Marine and land zones: the union of world country boundaries and EEZ's"
eez.us <- sf::st_read("C:/mydocs/whiteCAN/Data/EEZ_Land_v3_202030.shp")
eez.us <- eez.us %>% filter(UNION == "United States")
eez.us <- sf::st_transform(eez.us, 2960)

#Create a polygon to clip the Canadian and US EEZ areas to the Atlantic Coast#

cliped.eez <- sf::st_as_sf(data.frame(X = c(-90,-90,-50,-50),Y = c(23,49,49,23)), coords = c("X","Y"),crs = 4326)
cliped.eez <- sf::st_cast(sf::st_combine(cliped.eez),"POLYGON")
cliped.eez <- sf::st_transform(cliped.eez, 2960)

#Clip the EEZ files#
eez.CAN.clip <<- sf::st_intersection(eez.CAN,cliped.eez)
eez.US.clip <<- sf::st_intersection(eez.us,cliped.eez)

#EEZ's clipped. Now get the map data in and clipped to the same area#
#Canada
poly_ca <- raster::getData('GADM', country='CA', level=1)
sf_ca <- sf::st_as_sf(poly_ca)
atl_Can <<- sf::st_crop(sf_ca, sf::st_as_sfc("POLYGON((-90 23,-90 49,-50 49,-50 23))"))
#USA
poly_us <- raster::getData('GADM', country='USA', level=1)
sf_us <- sf::st_as_sf(poly_us)
atl_us <<- sf::st_crop(sf_us, sf::st_as_sfc("POLYGON((-90 23,-90 49,-50 49,-50 23))"))
#Saint Pierre and Miquelon
poly_SPM <- raster::getData('GADM', country='SPM', level=1)
sf_SPM <- sf::st_as_sf(poly_SPM)
atl_SPM <<- sf::st_crop(sf_SPM, sf::st_as_sfc("POLYGON((-90 23,-90 49,-50 49,-50 23))"))

#Adding a bathymetry element to the maps#

Download specific layers to the current directory

bathy <- sdmpredictors::load_layers("MS_bathy_5m",datadir = "C:/mydocs/whiteCAN/Data/")

https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1466-8238.2011.00656.x (citation for Bio-ORACLE data)

#Clip bathy to area in question
bathy.clip <- sf::st_as_sf(data.frame(X = c(-90,-90,-50,-50),Y = c(23,49,49,23)), coords = c("X","Y"),crs = 4326)
NATL.bathy <- raster::crop(bathy,bathy.clip)
#Get bathy into same projection as track data
NATL.bathy.nad83<-raster::projectRaster(NATL.bathy,crs=2960)
con <- raster::rasterToContour(NATL.bathy.nad83, levels = -200)
canus.bath <<- sf::st_as_sf(con)

#For only the Canadian EEZ#
#land map:
Can_only <- sf::st_crop(sf::st_make_valid(atl_Can), sf::st_as_sfc("POLYGON((-68 40,-68 49,-58 49,-58 40))"))
Can_only <- sf::st_transform(Can_only, 2960)
US_bit <- sf::st_crop(sf::st_make_valid(atl_us), sf::st_as_sfc("POLYGON((-68 40,-68 49,-58 49,-58 40))"))
US_bit <- sf::st_transform(US_bit, 2960)
#EEZ:
CanEEZclip <- sf::st_as_sf(data.frame(X = c(-68,-68,-50,-50),Y = c(39.5,49,49,39.5)), coords = c("X","Y"),crs = 4326)
CEC <- sf::st_cast(sf::st_combine(CanEEZclip),"POLYGON")
CEC <- sf::st_transform(CEC, 2960)
OCEC <<- sf::st_intersection(eez.can.plot,CEC)

#Clip bathy to same size#
CBC <- sf::st_intersection(canus.bath,CEC)

Canada.Plot <<- ggplot() +
geom_sf(data = CBC, color = "#85C1E9",lwd = 1)+
geom_sf(data = Can_only, color = "darkgreen", fill = "lightgreen") +
geom_sf(data = US_bit, color = "darkgrey", fill = "grey") +
geom_sf(data = OCEC, fill = NA)

#For when the Whole Eastern seaboard needs to be plotted we'll need a larger scale map#
CANUS <- rbind(atl_Can, atl_us,atl_SPM)
CANUS <<- sf::st_transform(CANUS, 2960)

to get this out of the function - note that you can't use a double assigned object within a function (can't double assign <<- it above)

eez.CAN<<-eez.CAN
return('NULL')
}

run function - warnings are just about spatial assumptions

create.eez()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment