The Marine-Analyst provides augmented data access and reproducible data analysis for marine data. The Marine Analyst is intended for the general public with a focus on educational content. It allows students, teachers and academia to access marine knowledge. Marine topics are developed under the web portal main menu. The Marine Analyst is an essential tool for oceanographers, scientists, engineers, managers and policy-makers who are analysing the state and dynamics of Europe's seas. The Marine Analyst Jupyter hub provides you with usefull data analysis codes. Alternatively you can manage your analyses and download links for your data selections via a personal dashboard. This Jupyter notebook is licensed under the MIT License.
# Edit the longitude and latitude coordinates to define the geographical area:
###################################
minlon=11.3 #minimum longitude
minlat=53.6 #minimum latitude
maxlon=15.5 #maximum longitude
maxlat=55.9 #maximum latitude
###################################
# This Jupyter Notebook has been automatically generated from R Markdown Notebook and may contain inconsistencies.
# Edit it, execute it and save it as HTML document.
# Anytime you have a question, a suggestion or an issue, you can contact us at my-beach@knowcean.eu.
# Your collaboration is highly appreciated!
###################################
wdpaid=paste(minlon,minlat,maxlon,maxlat,sep="_")
wdpaidsplit <- unlist(strsplit(wdpaid, "[_]"))
xmin <- as.numeric(wdpaidsplit[1])
ymin <- as.numeric(wdpaidsplit[2])
xmax <- as.numeric(wdpaidsplit[3])
ymax <- as.numeric(wdpaidsplit[4])
Sessionid<-"JupyterNotebook752"
source_provider <- "EMODnet Chemistry"
source_provider_url <- "https://www.emodnet.eu"
layer_title<-"Beach Litter - Mean number of Cigarette"
layer="bl_cigarette_cleaning"
wfs_url <- "https://www.ifremer.fr/services/wfs/emodnet_chemistry2?"
wms_url <- "https://www.ifremer.fr/services/wms/emodnet_chemistry2?"
wms_layer="bl_cigarette_cleaning"
layer_id<-752
map_legend <- "litterabundance"
map_label<-"beachname"
link_csv<-paste0("./Report-", layer_id, "_", Sessionid, "_", wdpaid, "-csvfile.csv",sep="")
csvfile_name = paste("Report-", layer_id, "_", Sessionid, "_", wdpaid, "-csvfile.csv",sep="")
link_geojson<-paste0("./Report-", layer_id, "_", Sessionid, "_", wdpaid, "-geojsonfile.geojson",sep="")
geojsonfile_name = paste("Report-", layer_id, "_", Sessionid, "_", wdpaid, "-geojsonfile.geojson",sep="")
temp_path<- "."
library(rgdal)
library(downloader)
library(ggplot2)
library(mapdata)
library(ggmap)
library(ggrepel)
library(httr)
library(sf)
library(rasterVis)
library(rgeos)
library(sp)
library(raster)
library(dplyr)
library(XML)
# Script for Wekeo environment
sr=SpatialPolygons(list(Polygons(list(Polygon(cbind(c(xmin, xmin, xmax, xmax),c(ymax, ymin, ymin, ymax)))),"1")))
mpa=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:1), row.names=c("1")))
proj4string(mpa)<-CRS("+proj=longlat +datum=WGS84")
bbox<-paste(xmin,ymin,xmax,ymax,sep=",")
# Link to Marine Analyst dataset page
link_marineanalyst <- paste0("http://marine-analyst.eu/dev.py?N=simple&O=",layer_id,"&maxlat=",ymax,"&maxlon=",xmax,"&minlon=",xmin,"&minlat=",ymin)
# Link to open the openlayer page for EMODnet HA
openlayer<-paste0("http://www.marine-analyst.eu/openlayers3/openlayer.py?wms_url=",wms_url,"/wms&wms_layer=",layer,"&bbox=",bbox)
print (paste("West-Longitude:",round(xmin,2)))
print (paste("South-Latitude:",round(ymin,2)))
print (paste("East-Longitude:",round(xmax,2)))
print (paste("North-Latitude:",round(ymax,2)))
value<-(xmax-xmin)*(ymax-ymin)
if (value > 100) {
zoom_value<-6
} else if (value > 1) {
zoom_value<-7
} else {
zoom_value<-8
}
base<-get_map(location=c(xmin-1,ymin-1,xmax+1,ymax+1), zoom=zoom_value, maptype="terrain-background", source = "stamen")
terrain <- ggmap(base)
map <- terrain + geom_polygon(data=mpa,aes(x=long,y=lat,group=group,fill="mpa"),colour="green",fill="blue",alpha=.1) +
ggtitle("")+xlab("Longitude")+ylab("Latitude")
plot(map)
r toString(layer_title)
¶# DescribeFeatureType request function
DescribeFeatureType<-function(layer){
layer<-as.character(layer)
con<-paste0(wfs_url,"service=WFS&version=1.1.0&request=DescribeFeatureType&typeName=",layer,"&outpuformat=XMLSCHEMA")
xml <- "file.xml"
xml <- tempfile(xml)
httr::GET(con,write_disk(xml))
xmldoc <- XML::xmlParse(xml)
xml_data <- XML::xmlToList(xmldoc)
data.catalog <- data.frame(t(xml_data$complexType$complexContent$extension$sequence),row.names=NULL)
return(data.catalog)
}
WFS_DescribeFeatureType <- DescribeFeatureType(layer)
WFS_Colnames<-c()
for (i in 1:ncol(WFS_DescribeFeatureType)) {
WFS_Colnames<-append(WFS_Colnames, WFS_DescribeFeatureType[i]$element$element[1][["name"]])
}
WFS_Colnames
rez_nblist<- c(1:length(WFS_Colnames))
getWFSgml3<-function(layer){
layer<-as.character(layer)
con<-paste0(wfs_url,"service=WFS&version=1.0.0&request=GetFeature&typeName=",layer,"&OUTPUTFORMAT=gml3&srsName=EPSG%3A4326")
pipo<-sf::st_read(con)
return(pipo)
}
wfs_data<-getWFSgml3(layer)
#Transform mpa in Simple feature collection to perform the subsetting because wfs_data contains the whole info - use bbox and tyname are exclusive (Mapserver)
mpaSP <- as(mpa, "SpatialPolygonsDataFrame")
wfs_data<-wfs_data[st_as_sf(mpaSP),]
# get the geometry as lat and lon cols
wfs_data <- wfs_data %>% dplyr::mutate(lat = sf::st_coordinates(.)[,1],lon = sf::st_coordinates(.)[,2])
st_write(wfs_data, file.path(temp_path,csvfile_name), layer = csvfile_name, driver = "csv", delete_dsn = TRUE)
st_write(wfs_data, file.path(temp_path,geojsonfile_name), layer = geojsonfile_name, driver = "GeoJSON", delete_dsn = TRUE)
if(nrow(wfs_data) > 0) {
wfs_data
} else {
print("No data available for the defined geographical extent")
}
if(nrow(wfs_data) > 0) {
map <- ggplot() +
borders("worldHires", fill = "gray", colour = "black", xlim = range(xmin,xmax), ylim = range(ymin,ymax), size = .25) +
theme(legend.position = "bottom") +
theme(panel.grid.minor.y= element_blank(), panel.grid.minor.x = element_blank()) +
geom_polygon(data=mpa,aes(x=long,y=lat,group=group,fill="mpa"),colour="green",fill="blue",alpha=.1) +
geom_sf() +
geom_point(data = wfs_data, aes(x = lat, y = lon, size =.data[[map_legend]]), fill = "red", color = "red", alpha = .4) +
coord_sf(xlim = c(xmin, xmax),ylim = c(ymin, ymax))+
ggtitle(layer_title)+xlab("Longitude (x)")+ylab("Latitude (y)")
map
} else {
print("No data available for the defined geographical extent")
}
if(nrow(wfs_data) > 0) {
centroid<- st_centroid(wfs_data)
centroid<- cbind(wfs_data, st_coordinates(st_centroid(wfs_data$geometry)))
map <- ggplot() +
borders("worldHires", fill = "gray", colour = "black", xlim = range(xmin,xmax), ylim = range(ymin,ymax), size = .25) +
theme(legend.position = "bottom") +
theme(panel.grid.minor.y= element_blank(), panel.grid.minor.x = element_blank()) +
geom_polygon(data=mpa,aes(x=long,y=lat,group=group,fill="mpa"),colour="green",fill="blue",alpha=.1) +
geom_sf() +
geom_point(data = wfs_data, aes(x = lat, y = lon, size = .data[[map_legend]]), fill = "red", color = "red", alpha = .4) +
geom_text(data=centroid,aes(x=lat, y=lon, label=.data[[map_label]]), color = "black", fontface = "bold", size=2, hjust= 0, vjust=2, check_overlap = TRUE) +
coord_sf(xlim = c(xmin, xmax),ylim = c(ymin, ymax))+
ggtitle(layer_title)+xlab("Longitude (x)")+ylab("Latitude (y)")
map
} else {
print("No data available for the defined geographical extent")
}
if(nrow(wfs_data) > 0) {
getWMSmap<-function (wms_layer,xmin,xmax,ymin,ymax)
{
width <- 960
height <- as.integer(width * (ymax-ymin) / (xmax-xmin))
wms_layer<-as.character(wms_layer)
bbox <- paste(xmin, ymin, xmax, ymax, sep = ",")
con<-paste0(wms_url,"/wms?SERVICE=WMS&VERSION=1.1.0&request=GetMap&layers=",wms_layer,"&format=image/jpeg&srs=EPSG:4326&bbox=",bbox,"&height=",height,"&width=",width,"")
wms <- "img.png"
wms <- tempfile(wms)
download(con, wms, quiet = TRUE, mode = "wb")
img <- brick(wms)
names(img) <- c("img.1", "img.2", "img.3")
img[img$img.1 == 255 & img$img.2 == 255 & img$img.3 == 255] <- NA
wms_basemap_url="http://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv"
wms_basemap_layer="gebco_latest"
con<-paste0(wms_basemap_url,"?SERVICE=WMS&VERSION=1.1.0&request=GetMap&layers=",wms_basemap_layer,"&format=image/png&srs=EPSG:4326&bbox=",bbox,"&height=",height,"&width=",width,"")
wms <- "img.png"
wms <- tempfile(wms)
download(con, wms, quiet = TRUE, mode = "wb")
basemap <- brick(wms)
names(basemap) <- c("img.1", "img.2", "img.3")
img <- merge(basemap,img)
img@extent@xmin <- xmin
img@extent@ymin <- ymin
img@extent@xmax <- xmax
img@extent@ymax <- ymax
proj4string(img)<-CRS("+proj=longlat +datum=WGS84")
return(img)
}
wms_img<-getWMSmap(wms_layer,xmin,xmax,ymin,ymax)
rggbplot <- function(inRGBRst,npix=NA,scale = 'lin'){
rgblinstretch <- function(rgbDf){
maxList <- apply(rgbDf,2,max)
minList <- apply(rgbDf,2,min)
temp<-rgbDf
for(i in c(1:3)){
temp[,i] <- (temp[,i]-minList[i])/(maxList[i]-minList[i])
}
return(temp)
}
rgbeqstretch<-function(rgbDf){
temp<-rgbDf
for(i in c(1:3)){
unique <- na.omit(temp[,i])
if (length(unique>0)){
ecdf<-ecdf(unique)
temp[,i] <- apply(temp[,i,drop=FALSE],2,FUN=function(x) ecdf(x))
}
}
return(temp)
}
npix <- ncell(inRGBRst)
x <- sampleRegular(inRGBRst, size=npix, asRaster = TRUE)
dat <- as.data.frame(x, xy=TRUE)
colnames(dat)[3:5]<-c('r','g','b')
if(scale=='lin'){
dat[,3:5]<- rgblinstretch(dat[,3:5])
} else if(scale=='stretch'){
dat[,3:5]<- rgbeqstretch(dat[,3:5])
}
p <- ggplot()+ geom_tile(data=dat, aes(x=x, y=y, fill=rgb(r,g,b))) + scale_fill_identity()
}
}
if(nrow(wfs_data) > 0) {
map <- rggbplot(wms_img)+
#borders("worldHires", fill = "gray", colour = "black", xlim = range(xmin,xmax), ylim = range(ymin,ymax), size = .25) +
coord_quickmap(xlim=range(xmin,xmax),ylim=range(ymin,ymax))+
ggtitle(layer_title)+xlab("Longitude")+ylab("Latitude")
plot(map)
} else {
print("No data available for the defined geographical extent")
}
if(nrow(wfs_data) > 0) {
litterabundance <- "litterabundance"
abundance <- ggplot() +
theme(legend.position = "bottom") +
geom_point(data = wfs_data, aes(x = year, y = litterabundance, color =.data[[map_label]]), size=4) +
ggtitle("Litter abundance per year and location")+xlab("Year")+ylab("Litter abundance")
plot(abundance)
} else {
print("No data available for the defined geographical extent")
}