Commit 0c650664 authored by Boulangeat Isabelle's avatar Boulangeat Isabelle
Browse files

extracted climate data

parent 2354e5ed
File added
File added
install.packages("remotes")
remotes::install_github("jalvesaq/colorout")
??inout
install.packages("splancs")
?reclass
?breaks
?break
?addPolygons
library(leaflet)
?addPolygons
?colorBins
?colorBin
?addLegend
library(sf)
install.packages("sf")
library(sf)
?sf_interset
?sf_intersect
?intersect
?st_intersection
install.packages("rgeos")
library(rgeos)
?gIntersect
?gOverlap
?gOverlaps
?pairs
install.packages("queyras")
shiny::runApp('PROJETS/SOCBIO_Yves/ASONB/gitlab-folder/exploHumanImpactData')
?write.csv
library(FD)
?divc
install.pacakges("rgdal")
install.packages("rgdal")
library(rgdal)
install.packages("adespatial")
library(adespatial)
library(adespatial)
library(sf)
setwd("~/Documents/PROJETS/HISTFUNC/H3")
library(remake)
remake::make("herbivores")
remake::make("herbivores")
?pairs
?extract
install.packages("tidyverse")
shiny::runApp('~/Documents/PROJETS/SOCBIO_Yves/ASONB/gitlab-folder/exploHumanImpactData')
install.packages("PostgreSQL")
install.packages("RPostgres")
install.packages("rpostgis")
install.packages("RPostgreSQL")
sapply(list.files("R_fct"), function(x)source(paste0("R_fct/",x)))
setwd("~/Documents/PROJETS/ZAA/AS/git_as_biomass")
sapply(list.files("R_fct"), function(x)source(paste0("R_fct/",x)))
require(RPostgreSQL)
data = readRDS("data.rds")
head(data$sites)
dat_h_veg = merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site")
str(dat_h_veg)
summary(dat_h_veg$hmean)
library(dplyr)
# library(reshape)
library(tidyr)
library(ggplot2)
dat_long <- dat_h_veg %>%
select(hmean, hmad, ref_typoveg) %>%
drop_na() %>%
group_by(ref_typoveg) %>%
summarize(hauteur=mean(hmean), var_inter=mad(hmean), var_intra = mean(hmad)) %>%
gather(stat, value, hauteur:var_intra)
ggplot(dat_long, aes(fill=stat, x=ref_typoveg, y=value)) +
geom_bar(position = "dodge", stat = "identity") +
facet_wrap(~stat, ncol = 1, scales = "free")
head(dat_long)
head(dat_h_veg)
pl <- ggplot(dat_h_veg, aes(x = date_releve, y = hmean)) +
geom_line(aes(color = ref_site), show.legend = FALSE) +
facet_wrap(~ref_typoveg)
pl + theme(legend.position = "none")
head(dat_h_veg)
ggplot(dat_h_veg, aes(x = ref_typoveg, y = hmean)) +
geom_boxplot(stat = "identidy")
ggplot(dat_h_veg, aes(x = ref_typoveg, y = hmean)) +
geom_boxplot(stat = "identity")
ggplot(dat_h_veg, aes(x = ref_typoveg, y = hmean)) +
geom_boxplot()
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean), y = hmean)) +
geom_boxplot()
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
geom_boxplot()
summary(dat_h_veg)
dat_h_veg[which(is.na(dat_h_veg$hmad)),]
dat_h_veg[which(is.na(dat_h_veg$hmad)),"ref_site"]
unique(dat_h_veg[which(is.na(dat_h_veg$hmad)),"ref_site"])
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
geom_boxplot() +
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation",
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
geom_boxplot() +
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")
head(data)
head(data$h)
head(data$sites)
tail(data$h)
summary(data$h$hmean)
data$h$ref_releve[which(is.na(data$h$hmean))]
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hsd, na.rm=TRUE), y = hsd)) +
geom_boxplot() +
labs(y="Hétérogénéité", x="Type vegetation")
summary(data_h_veg$ref_typoveg)
summary(dat_h_veg$ref_typoveg)
table(dat_h_veg$ref_typoveg)
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hsd, na.rm=TRUE), y = hsd)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hsd, na.rm=TRUE), y = hsd/hmean)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hsd)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
geom_points()
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
geom_point()
ggplot(dat_h_veg, aes(x = date_releve, y = hsd)) +
geom_line(aes(color = ref_site), show.legend = FALSE) +
facet_wrap(~ref_typoveg)
ggplot(dat_h_veg, aes(x = date_releve, y = hmad)) +
geom_line(aes(color = ref_site), show.legend = FALSE) +
facet_wrap(~ref_typoveg)
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad, na.rm=TRUE), y = hmad)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad/hmean)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad/hmean)) +
geom_boxplot() +
labs(y="Variation intra de hauteur (stand. moy.)", x="Type vegetation")
head(dat_h_veg)
head(dat_h_veg)
dat_h_veg[which(dat_h_veg$ref_site=="BELRIV01"),]
sapply(dat_h_veg$ref_releve, function(x){
paste(strplit(x, "_")[-2])
})
sapply(dat_h_veg$ref_releve, function(x){
paste(strsplit(x, "_")[-2])
})
sapply(dat_h_veg$ref_releve, function(x){
paste(unlist(strsplit(x, "_")[-2]))
})
sapply(dat_h_veg$ref_releve, function(x){
paste(unlist(strsplit(x, "_"))[-2])
})
lapply(dat_h_veg$ref_releve, function(x){
paste(unlist(strsplit(x, "_"))[-2])
})
unlist(lapply(dat_h_veg$ref_releve, function(x){
paste(unlist(strsplit(x, "_"))[-2], sep="_")
}))
sapply(dat_h_veg$ref_releve, function(x){
paste(unlist(strsplit(x, "_"))[-2], sep="_")
})
x = dat_h_veg$ref_releve[1]
strsplit(x, "_")
strsplit(x, "_")[-2]
strsplit(x, "_")[[1]][-2]
sapply(dat_h_veg$ref_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], sep="_")
})
lapply(dat_h_veg$ref_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], sep="_")
})
x
strsplit(x, "_")[[1]][-2]
paste(strsplit(x, "_")[[1]][-2], sep="_")
paste(strsplit(x, "_")[[1]][-2], sep="_", collapse = TRUE)
x1 = strsplit(x, "_")[[1]][-2]
paste(x1, sep="_")
x1
paste(x1, collapse="_")
x1 =
paste(strsplit(x, "_")[[1]][-2], collapse="_")
x1 =
paste(strsplit(x, "_")[[1]][-2], collapse="_")
paste(strsplit(x, "_")[[1]][-2], collapse="_")
lapply(dat_h_veg$ref_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], collapse="_")
})
dat_h_veg$releveL1L2 =
unlist(lapply(dat_h_veg$ref_releve, function(x){
paste(strsplit(x, "_")[[1]][-2], collapse="_")
}))
head(dat_h_veg)
data = readRDS("data.rds")
#==============================
## dataframe ##
#===============================
head(data$sites)
head(data$h)
dat_h_veg = merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site")
str(dat_h_veg)
data = readRDS("data.rds")
#==============================
## dataframe ##
#===============================
head(data$sites)
head(data$h)
data = readRDS("data.rds")
#==============================
## dataframe ##
#===============================
head(data$sites)
head(data$h)
dat_h_veg = merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site")
str(dat_h_veg)
summary(dat_h_veg$hmean)
dat_h_veg = unique(merge(data$h, data$sites[,c("id_site", "ref_typoveg")], by.x = "ref_site", by.y ="id_site"))
length(unique(dat_h_veg$releve))
dim(dat_h_veg)
summary(dat_h_veg$hmean)
head(data$otable)
data$otable[which(is.na(data$otable$hauteur)),]
summary(dat_h_veg$hmean)
dat_long <- dat_h_veg %>%
select(hmean, hmad, ref_typoveg) %>%
drop_na() %>%
group_by(ref_typoveg) %>%
summarize(hauteur=mean(hmean), var_inter=mad(hmean), var_intra = mean(hmad)) %>%
gather(stat, value, hauteur:var_intra)
ggplot(dat_long, aes(fill=stat, x=ref_typoveg, y=value)) +
geom_bar(position = "dodge", stat = "identity") +
facet_wrap(~stat, ncol = 1, scales = "free")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
geom_boxplot() +
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")
pl <- ggplot(dat_h_veg, aes(x = date_releve, y = hmean)) +
geom_line(aes(color = ref_site), show.legend = FALSE) +
facet_wrap(~ref_typoveg)
pl + theme(legend.position = "none")
head(dat_h_veg)
#==============================
## effet obs ##
#===============================
dat_h_veg[which(is.duplicated(dat_h_releve)),]
#==============================
## effet obs ##
#===============================
dat_h_veg[which(is.duplicated(dat_h_veg$releve)),]
?is.duplicated
??is.duplicated
#==============================
## effet obs ##
#===============================
dat_h_veg[which(duplicated(dat_h_veg$releve)),]
duplicated(dat_h_veg$releve
duplicated(dat_h_veg$releve)
duplicated(dat_h_veg$releve)
### sd vs mad
ggplot(dat_h_veg, aes(x = hmad, y = hsd)) +
geom_point()
### dans le temps
ggplot(dat_h_veg, aes(x = date_releve, y = hmad)) +
geom_line(aes(color = ref_site), show.legend = FALSE) +
facet_wrap(~ref_typoveg)
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad, na.rm=TRUE), y = hmad)) +
geom_boxplot() +
labs(y="Variation intra de hauteur", x="Type vegetation")
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmad/hmean, na.rm=TRUE), y = hmad/hmean)) +
geom_boxplot() +
labs(y="Variation intra de hauteur (stand. moy.)", x="Type vegetation")
calage = read.csv("../PNE_calage_reformat.csv")
calage = read.csv("../PNE_calage_reformat.csv")
head(calage)
calage = read.csv("../PNE_calage_reformat.csv", sep= ";")
head(calage)
str(calage)
calage %>% group_by(site) %>% summarize(hmean = mean(hauteur))
str(calage)
calage = read.csv("../PNE_calage_reformat.csv", sep= ";", dec = ",")
str(calage)
calage %>% group_by(site) %>% summarize(hmean = mean(hauteur))
calage %>% group_by(site) %>% summarize(hmean = mean(hauteur, na.rm=TRUE))
cal_mean = calage %>% group_by(site) %>% summarize(hmean = mean(hauteur, na.rm=TRUE))
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_boxplot()
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_point()
head(calage)
calage$releve = paste0(calage$site, calage$Releveur)
cal_mean = calage %>% group_by(releve) %>% summarize(hmean = mean(hauteur, na.rm=TRUE))
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_point()
ggplot(cal_mean, aes(x = releve, y = hmean)) +
geom_point()
ggplot(cal_mean, aes(x = releve, y = hmean)) +
geom_point() +
geom_boxplot()
ggplot(cal_mean, aes(x = releve, y = hmean)) +
geom_point()
cal_mean
cal_mean = calage %>% group_by(releve) %>% summarize(site = site, hmean = mean(hauteur, na.rm=TRUE))
cal_mean
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_point() +
geom_boxplot()
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_boxplot() +
geom_point()
cal_mean
cal_mean %>% group_by(site) %>% summarize(hmad = hmad(hmean))
cal_mean %>% group_by(site)
calage %>% group_by(releve)
cal_mean = unique(calage %>% group_by(releve) %>% summarize(site = site, hmean = mean(hauteur, na.rm=TRUE)))
ggplot(cal_mean, aes(x = site, y = hmean)) +
geom_boxplot() +
geom_point()
cal_mean %>% group_by(site)
cal_mean %>% group_by(site) %>% summarize(hmad = hmad(hmean))
cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean))
cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean)/mean(hmean))
cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean))
cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean))
ggplot(cal_mean, aes(x = reorder(site, hmean), y = hmean)) +
geom_boxplot() +
geom_point()
cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean))
obs_error = cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean))
ggplot(obs_error, aes(x = site, y = hmad_st)) +
geom_bar()
ggplot(obs_error, aes(x = site, y = hmad_st)) +
geom_bar(stat = "identity")
obs_error = cal_mean %>% group_by(site) %>% summarize(hmad = mad(hmean), hmad_st = mad(hmean)/mean(hmean), hm = mean(hmean))
ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
geom_bar(stat = "identity")
ggplot(cal_mean, aes(x = reorder(site, hmean), y = hmean)) +
geom_boxplot() +
geom_point()
ggplot(obs_error, aes(x = reorder(site, hm), y = hmad_st)) +
geom_bar(stat = "identity")
ggplot(cal_mean, aes(x = reorder(site, hmean), y = hmean)) +
geom_boxplot() +
geom_point()
ggplot(dat_h_veg, aes(x = reorder(ref_typoveg, hmean, na.rm=TRUE), y = hmean)) +
geom_boxplot() +
labs(y="Hauteur moyenne (est. biomasse)", x="Type vegetation")
head(dat_h_veg)
table(ref_typoveg)
table(dat_h_veg$ref_typoveg)
dat_h_veg[which(dat_h_veg$ref_typoveg=="EBOU"),]
dat_h_veg[which(dat_h_veg$ref_typoveg=="HUM"),]
head(cal_mean)
cal_mean
extract_ncdat <- function(path, variables, NoPs){
nc.read = nc_open(path)
time<-ncvar_get(nc.read,"time")
tunits<-ncatt_get(nc.read,"time",attname="units")
tustr<-strsplit(tunits$value, " ")
dates<-as.Date(time,origin=unlist(tustr)[3])
nc.data = lapply(variables, function(x) ncvar_get(nc.read,x)[unique(NoPs),])
nc_close(nc.read)
return(list(data = abind(nc.data, along = 3), dates = dates))
}
calc_meteo_variables <- function(path_data_allslopes, year, days, NoPs, periods, tbase = 5.5){
# year = as.numeric(test$year[1])
# days = test$date_releve
# NoPs = test$NoP
variables = c("RN_ISBA", "TS_ISBA", "WSN_T_ISBA", "DSN_T_ISBA", "RAINF_ISBA", "TALB_ISBA")
# periods = c(30,60)
# library required
require(tidync)
require(ncmeta)
require(tidyr)
require(ncdf4)
require(ncdf4.helpers)
require(abind)
# check filepath dataset
if(!(length(list.files(path_data_allslopes)) > 0) ) stop("not able to connect to data path or empty path")
dat1 = extract_ncdat(paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(year),"080106_",as.character(year+1),"080106.nc"), variables, NoPs)
dat2 = extract_ncdat(paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(year-1),"080106_",as.character(year),"080106.nc"), variables, NoPs)
datAll = abind(dat2$data, dat1$data, along = 2)
datesAll = c(dat2$dates, dat1$dates)
output = array(data = NA, dim = list(length(days),8, length(periods)), dimnames = list(1:length(days), c("snowfree", "cumgdd", "frost", "gddspeed", "rainfall", "radiations", "albedo", "snowstock"), periods))
# releves loop
for (rel in 1:length(days)) {
# rel = 1
nop = which(NoPs[rel]%in%unique(NoPs))
d = which(days[rel]==datesAll)
for (p in 1:length(periods)){
# p = 1
temp = datAll[nop,(d-periods[p]):d, which(variables == "TS_ISBA")]-273.15
snow = datAll[nop,(d-periods[p]):d, which(variables == "DSN_T_ISBA")]
snowfree = snow==0
gdd = (temp - tbase)
gdd[gdd<0]=0
#### snowfree
output[rel, "snowfree", p] = sum(snowfree)
#### cumgdd (snowfree)
output[rel, "cumgdd", p] = sum(gdd[snowfree])
#### snowfreefrost
output[rel, "frost", p] = sum(which(temp<0 & snowfree))
#### gddspeed
output[rel, "gddspeed", p] = mean(gdd[gdd>0])
#### rainfall
output[rel, "rainfall", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "RAINF_ISBA")])
#### radiations
output[rel, "radiations", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "RN_ISBA")])
#### albedo
output[rel, "albedo", p] = sum(datAll[nop,(d-periods[p]):d, which(variables == "TABL_ISBA")])
#### snowstock
output[rel, "snowstock", p] = max(datAll[nop,(d-periods[p]):d, which(variables == "WSN_T_ISBA")])
}
}
return(output)
}
#=============================================
#--- build NoP raster functions
#=============================================
NoP_from_values <- function(x, safran_classes_complet, alti_field="ZS", slope_field="slope", aspect_field="aspect", massif_field="massif_num"){
num_point = safran_classes_complet[which(((safran_classes_complet[,alti_field] == x[alti_field] & safran_classes_complet[,slope_field] == x[slope_field]) & safran_classes_complet[,aspect_field]==x[aspect_field]) & safran_classes_complet[,massif_field]==x[massif_field]), "Number_of_points"]
if(length(num_point)==0) num_point=9999
return(num_point)
}
NoP_from_rasters <- function(r_alti, r_slope, r_aspect, r_massif, safran_classes_complet, parallel=FALSE, nclust=2){
if(resolution(r_alti)!=resolution(r_slope) || resolution(r_alti)!=resolution(r_aspect) || resolution(r_alti)!=resolution(r_massif)) stop("pb resolution")
if(projection(r_alti)!=projection(r_slope) || projection(r_alti)!=projection(r_aspect) || projection(r_alti)!=projection(r_aspect)) stop("pb proj")
crocus_alti <- cut(r_alti, breaks = seq(0, 5200, by = 300) , include.lowest = TRUE)
crocus_aspect <- cut(r_aspect, breaks = c(-1,seq(0, 360, by = 45)) , include.lowest = TRUE)
crocus_slope <- cut(r_slope, breaks = c(0, 10, 30, 90) , include.lowest = TRUE)
crocus_stack <- stack(crocus_alti, crocus_slope, crocus_aspect, crocus_massif_raster)
names(crocus_stack) <- c("altitude", "slope", "aspect", "massif")
corresp <- function(x, na.rm){
if(sum(is.na(x))==0){
num_point = safran_classes_complet[which(((sf4$cat_raster_alti == x[1] & safran_classes_complet$cat_raster_aspect == x[3]) & safran_classes_complet$cat_raster_slope==x[2]) & safran_classes_complet$massif_num==x[4]), "Number_of_points"]
}else num_point=NA
if(length(num_point)==0) num_point=9999
return(num_point)
}
if(parallel==TRUE){
library(doParallel)
beginCluster(n=nclust, type="SOCK")
cl = getCluster()
clusterExport(cl, "safran_classes_complet", envir=environment())
NoP = clusterR(crocus_stack, stackApply, args = list(indices=rep(1, 4), fun = corresp, na.rm=TRUE))
endCluster()
}else{
NoP <- stackApply(crocus_stack, indices=rep(1, 4),fun = corresp, na.rm=TRUE)
}
return(NoP)
}
library(tidync)
library(ncmeta)
library(tidyr)
library(raster)
library(sp)
#=============================================
#--- Data input
#=============================================
sites <- readRDS("data_sites.rds
sites <- sites[which(!is.na(sites$x)),]
coordinates(sites)= cbind(sites$y, sites$x)
proj4string(sites) = CRS("+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
NoP_raster = raster()
plot(NoP_raster)
points(sites)
spTransform()
extract()
# path_data_allslopes = "/Volumes/infogeo/Meteo_France/SAFRAN_montagne-CROCUS_2019/alp_allslopes"
path_data_allslopes = "/Volumes/ISA-RESEARCH/alp_allslopes"
safran_classes <- read.table("safran_classes.txt", h=T, sep="\t")
#=============================================
#--- Extract function
#=============================================
# list of variables:
#
extract_meteo <- function(path_data_allslopes, years, variables, NoP){
# library required
require(tidync)
require(ncmeta)
require(tidyr)
# check filepath dataset
if(!(length(list.files(path_data_allslopes)) > 0) ) stop("not able to connect to data path or empty path")
# temporal loop
for (y in years){
# read data and stock for a year y
pro_file = paste0(path_data_allslopes, "/reanalysis/pro/PRO_",as.character(y),"080106_",as.character(y+1),"080106.nc")
src <- tidync(pro_file)
data_y <- src %>% activate("D0,D3") %>% hyper_filter(Number_of_points = Number_of_points %in% NoP) %>% hyper_tibble(select_var = variables)
# TODO: compute values from data_y
}
# build output
# TODO!!
return(output)
}
File added
File added
#=============================================
#--- Sites, extration NoP_raster
# inputs: sites, NoP_raster
# outputs: NoP for each site
#=============================================
library(tidyr)
library(raster)
library(sp)
sites <- readRDS("data_sites.rds")
head(sites)
summary(sites$x)
sites1 <- sites[which(sites$x<10),]
coordinates(sites1)= c("x","y")
proj4string(sites1) = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
sites2 <- sites[which(sites$x>10),]
coordinates(sites2)= c("x","y")