Commit e3ec61bb authored by peguerin's avatar peguerin
Browse files

Initial commit

parents
# exclude everything
01-infos/*
*.tif
# exception to the rule
!01-infos/.gitkeep
###############################################################################
## Librairies
library("raster")
library("rgdal")
###############################################################################
## load data
coo=read.csv("01-sample_coords/all_coords.csv",sep="\t")
coo=na.omit(coo)
names(coo)=c("id","lat","lon")
latlon = data.frame(lon=coo$lon,lat=coo$lat)
coordinates(latlon) = c("lon", "lat")
### indiquer que les lat/lon est en projection google map
proj4string(latlon) = CRS("+init=epsg:4326")
## load raster of "chlorophyll benthic"
surface_files=c("01-infos/chlorophyll/surface/Present.Surface.Chlorophyll.Max.tif.zip",
"01-infos/chlorophyll/surface/Present.Surface.Chlorophyll.Mean.tif.zip",
"01-infos/chlorophyll/surface/Present.Surface.Chlorophyll.Min.tif.zip")
benthic_files=c("01-infos/chlorophyll/benthic/Present.Benthic.Mean.Depth.Chlorophyll.Max.tif.zip",
"01-infos/chlorophyll/benthic/Present.Benthic.Mean.Depth.Chlorophyll.Mean.tif.zip",
"01-infos/chlorophyll/benthic/Present.Benthic.Mean.Depth.Chlorophyll.Min.tif.zip")
chloBenthicMeanRa=raster(unzip(benthic_files[1]))
chloBenthicMaxRa=raster(unzip(benthic_files[2]))
chloBenthicMinRa=raster(unzip(benthic_files[3]))
chloSurfaceMeanRa=raster(unzip(surface_files[1]))
chloSurfaceMaxRa=raster(unzip(surface_files[2]))
chloSurfaceMinRa=raster(unzip(surface_files[3]))
##convertir les lat/lon dans la projection du raster "rp"
latlon = spTransform(latlon,proj4string(chloBenthicMeanRa))
###############################################################################
## extract habitat information from coordinates
#radius(rayon) du cercle du buffer ayant pour centre le point d'echantillonnage
#cercle de 2 km de diametre
buffer_distance=1000
chloBenthicMeanExtract=extract(chloBenthicMeanRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chloBenthicMaxExtract=extract(chloBenthicMaxRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chloBenthicMinExtract=extract(chloBenthicMinRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chloSurfaceMeanExtract=extract(chloSurfaceMeanRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chloSurfaceMaxExtract=extract(chloSurfaceMaxRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chloSurfaceMinExtract=extract(chloSurfaceMinRa,latlon,buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
chlotab=cbind(coo,as.numeric(chloBenthicMeanExtract),
as.numeric(chloBenthicMaxExtract),
as.numeric(chloBenthicMinExtract),
as.numeric(chloSurfaceMeanExtract),
as.numeric(chloSurfaceMaxExtract),
as.numeric(chloSurfaceMinExtract))
names(chlotab)=c("id","latitude","longitude",
"chloBenthicMean",
"chloBenthicMax",
"chloBenthicMin",
"chloSurfaceMean",
"chloSurfaceMax",
"chloSurfaceMin")
write.table(chlotab, file="03-results/all_chlorophyll.tsv",quote=F,sep="\t",col.names=T,row.names=F)
## Carte Salinit? Max
library(raster)
library(ggplot2)
library(ggmap)
library(rasterVis)
library(rgdal)
#### Download data
load("data_sal_Max")
data_sal_Max = data.frame(data_sal_Max)
crs.geo <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") # geographical, datum WGS84
coordinates(data_sal_Max) = c("long", "lat")
proj4string(data_sal_Max) <- crs.geo # define projection system of our data
Smax = SpatialPixelsDataFrame(as(data_sal_Max, "SpatialPoints"), data=as(data_sal_Max, "data.frame"), tolerance=0.9)
Smax.df <- as.data.frame(Smax)
gridded(Smax)=TRUE
plot(Smax)
r_Smax=raster(Smax, layer="res", values=T)
s <- stack(r_Smax, r_Smax*2)
#layerNames(s) <- c('meuse', 'meuse x 2')
projection(r_Smax) = CRS(proj4string(Smax))
OR<-readOGR(dsn="/Users/laura/Desktop/", layer="ne_110m_land.shp")
ggplot(Smax.df) +
geom_tile(aes(x=long, y=lat,fill=factor(res),alpha=0.8)) +
geom_polygon(data=OR, aes(x=long, y=lat, group=group),
fill=NA,color="grey50", size=1)
#open ASCII file using ‘raster’ command, which converts the ASCII to a raster object
load("01-infos/data_sal_Max")
map = data.frame(data_sal_Max)
#convert the raster to points for plotting
map.p <- rasterToPoints(map)
#Make the points a dataframe for ggplot
df <- data.frame(map.p)
#Make appropriate column headings
colnames(df) <- c("Longitude", "Latitude")
#Call in point data, in this case a fake transect (csv file with lat and lon coordinates)
sites <- data.frame(read.csv(/your/path/to/pointfile.csv))
#Now make the map
ggplot(data=map, aes(y=lat, x=long)) +
geom_raster(aes(fill=res)) +
geom_point(data=sites, aes(x=x, y=y), color=white, size=3, shape=4) +
theme_bw() +
coord_equal() +
scale_fill_gradient(MAP (mm/yr), limits=c(0,2500)) +
theme(axis.title.x = element_text(size=16),
axis.title.y = element_text(size=16, angle=90),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = right,
legend.key = element_blank()
)
\ No newline at end of file
## Librairies
library("raster")
library("rgdal")
## le nom du fichier geodatabase
fgdb="EUSeaMap2016.gdb"
# afficher le nom des layers dans le fichier geodatabase
fc_list <- ogrListLayers(fgdb)
print(fc_list)
#le nom du layer est "EUSM2016"
# combien de layers dans le fichier geodatabase
attr(,"nlayers")
#il y a 1 layer
## charger le layer EUSM2016 du fichier geodatabase comme objet SpatialPolygon
eusm=readOGR(dsn="EUSeaMap2016.gdb",layer="EUSM2016")
## charger la liste des latitudes/longitudes des echantillons (alicia)
sar.geo=read.table("~/working/seaconnect/environmental/donnees/geolocation_sample_pop_seaconnect2_sar.csv",sep="\t",header=T)
# une dataframe des longitude latitude pour chaque echantillon
latlon = data.frame(cbind(sar.geo$Longitude, sar.geo$Latitude))
colnames(latlon) = c("lon","lat")
coordinates(latlon) = c("lon", "lat")
### indiquer que les lat/lon est en projection google map
proj4string(latlon) = CRS("+init=epsg:4326")
##convertir les lat/lon dans la projection du layer EUSM2016
latlon = spTransform(latlon,proj4string(eusm))
#sur quel polygon se situe chaque point d'echantillonnage
res=intersect(latlon,eusm)
#la valeur substrate sur le polygon ou se situe le point d'echantillonnage
res.substrate=res@data$Substrate
res.substrate=data.frame(as.character(res@data$Substrate), stringsAsFactors=FALSE)
coords=data.frame(res@coords)
#tableau final coordonnées et substrat
final=cbind(res.coords,res.substrate)
#reconvertir lat/lon en projection google map
colnames(final)=c("lon","lat","substrate")
coordinates(final)=c("lon","lat")
proj4string(final)=CRS(proj4string(eusm))
finalg=spTransform(final,CRS("+init=epsg:4326"))
finalg.df=cbind(data.frame(finalg@coords),data.frame(finalg@data$substrate))
colnames(finalg.df)=c("lon","lat","substrate")
### buffer
#correspondance ID/levels(Substrate)
levels.substrate=levels(eusm$Substrate)
#radius(rayon) du cercle du buffer ayant pour centre le point d'echantillonnage
#cercle de 100 km de diametre
buffer_distance=50000
#resolution pixel ~10km
r <- raster(ncol=860,nrow=98)
#resolution pixel ~5km
r <- raster(ncol=1720,nrow=196)
#resolution pixel ~2km
r <- raster(ncol=3440,nrow=392)
extent(r) <- extent(eusm)
#convertir la carte gdb en raster (grille de pixels)
rp <- rasterize(eusm, r, 'Substrate')
#initialiser le dataframe qui contient coordonnes et ratio substrats
substrat.frame <- as.data.frame(matrix(,ncol=length(levels.substrate)+2,nrow=0))
#pour chaque point d'echantillonnage, relever les valeurs des pixels dans le buffer
for(i in seq(1,length(latlon))){
reb=extract(rp,latlon[i],buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
pixs.valeur.occurence=table(reb)
total.pixs=sum(pixs.valeur.occurence)
pixs.valeur.ratio=data.frame(pixs.valeur.occurence/total.pixs)
i.substrats=c()
for(j in seq(1,length(levels.substrate))){
if(j %in% pixs.valeur.ratio$reb){
i.substrats=cbind(i.substrats,pixs.valeur.ratio[which(pixs.valeur.ratio==j),]$Freq)
}else{
i.substrats=cbind(i.substrats,0)
}
}
icoords=data.frame(latlon[i]@coords)
isubstrats=data.frame(i.substrats)
i.frame=cbind(icoords,isubstrats)
names(i.frame) = c("lon","lat",levels.substrate)
substrat.frame=rbind(substrat.frame,i.frame)
}
##ajouter les coordonnees en projection google map
eusm_coords=data.frame(cbind(substrat.frame$lon,substrat.frame$lat))
colnames(eusm_coords)=c("gm_lon","gm_lat")
coordinates(eusm_coords)=c("gm_lon","gm_lat")
proj4string(eusm_coords)=CRS(proj4string(eusm))
gm_coords=spTransform(eusm_coords,CRS("+init=epsg:4326"))
substrat.final=cbind(substrat.frame,data.frame(gm_coords@coords))
#ecrire le tableau avec les coordonnees dans les deux projections et le ratio de chaque substrat
write.table(substrat.final,file="ratio_substrat_coords.csv",sep=";",row.names=F,quote=F)
##############
##### extract environmental variable frome GIS data
##############
setwd("~/Documents/project_SEACONNECT/sub_DSARGUS/data_analyse/environmental_variables/")
## Librairies
library("raster")
library("rgdal")
## le nom du fichier geodatabase
fgdb="EMODnet/R20170615_EUSeaMap2016/EUSeaMap2016.gdb"
# afficher le nom des layers dans le fichier geodatabase
fc_list <- ogrListLayers(fgdb) # read OGR vector maps into spatial objects
print(fc_list)
#le nom du layer est "EUSM2016"
# combien de layers dans le fichier geodatabase
attr(fc_list,"nlayers")
#il y a 1 layer
## charger le layer EUSM2016 du fichier geodatabase comme objet SpatialPolygon
eusm=readOGR(dsn="EMODnet/R20170615_EUSeaMap2016/EUSeaMap2016.gdb",layer="EUSM2016")
str(eusm@data)
## charger la liste des latitudes/longitudes des echantillons (alicia)
seaconnect_geo=read.table("~/Documents/project_SEACONNECT/sub_DSARGUS/data/coord_seaconnect_tous.txt",sep="\t",header=T)
# une dataframe des longitude latitude pour chaque echantillon
latlon = data.frame(cbind(seaconnect_geo$Longitude, seaconnect_geo$Latitude))
colnames(latlon) = c("lon","lat")
coordinates(latlon) = c("lon", "lat")
### indiquer que les lat/lon est en projection google map
proj4string(latlon) = CRS("+init=epsg:4326")
##convertir les lat/lon dans la projection du layer EUSM2016
latlon = spTransform(latlon,proj4string(eusm))
#sur quel polygon se situe chaque point d'echantillonnage
res=intersect(latlon,eusm)
#la valeur substrate sur le polygon ou se situe le point d'echantillonnage
res.substrate=res@data$Substrate
res.substrate=data.frame(as.character(res@data$Substrate), stringsAsFactors=FALSE)
res.coords=data.frame(res@coords)
#tableau final coordonnees et substrat
final=cbind(res.coords,res.substrate)
#reconvertir lat/lon en projection google map
colnames(final)=c("lon","lat","substrate")
coordinates(final)=c("lon","lat")
proj4string(final)=CRS(proj4string(eusm))
finalg=spTransform(final,CRS("+init=epsg:4326"))
finalg.df=cbind(data.frame(finalg@coords),data.frame(finalg@data$substrate))
colnames(finalg.df)=c("lon","lat","substrate")
### get substrate types within buffer
#correspondance ID/levels(Substrate)
levels.substrate=levels(eusm$Substrate)
#radius(rayon) du cercle du buffer ayant pour centre le point d'echantillonnage
#cercle de 100 km de diametre
buffer_distance=50000
## rasterize shapefile for subsequent extraction
# crop shapefile to mediterranean extent
eusm_med <- crop(eusm, extent(-955183, 4151749, 3296265, 5822348))
plot(eusm_med)
#writeOGR(obj=eusm_med, dsn="EMODnet/R20170615_EUSeaMap2016/EUSeaMap2016.gdb",)
# order data by polygon surface area
# https://gis.stackexchange.com/questions/200556/rasterrasterize-set-wrong-values-with-enclosed-polygons
area <- sapply(1:nrow(eusm_med), function(x) {eusm_med@polygons[[x]]@Polygons[[1]]@area})
eusm_med_area <- eusm_med[order(area), ]
# same but without previous cropping
area_t <- sapply(1:nrow(eusm), function(x) {eusm@polygons[[x]]@Polygons[[1]]@area})
eusm_t_area <- eusm[order(area_t), ]
# define a raster with certain resolution & copy shapefile extent
#resolution pixel ~10km: r <- raster(ncol=860,nrow=98)
#resolution pixel ~5km: r <- raster(ncol=1720,nrow=196)
#resolution pixel ~2km: r <- raster(ncol=3440,nrow=392)
#copy the resolution of PAR raster
#r <- raster(ncol=PAR@ncols, nrow=PAR@nrows) #dimensions : 24966, 42246, 1054713636 (nrow, ncol, ncell)
# after copying extent: resolution : 204.9172, 394.3797 (x, y)
# takes a long time! > 48 hrs
#resolution pixel ~1km: r <- raster(ncol=8600,nrow=980)
#resolution pixel ~500m: r <- raster(ncol=17200,nrow=1960)
#resolution pixel ~200m (frome PAR dimensions)
r <- raster(ncol=24966,nrow=42246)
extent(r) <- extent(eusm_med_area)
extent(r) <- extent(eusm_t_area)
# convert shapefile to raster
#give priority to polygons with smallest area (fun = "first")
rp <- rasterize(eusm_med_area, r, 'Substrate', fun = "first")
rp_t <- rasterize(eusm_t_area, r, 'Substrate', fun = "first")
#export raster
writeRaster(rp, "EUSeaMap2016_substrate_raster_med_res200m", format="GTiff")
writeRaster(rp_t, "EUSeaMap2016_substrate_raster_res200m", format="GTiff")
#extract environmental variables
#initialiser le dataframe qui contient coordonnes et ratio substrats
substrat.frame <- as.data.frame(matrix(,ncol=length(levels.substrate)+2,nrow=0))
#pour chaque point d'echantillonnage, relever les valeurs des pixels dans le buffer
for(i in seq(1,length(latlon))){
reb=extract(rp_t,latlon[i],buffer=buffer_distance,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
pixs.valeur.occurence=table(reb)
total.pixs=sum(pixs.valeur.occurence)
pixs.valeur.ratio=data.frame(pixs.valeur.occurence/total.pixs)
i.substrats=c()
for(j in seq(1,length(levels.substrate))){
if(j %in% pixs.valeur.ratio$reb){
i.substrats=cbind(i.substrats,pixs.valeur.ratio[which(pixs.valeur.ratio==j),]$Freq)
}else{
i.substrats=cbind(i.substrats,0)
}
}
icoords=data.frame(latlon[i]@coords)
isubstrats=data.frame(i.substrats)
i.frame=cbind(icoords,isubstrats)
names(i.frame) = c("lon","lat",levels.substrate)
substrat.frame=rbind(substrat.frame,i.frame)
}
##ajouter les coordonnees en projection google map
eusm_coords=data.frame(cbind(substrat.frame$lon,substrat.frame$lat))
colnames(eusm_coords)=c("gm_lon","gm_lat")
coordinates(eusm_coords)=c("gm_lon","gm_lat")
proj4string(eusm_coords)=CRS(proj4string(eusm))
gm_coords=spTransform(eusm_coords,CRS("+init=epsg:4326"))
substrat.final=cbind(substrat.frame,data.frame(gm_coords@coords))
#ecrire le tableau avec les coordonnees dans les deux projections et le ratio de chaque substrat
write.table(substrat.final,file="ratio_substrat_coords.csv",sep=";",row.names=F,quote=F)
#add sampling cell info
sc_substr <- cbind(seaconnect_geo[,c(1:2,6:7)], substrat.final[,3:18])
write.table(sc_substr, file="seaconnect_substrate_res200m.csv",sep=";",row.names=F)
# plot everything
cuts<-seq(rp@data@min-1, rp@data@max, by = 1)
library(colorRamps)
library(maps)
jpeg(file="Substrate_samplingcells_buffer_res200m.jpeg", width = 1700, height = 900, units = "px")
plot(rp_t,
axes = FALSE,
#ext = c(-915183, 4151749, 3296265, 5822348),
ext = c(-955183, 4151749, 3296265, 5822348),
breaks=cuts,
#col=primary.colors(16),
col=c("black","darkslategray","beige",
"bisque","lightgreen","burlywood3",
"orange","chocolate","darkgoldenrod",
"darkgreen","dodgerblue","gold",
"yellow","darkorange","azure",
"purple"),
legend=FALSE,
par(bty = "n"),
colNA="gray")
title(main= "Substrate", adj = 0, line = 1, cex.main = 2)
#legend("bottomleft",legend=levels.substrate, fill = primary.colors(16), ncol = 3, cex = 0.45, bty = "n")
#legend("bottomleft",legend=levels.substrate, fill = primary.colors(16), ncol = 7, cex = 1, bty = "n") # legend not correct, need to figure out but values of raster are correct
legend("bottomleft",legend=levels.substrate,
fill = c("black","darkslategray","beige",
"bisque","lightgreen","burlywood3",
"orange","chocolate","darkgoldenrod",
"darkgreen","dodgerblue","gold",
"yellow","darkorange","azure",
"purple"),
ncol = 7, cex = 1, bty = "n", inset = c(0,0))
points(latlon, pch = 19, cex = 0.4, col ="black")
bf <- buffer(latlon, width = buffer_distance, dissolve = F)
plot(bf, add = TRUE)
legend("bottomright", inset = c(0.12,0.02), legend = c("sampling cell centroid","buffer"), col= "black", pch = c(20, 1), bty = "n")
dev.off()
#plotting test
# Set margin: Bottom, left, top & right
jpeg(file="Substrate_samplingcells_buffer_res200m_med2.jpeg", width = 1700, height = 900, units = "px")
#par(mar=c(1,1,1,1),lwd=0.25)
#split.screen(rbind(c(0,1,0.25,1), c(0,1,0,0.25)))
split.screen(c(2,1))
screen(1)
plot(rp,
axes = FALSE,
#ext = c(-915183, 4151749, 3296265, 5822348),
#ext = c(-955183, 4151749, 3296265, 5822348),
breaks=cuts,
#col=primary.colors(16),
col=c("black","darkslategray","beige",
"bisque","lightgreen","burlywood3",
"orange","chocolate","darkgoldenrod",
"darkgreen","dodgerblue","gold",
"yellow","darkorange","azure",
"purple"),
legend=FALSE,
par(bty = "n"),
colNA="gray")
title(main= "Substrate", adj = 0, cex.main = 2)
points(latlon, pch = 19, cex = 0.4, col ="black")
bf <- buffer(latlon, width = buffer_distance, dissolve = F)
plot(bf, add = TRUE)
# Legend
screen(2)
legend("bottomleft",legend=levels.substrate,
fill = c("black","darkslategray","beige",
"bisque","lightgreen","burlywood3",
"orange","chocolate","darkgoldenrod",
"darkgreen","dodgerblue","gold",
"yellow","darkorange","azure",
"purple"),
ncol = 7, cex = 1, bty = "n")
legend("bottomright", legend = c("sampling cell centroid","buffer"), col= "black", pch = c(20, 1), bty = "n")
dev.off()
#############################################################################
#### mean seafloor depth each sampling area (= buffer)
#import rasters
# doesn't import CRS, but when plot: same crs as seaconnect_geo (?)
B3 <- raster("EMODnet/Bathymetry/B3.asc")
B4 <- raster("EMODnet/Bathymetry/B4.asc")
C3 <- raster("EMODnet/Bathymetry/C3.asc")
C4 <- raster("EMODnet/Bathymetry/C4.asc")
D3 <- raster("EMODnet/Bathymetry/D3.asc")
D4 <- raster("EMODnet/Bathymetry/D4.asc")
# merge rasters to get full extent of mediterranean in one raster
bat <- merge(B3,B4,C3,C4,D3,D4)
#export raster
#writeRaster(bat, "EMODnet_Bathymetry_Med", format = "GTiff")
bat <- raster("EMODnet_Bathymetry_Med.tif")
projection(bat)
#create dataframe of coordinates
coord = data.frame(cbind(seaconnect_geo$Longitude, seaconnect_geo$Latitude))
colnames(coord) = c("lon","lat")
coordinates(coord) = c("lon", "lat")
### indiquer que les lat/lon est en projection google map
proj4string(coord) = CRS("+init=epsg:4326")
# faire parreil pour rasters
proj4string(bat) = CRS("+init=epsg:4326")
plot(bat)
points(coord)
# extract mean bathymetry of buffer zone around sampling cell coordinate
bat_sc_mean <- extract(bat, # raster layer
coord, # centroid coordinates of sampling cells
buffer = buffer_distance, # buffer size, units depend on CRS
fun = mean, # what value to extract
df=TRUE) # return a dataframe?
sc_bat <- cbind(seaconnect_geo[,c(1:2,6:7)], bat_sc_mean$layer)
names(sc_bat) <- c("ID","SamplingCell","Longitude","Latitude","mean_bat")
#export
write.table(sc_bat, file="seaconnect_bathymetry_tous.csv",sep=";",row.names=F)
# plot everything
library(RColorBrewer)
jpeg(file="Bathymetry_samplingcells_buffer.jpeg", width = 1700, height = 900, units = "px")
plot(bat,
col = rev(colorRampPalette(brewer.pal(9,"Blues"))(100)), # apply RColorBrewer "Blues" palette
axes = FALSE,
ext = c(-8,37, 29.5,47), # set extent, vector = c(xmin, xmax, ymin, ymax)
par(bty = "n"),
colNA="gray")
points(coord, pch = 19, cex = 0.1, col ="black")
# add buffer to map
buf <- buffer(coord, width = 0.5, dissolve =)
plot(buf, add = TRUE)
dev.off()
#"problem": I'm confused by the different CRS and corresponding units
# apply other projection system to raster, and see how this changes distances
bat <- raster("EMODnet_Bathymetry_Med.tif")
projection(bat) #NA
proj4string(bat) = proj4string(eusm)
projection(bat) #"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
projection(latlon) #idem
bat_sc_mean_t <- extract(bat, # raster layer
latlon, # centroid coordinates of sampling cells
buffer = buffer_distance, # buffer size, units depend on CRS
fun = mean, # what value to extract
df=TRUE) # return a dataframe?
sc_bat_t <- cbind(seaconnect_geo[,c(1:2,6:7)], bat_sc_mean$layer)
names(sc_bat_t) <- c("ID","SamplingCell","Longitude","Latitude","mean_bat")
jpeg(file="Bathymetry_samplingcells_buffer_t.jpeg", width = 1700, height = 900, units = "px")
plot(bat,
col = rev(colorRampPalette(brewer.pal(9,"Blues"))(100)), # apply RColorBrewer "Blues" palette
axes = FALSE,
ext = c(-8,37, 29.5,47), # set extent, vector = c(xmin, xmax, ymin, ymax)
par(bty = "n"),
colNA="gray")
title(main= "Bathymetry", adj = 0, line = -3.5 ,cex.main = 2)
points(coord, pch = 19, cex = 0.1, col ="black")
points(latlon, pch = 19, cex = 0.1, col ="red")
# add buffer to map
buf <- buffer(coord, width = buffer_distance, dissolve = F)
plot(buf, add = TRUE)
dev.off()
# gives same results! dus bathymetry values hierboven = OK. was gewoon bij het plotten dat de buffer afstanden aangepast moesten worden
#############################################################################
#### mean Photosynthetic Active Radiation (PAR) at seabed
PAR <- raster("EMODnet/PAR/T_MOSAIC_SURF_PAR_2005_2009.tif")
plot(PAR)
#proj4string(PAR) = proj4string(eusm)
projection(PAR) #"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
PAR_sc_mean <- extract(PAR, # raster layer
coord, # centroid coordinates of sampling cells
buffer = buffer_distance, # buffer size, units depend on CRS
fun = mean, # what value to extract
df=TRUE) # return a dataframe?
sc_PAR <- cbind(seaconnect_geo[,c(1:2,6:7)], PAR_sc_mean$T_MOSAIC_SURF_PAR_2005_2009)
names(sc_PAR) <- c("ID","SamplingCell","Longitude","Latitude","mean_PAR")
# add to bathymetry df
sc_bat_PAR <- cbind(sc_bat, PAR_sc_mean$T_MOSAIC_SURF_PAR_2005_2009)
names(sc_bat_PAR) <- c("ID","SamplingCell","Longitude","Latitude","mean_bat","mean_PAR")
write.table(sc_bat_PAR, "Seaconnect_bathymetry_PAR.csv", sep=";",row.names=F)
PAR_med <- crop(PAR, c(-8,37, 29.5,47))
jpeg(file="PAR_samplingcells_buffer_t.jpeg", width = 1700, height = 900, units = "px")
plot(PAR,
#col = rev(colorRampPalette(brewer.pal(9,"Blues"))(100)), # apply RColorBrewer "Blues" palette
axes = FALSE,
ext = c(-8,37, 29.5,47), # set extent, vector = c(xmin, xmax, ymin, ymax)
par(bty = "n"),
colNA="gray")
title(main= "Photosynthetic Active Radiation", adj = 0, line = -3.5 ,cex.main = 2)
points(coord, pch = 19, cex = 0.1, col ="black")
points(latlon, pch = 19, cex = 0.1, col ="red")
# add buffer to map
buf <- buffer(coord, width = buffer_distance, dissolve = F)
plot(buf, add = TRUE)
legend("bottomleft", inset = c(0,0.11),
legend = c("sampling cell centroid","buffer"), col= c("black", "black"), pch = c(20, 1), bty = "n", cex = 1.5)
dev.off()
\ No newline at end of file
###############################################################################
## Librairies
library("raster")