Commit de25808c authored by peguerin's avatar peguerin
Browse files

clean scripts folder

parent 572f747b
###############################################################################
## Librairies
library("raster")
library("rgdal")
library("ncdf4")
## 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")
fichier="01-infos/salinity/MEDSEA_PHY_006_004_SSS_jan1987-nov2016_monthly.nc"
br <- brick(fichier)
##convertir les lat/lon dans la projection du raster "rp"
latlon = spTransform(latlon,proj4string(br))
buffer_distance=13000
tt=extract(br,latlon, buffer=buffer_distance,fun=mean,normalizeWeights=FALSE,weigths=TRUE,factors=TRUE)
january=tt[,seq(1,(359-10),12)]
february=tt[,seq(2,(359-9),12)]
march=tt[,seq(3,(359-8),12)]
april=tt[,seq(4,(359-7),12)]
may=tt[,seq(5,(359-6),12)]
june=tt[,seq(6,(359-5),12)]
july=tt[,seq(7,(359-4),12)]
august=tt[,seq(8,(359-3),12)]
september=tt[,seq(9,(359-2),12)]
october=tt[,seq(10,(359-1),12)]
november=tt[,seq(11,359,12)]
december=tt[,seq(12,359-11,12)]
year=list(january,february,march,april,may,june,july,august,september,october,november,december)
monthlysalinity.frame <- as.data.frame(matrix(,ncol=(36),nrow=0))
for(i in seq(1,dim(tt)[1])){
yeartemp=c()
for(mois in seq(1,12)) {
month_min=min(as.numeric(as.data.frame(year[mois])[i,]),na.rm=TRUE)
month_max=max(as.numeric(as.data.frame(year[mois])[i,]),na.rm=TRUE)
month_mean=mean(as.numeric(as.data.frame(year[mois])[i,]),na.rm=TRUE)
yeartemp=c(yeartemp,month_min,month_max,month_mean)
}
monthlysalinity.frame=rbind(monthlysalinity.frame,yeartemp)
}
mois_nom=c("january","february","march","april","may","june","july","august","september","october","november","december")
value_nom=c("min","max","mean")
col_nom=c()
for(m in mois_nom){
for(n in value_nom) {
mois_value=paste(m,n,sep="_")
nom_value=paste("SSS",as.character(mois_value),sep="_")
col_nom=c(col_nom,as.character(nom_value))
}
}
colnames(monthlysalinity.frame)=col_nom
monthlysalinity.coo=cbind(coo,monthlysalinity.frame)
write.table(monthlysalinity.coo, file="03-results/all_monthly_salinity.tsv",quote=F,sep="\t",col.names=T,row.names=F)
\ No newline at end of file
###############################################################################
## Librairies
library("raster")
library("rgdal")
library("ncdf4")
###############################################################################
###############################################################################
## load data
##salinity
fichier="01-infos/MEDSEA_REANALYSIS_PHYS_006_004_6Wto17E_0to210m_1m_vosaline_1987-2015.nc"
fichier="01-infos/MEDSEA_REANALYSIS_PHYS_006_004_6Wto17E_0to210m_1m_votemper_1987-2015.nc"
## coordinates
coof="01-sample_coords/GeneticData_usedInAnalyses_adjustedGPS.txt"
coo=read.table(coof,sep="\t",header=T,stringsAsFactors = FALSE)
labels=as.character(coo[,1])
latlon = data.frame(lon=as.numeric(coo$long_adjusted) ,lat=as.numeric(coo$lat_adjusted))
## shit lattitude to get a grid point
latlon[c(24:34,39, 40),]$lat=latlon[c(24:34,39, 40),]$lat-0.5
coordinates(latlon) = c("lon", "lat")
### indiquer que les lat/lon est en projection google map
proj4string(latlon) = CRS("+init=epsg:4326")
br <- nc_open(fichier)
brdepth = ncvar_get(br, "depth")
january=seq(1,(348-11),12)
february=seq(2,(348-10),12)
march=seq(3,(348-9),12)
april=seq(4,(348-8),12)
may=seq(5,(348-7),12)
june=seq(6,(348-6),12)
july=seq(7,(348-5),12)
august=seq(8,(348-4),12)
september=seq(9,(348-3),12)
october=seq(10,(348-2),12)
november=seq(11,348-1,12)
december=seq(12,348,12)
## salinities ALL
repro_period=cbind(january,february,march,april,may,june,july,august,september,october,november,december)
depth_range=which(brdepth <300)
salinities=array(, dim=c(dim(latlon)[1],348,26))
for(LVL in depth_range){
br <- brick(fichier,level=LVL)
LVL_ma_var=extract(br,latlon)
salinities[,,LVL]=LVL_ma_var[,repro_period]
}
temperatures=array(, dim=c(dim(latlon)[1],348,26))
for(LVL in depth_range){
br <- brick(fichier,level=LVL)
LVL_ma_var=extract(br,latlon)
temperatures[,,LVL]=LVL_ma_var[,repro_period]
}
write.table(salt_temp, file="03-results/ADJUSTED_depth_repro_temperature_salinity.tsv",quote=F,sep="\t",col.names=T,row.names=F)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment