data-raw/generate_package_data.R
In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis
library(readxl)
library(IFNdyn)
library(maptools)
library(RODBC)
dir_devel = "~/Documents/KnowledgeTransfer/2018_2019_Modelitzacio_INIA_CTFC"
plotCodes = c("80001", "80002", "81073")
#Read IFN2 tree records
DBFdirIFN2 = "D:/Datasets/IFN/Sources/IFN2/DBF/"
treeDataIFN2 = IFNread::readPiesMayoresIFN2(prov = "08", DBFdir = DBFdirIFN2, subsetVars=T)
exampleTreeData = treeDataIFN2[treeDataIFN2$ID %in% plotCodes, c("ID", "Species", "N", "DBH", "H")]
usethis::use_data(exampleTreeData, overwrite = TRUE)
#example management data
examplePrescriptionData = data.frame(type = c("irregular", "irregular", "regular"),
thinning = c("above", "15-60/30-30/40-10", "below"),
thinningBA = c(10,10,10),
thinningHB = c(NA, NA, NA),
thinningPerc = c(20,20,20),
finalMeanDBH = c(NA, NA, 20),
finalPerc = c(NA, NA, "40-60-100"),
plantingSpecies = c(NA, NA, 24),
plantingDBH = c(NA,NA,5),
plantingHeight = c(NA,NA,5),
plantingDensity = c(NA,NA, 200),
finalPreviousStage = c(NA, NA, 0),
row.names =plotCodes, stringsAsFactors = F)
usethis::use_data(examplePrescriptionData, overwrite = TRUE)
df=as.data.frame(read_xls(paste0(dir_devel,"/ProvinciasCCAA.xls")))
provincias = df[,-1]
row.names(provincias) = df$CODIGO
# provincias = provincias[order(as.numeric(row.names(provincias))),]
usethis::use_data(provincias, overwrite = TRUE)
biomass_species_match = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Biomass models def.xlsx"), sheet = "Equiv_species"))
biomass_parameters = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Biomass models def.xlsx"), sheet = "Biomass_IFNdyn"))
volume_parameters = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Volume models def.xlsx"), sheet = "TODO", na=c("","-")))
#Remove column 'Obs'
volume_parameters = volume_parameters[,-which(names(volume_parameters)=="Obs")]
#Remove records with missing species (improve!)
volume_parameters = volume_parameters[!is.na(volume_parameters$ID_TAXON),]
usethis::use_data(biomass_species_match, biomass_parameters,
volume_parameters,
internal = TRUE, overwrite = TRUE)
## Example tree data
treeDataIFN23 = readRDS(paste0(dir_devel,"/Rdata/treeDataIFN23.rds"))
df = treeDataIFN23[,c("ID","Species","DBH2", "H2", "N2")]
names(df)<-c("ID","Species","DBH","H","N")
df = df[!is.na(df$DBH),]
df = df[!is.na(df$Species),]
a= basalArea(df, bySpecies=TRUE)
arel = 100*sweep(as.matrix(a), STATS= rowSums(a, na.rm=T), MARGIN=1, FUN="/")
species_selection = readRDS(paste0(dir_devel,"/Rdata/species_selection.rds"))
source(paste0(dir_devel,"/R scripts/maximum_values.R"))
species_selection$max_dinc10 = max_dinc10
species_selection$max_diameter = max_diameter
species_selection$max_height = max_height
usethis::use_data(species_selection, overwrite = TRUE)
RP_mapping = list("RP_Pinuspinea","RP_Pinussylvestris", "RP_Pinusuncinata", "RP_Pinuspinaster",
"RP_Pinushalepensis", "RP_Pinusnigra", NA,NA,
c("RP_Abiesalba","RP_Abiespinsapo"), NA, NA, NA,
NA, "RP_Quercusilex", "RP_Quercussuber", "RP_Quercusfaginea",
c("RP_Quercusrobur", "RP_Quercuspetraea"), c("RP_Quercuspyrenaica", "RP_Quercushumilis", "RP_Quercuscanariensis"),
NA,NA,NA,NA,
NA, "RP_Fagussylvatica",NA,NA,NA)
names(RP_mapping) = species_selection$names
usethis::use_data(RP_mapping, overwrite = TRUE)
df_list = vector("list", nrow(species_selection))
names(df_list) = as.character(species_selection$names)
for(i in 1:nrow(species_selection)) {
spp = unlist(strsplit(as.character(species_selection$species[i]), split = ","))
ccs = which(colnames(arel) %in% spp)
ids = row.names(arel)[which(rowSums(arel[,ccs, drop=FALSE])>99)]
if(i==5) {
ids = ids[substr(ids, 1,nchar(ids)-4)==17]
}
dfi = df[df$ID %in% ids,]
sd.D = tapply(dfi$DBH, INDEX=dfi$ID, FUN=length)
idj = names(which.max(sd.D)) # Choose plot with maximum SD(D)
dfout = dfi[dfi$ID==idj,]
# dfout = df[df$ID %in% sample(ids,1),]
row.names(dfout) = 1:nrow(dfout)
df_list[[i]] = dfout
}
usethis::use_data(df_list, overwrite = TRUE)
# Surface by SIG estrato (to calculate expansion factors for volume/biomass)
estratosIFN3 = vector("list", 50)
names(estratosIFN3) = provincias$Provincia
for(i in 1:50) {
ich = as.character(i)
if(i<10) ich = paste0("0",ich)
cat(".")
ch<-RODBC::odbcConnectAccess2007(paste0("D:/Datasets/IFN/Sources/IFN3/Sig/Sig_", ich,".accdb"))
st<-RODBC::sqlFetch(ch,"Estratos", as.is=TRUE)
estratosIFN3[[i]]<-data.frame(Area = st$Superficie, NPar = st$NPar, row.names = st$Estrato)
close(ch)
}
usethis::use_data(estratosIFN3, overwrite = TRUE)
#Load default demand by species and province (in annual m3 vcc)
defaultDemand = vector("list", 50)
names(defaultDemand) = provincias$Provincia
demand = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/DemandaAnual.xlsx"), sheet = "m3cc"))
demand$Espcodes[demand$Espcodes=="31.32"] ="31,32"
for(i in 1:50) {
df <-demand[demand$Provincia==as.character(i),-c(1,2)]
names(df)[c(1,2)]<-c("Species", "Name")
df$AnnualDemand = df$D2005_2015
df = df[,c("Species", "Name", "AnnualDemand")]
row.names(df)<-1:nrow(df)
defaultDemand[[i]] = df
}
usethis::use_data(defaultDemand, overwrite = TRUE)
defaultDemandCAT = defaultDemand[c("Barcelona", "Lleida","Girona","Tarragona")]
dx = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/DemandCatalonia.xlsx")))
dx = dx[-20,] #Remove total
row.names(dx) = dx$Species
defaultDemandCAT$Barcelona$AnnualDemand = dx[defaultDemandCAT$Barcelona$Name, "Barcelona"]
defaultDemandCAT$Lleida$AnnualDemand = dx[defaultDemandCAT$Lleida$Name, "Lleida"]
defaultDemandCAT$Girona$AnnualDemand = dx[defaultDemandCAT$Girona$Name, "Girona"]
defaultDemandCAT$Tarragona$AnnualDemand = dx[defaultDemandCAT$Tarragona$Name, "Tarragona"]
usethis::use_data(defaultDemandCAT, overwrite = TRUE)
#Load default prescription by species
defaultPrescription = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Prescription.xlsx"), sheet = "DefaultPrescription"))
defaultPrescription$Espcodes[defaultPrescription$Espcodes=="45.46"] ="45,46"
usethis::use_data(defaultPrescription, overwrite = TRUE)
defaultPrescriptionCAT = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Prescription.xlsx"), sheet = "DefaultPrescription_CAT"))
defaultPrescriptionCAT$Espcodes[defaultPrescriptionCAT$Espcodes=="45.46"] ="45,46"
usethis::use_data(defaultPrescriptionCAT, overwrite = TRUE)
#Load default product destiny by species
defaultProductsCAT = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Product_destination_by_species_CAT.xlsx"), sheet = "Destination"))
usethis::use_data(defaultProductsCAT, overwrite = TRUE)
#Creates a table with climate, topography and estratoIFN3 for plot data to be used in simulations
CCWD_year = readRDS(paste0(dir_devel,"/Rdata/CCWD_year_simulacion_FCC10_xy.rds"))
CSWD_year = readRDS(paste0(dir_devel,"/Rdata/CSWD_year_simulacion_FCC10_xy.rds"))
Rad_year = readRDS(paste0(dir_devel,"/Rdata/Rad_year_simulacion_FCC10_xy.rds"))
Temp_year = readRDS(paste0(dir_devel,"/Rdata/Temp_year_simulacion_FCC10_xy.rds"))
PET_year = readRDS(paste0(dir_devel,"/Rdata/PET_year_simulacion_FCC10_xy.rds"))
Prec_year = readRDS(paste0(dir_devel,"/Rdata/Prec_year_simulacion_FCC10_xy.rds"))
coords_topo_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/coords_topo_simulacion_FCC10_xy.rds"))
ID_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/ID_simulacion_FCC10_xy.rds"))
nplots = length(ID_simulacion_FCC10_xy)
plotDataIFN = data.frame(row.names=ID_simulacion_FCC10_xy)
plotDataIFN$Province = substr(ID_simulacion_FCC10_xy,1,nchar(ID_simulacion_FCC10_xy)-4)
plotDataIFN$CCWD = NA
plotDataIFN$CSWD = NA
plotDataIFN$Rad = NA
plotDataIFN$Temp = NA
plotDataIFN$Prec = NA
plotDataIFN$PET = NA
plotDataIFN$SWHC = NA
plotDataIFN$elevation = coords_topo_simulacion_FCC10_xy$elevation
plotDataIFN$slope = coords_topo_simulacion_FCC10_xy$slope
pb=txtProgressBar(1, length(ID_simulacion_FCC10_xy), style=3)
for(i in 1:length(ID_simulacion_FCC10_xy)) {
setTxtProgressBar(pb, i)
id = ID_simulacion_FCC10_xy[i]
plotDataIFN[i, "CCWD"] = mean(CCWD_year[id,])
plotDataIFN[i, "CSWD"] = mean(CSWD_year[id,])
plotDataIFN[i, "Temp"] = mean(Temp_year[id,])
plotDataIFN[i, "PET"] = mean(PET_year[id,])
plotDataIFN[i, "Rad"] = mean(Rad_year[id,])
plotDataIFN[i, "Prec"] = mean(Prec_year[id,])
}
rm(CCWD_year)
rm(CSWD_year)
rm(Temp_year)
rm(PET_year)
rm(Prec_year)
rm(Rad_year)
# Add SWHC to plot data
swhc_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/swhc_simulacion_FCC10_xy.rds"))
plotDataIFN[row.names(swhc_simulacion_FCC10_xy), "SWHC"] = swhc_simulacion_FCC10_xy$L13
# Add FCC to subset simulation plots
plotDataIFN3_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/plotDataIFN3_simulacion_FCC10_xy.rds"))
plotDataIFN$FCCARB = plotDataIFN3_simulacion_FCC10_xy$FCCARB
plotDataIFN$FCCTOT = plotDataIFN3_simulacion_FCC10_xy$FCCTOT
usethis::use_data(plotDataIFN, overwrite = TRUE)
# Read estrato from IFN3 sig (needed to scale output variables from plot to surface)
plotDataIFN$Estrato = NA
library(RODBC)
for(i in 1:50) {
cat(".")
ich = as.character(i)
if(i<10) ich = paste0("0",ich)
ch<-RODBC::odbcConnectAccess2007(paste0("D:/Datasets/IFN/Sources/IFN3/Sig/Sig_", ich,".accdb"))
parcPoly<-RODBC::sqlFetch(ch,"ParcPoly", as.is = TRUE)
parcPoly$IDs = paste0(i,parcPoly$Parcela)
sel = parcPoly$IDs %in% ID_simulacion_FCC10_xy
parcPolySel = parcPoly[sel, ]
plotDataIFN[parcPolySel$IDs, "Estrato"] = parcPolySel$Estrato
close(ch)
}
#NAs
sum(is.na(plotDataIFN$Estrato)) # Missing for 384 plots
usethis::use_data(plotDataIFN, overwrite = TRUE)
## Example neigbour data
lonlat<-coords_topo_simulacion_FCC10_xy[,c("lon","lat")]
lonlat<-as.matrix(lonlat)
examplePlotCoords = lonlat[plotCodes,]
usethis::use_data(examplePlotCoords, overwrite = TRUE)
#Identify protected areas
library(sf)
sp_latlon = SpatialPoints(cbind(coords_topo_simulacion_FCC10_xy$lon, coords_topo_simulacion_FCC10_xy$lat),
CRS("+proj=longlat +datum=WGS84"))
sf_latlon = sf::st_as_sf(sp_latlon)
nacionales_pb = sf::st_read("D:/Datasets/ProtectedAreas/Spain/ParquesNacionales/ParquesNacionales_P_B.shp")
nacionales_c = sf::st_read("D:/Datasets/ProtectedAreas/Spain/ParquesNacionales/ParquesNacionales_c.shp")
a=over(as(st_transform(sf_latlon,st_crs(nacionales_pb)), "Spatial"), as(nacionales_pb, "Spatial"))
b=over(as(st_transform(sf_latlon,st_crs(nacionales_c)), "Spatial"), as(nacionales_c, "Spatial"))
plotDataIFN$NationalPark = (!is.na(a$Shape_Area)) | (!is.na(b$Shape_Area))
enp = sf::st_read("D:/Datasets/ProtectedAreas/Spain/EspaciosNaturalesProtegidos/ENP_ES.shp")
c=over(as(st_transform(sf_latlon,st_crs(enp)), "Spatial"), as(enp, "Spatial"))
plotDataIFN$ProtectedArea = !is.na(c$area_ha)
table(plotDataIFN$ProtectedArea, plotDataIFN$NationalPark)
plotDataIFN$Managed = !plotDataIFN$NationalPark
usethis::use_data(plotDataIFN, overwrite = TRUE)
### Add regiones de Procedencia
library(sp)
library(sf)
coords_topo_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/coords_topo_simulacion_FCC10_xy.rds"))
sp_latlon = SpatialPoints(cbind(coords_topo_simulacion_FCC10_xy$lon, coords_topo_simulacion_FCC10_xy$lat),
CRS("+proj=longlat +datum=WGS84"))
plotRPIFN = getRegionesProcedencia(sp_latlon, RPfile = paste0(dir_devel,"/Rdata/RegionesProcedencia/All_species.rds"))
plotRPIFN = plotRPIFN[,-1] # remove codigo INE municipio
row.names(plotRPIFN) = row.names(coords_topo_simulacion_FCC10_xy)
usethis::use_data(plotRPIFN, overwrite = TRUE)
# provincias_pol = readShapePoly("/media/miquel/Elements/UNITATD/Recerca/Datasets/Limits/Sources/Spain/Provincias_ETRS89_30N/Provincias_ETRS89_30N.shp")
# provincias_pol@data = provincias_pol@data[,1,drop=FALSE]
# provincias_pol@data$Codigo = as.character(as.numeric(provincias_pol@data$Codigo))
# usethis::use_data(provincias_pol, overwrite = TRUE)
miquelcaceres/IFNdyn documentation built on Feb. 1, 2021, 10:55 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(readxl)
library(IFNdyn)
library(maptools)
library(RODBC)
dir_devel = "~/Documents/KnowledgeTransfer/2018_2019_Modelitzacio_INIA_CTFC"
plotCodes = c("80001", "80002", "81073")
#Read IFN2 tree records
DBFdirIFN2 = "D:/Datasets/IFN/Sources/IFN2/DBF/"
treeDataIFN2 = IFNread::readPiesMayoresIFN2(prov = "08", DBFdir = DBFdirIFN2, subsetVars=T)
exampleTreeData = treeDataIFN2[treeDataIFN2$ID %in% plotCodes, c("ID", "Species", "N", "DBH", "H")]
usethis::use_data(exampleTreeData, overwrite = TRUE)
#example management data
examplePrescriptionData = data.frame(type = c("irregular", "irregular", "regular"),
thinning = c("above", "15-60/30-30/40-10", "below"),
thinningBA = c(10,10,10),
thinningHB = c(NA, NA, NA),
thinningPerc = c(20,20,20),
finalMeanDBH = c(NA, NA, 20),
finalPerc = c(NA, NA, "40-60-100"),
plantingSpecies = c(NA, NA, 24),
plantingDBH = c(NA,NA,5),
plantingHeight = c(NA,NA,5),
plantingDensity = c(NA,NA, 200),
finalPreviousStage = c(NA, NA, 0),
row.names =plotCodes, stringsAsFactors = F)
usethis::use_data(examplePrescriptionData, overwrite = TRUE)
df=as.data.frame(read_xls(paste0(dir_devel,"/ProvinciasCCAA.xls")))
provincias = df[,-1]
row.names(provincias) = df$CODIGO
# provincias = provincias[order(as.numeric(row.names(provincias))),]
usethis::use_data(provincias, overwrite = TRUE)
biomass_species_match = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Biomass models def.xlsx"), sheet = "Equiv_species"))
biomass_parameters = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Biomass models def.xlsx"), sheet = "Biomass_IFNdyn"))
volume_parameters = as.data.frame(read_xlsx(paste0(dir_devel,"/Spanish Volume models def.xlsx"), sheet = "TODO", na=c("","-")))
#Remove column 'Obs'
volume_parameters = volume_parameters[,-which(names(volume_parameters)=="Obs")]
#Remove records with missing species (improve!)
volume_parameters = volume_parameters[!is.na(volume_parameters$ID_TAXON),]
usethis::use_data(biomass_species_match, biomass_parameters,
volume_parameters,
internal = TRUE, overwrite = TRUE)
## Example tree data
treeDataIFN23 = readRDS(paste0(dir_devel,"/Rdata/treeDataIFN23.rds"))
df = treeDataIFN23[,c("ID","Species","DBH2", "H2", "N2")]
names(df)<-c("ID","Species","DBH","H","N")
df = df[!is.na(df$DBH),]
df = df[!is.na(df$Species),]
a= basalArea(df, bySpecies=TRUE)
arel = 100*sweep(as.matrix(a), STATS= rowSums(a, na.rm=T), MARGIN=1, FUN="/")
species_selection = readRDS(paste0(dir_devel,"/Rdata/species_selection.rds"))
source(paste0(dir_devel,"/R scripts/maximum_values.R"))
species_selection$max_dinc10 = max_dinc10
species_selection$max_diameter = max_diameter
species_selection$max_height = max_height
usethis::use_data(species_selection, overwrite = TRUE)
RP_mapping = list("RP_Pinuspinea","RP_Pinussylvestris", "RP_Pinusuncinata", "RP_Pinuspinaster",
"RP_Pinushalepensis", "RP_Pinusnigra", NA,NA,
c("RP_Abiesalba","RP_Abiespinsapo"), NA, NA, NA,
NA, "RP_Quercusilex", "RP_Quercussuber", "RP_Quercusfaginea",
c("RP_Quercusrobur", "RP_Quercuspetraea"), c("RP_Quercuspyrenaica", "RP_Quercushumilis", "RP_Quercuscanariensis"),
NA,NA,NA,NA,
NA, "RP_Fagussylvatica",NA,NA,NA)
names(RP_mapping) = species_selection$names
usethis::use_data(RP_mapping, overwrite = TRUE)
df_list = vector("list", nrow(species_selection))
names(df_list) = as.character(species_selection$names)
for(i in 1:nrow(species_selection)) {
spp = unlist(strsplit(as.character(species_selection$species[i]), split = ","))
ccs = which(colnames(arel) %in% spp)
ids = row.names(arel)[which(rowSums(arel[,ccs, drop=FALSE])>99)]
if(i==5) {
ids = ids[substr(ids, 1,nchar(ids)-4)==17]
}
dfi = df[df$ID %in% ids,]
sd.D = tapply(dfi$DBH, INDEX=dfi$ID, FUN=length)
idj = names(which.max(sd.D)) # Choose plot with maximum SD(D)
dfout = dfi[dfi$ID==idj,]
# dfout = df[df$ID %in% sample(ids,1),]
row.names(dfout) = 1:nrow(dfout)
df_list[[i]] = dfout
}
usethis::use_data(df_list, overwrite = TRUE)
# Surface by SIG estrato (to calculate expansion factors for volume/biomass)
estratosIFN3 = vector("list", 50)
names(estratosIFN3) = provincias$Provincia
for(i in 1:50) {
ich = as.character(i)
if(i<10) ich = paste0("0",ich)
cat(".")
ch<-RODBC::odbcConnectAccess2007(paste0("D:/Datasets/IFN/Sources/IFN3/Sig/Sig_", ich,".accdb"))
st<-RODBC::sqlFetch(ch,"Estratos", as.is=TRUE)
estratosIFN3[[i]]<-data.frame(Area = st$Superficie, NPar = st$NPar, row.names = st$Estrato)
close(ch)
}
usethis::use_data(estratosIFN3, overwrite = TRUE)
#Load default demand by species and province (in annual m3 vcc)
defaultDemand = vector("list", 50)
names(defaultDemand) = provincias$Provincia
demand = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/DemandaAnual.xlsx"), sheet = "m3cc"))
demand$Espcodes[demand$Espcodes=="31.32"] ="31,32"
for(i in 1:50) {
df <-demand[demand$Provincia==as.character(i),-c(1,2)]
names(df)[c(1,2)]<-c("Species", "Name")
df$AnnualDemand = df$D2005_2015
df = df[,c("Species", "Name", "AnnualDemand")]
row.names(df)<-1:nrow(df)
defaultDemand[[i]] = df
}
usethis::use_data(defaultDemand, overwrite = TRUE)
defaultDemandCAT = defaultDemand[c("Barcelona", "Lleida","Girona","Tarragona")]
dx = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/DemandCatalonia.xlsx")))
dx = dx[-20,] #Remove total
row.names(dx) = dx$Species
defaultDemandCAT$Barcelona$AnnualDemand = dx[defaultDemandCAT$Barcelona$Name, "Barcelona"]
defaultDemandCAT$Lleida$AnnualDemand = dx[defaultDemandCAT$Lleida$Name, "Lleida"]
defaultDemandCAT$Girona$AnnualDemand = dx[defaultDemandCAT$Girona$Name, "Girona"]
defaultDemandCAT$Tarragona$AnnualDemand = dx[defaultDemandCAT$Tarragona$Name, "Tarragona"]
usethis::use_data(defaultDemandCAT, overwrite = TRUE)
#Load default prescription by species
defaultPrescription = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Prescription.xlsx"), sheet = "DefaultPrescription"))
defaultPrescription$Espcodes[defaultPrescription$Espcodes=="45.46"] ="45,46"
usethis::use_data(defaultPrescription, overwrite = TRUE)
defaultPrescriptionCAT = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Prescription.xlsx"), sheet = "DefaultPrescription_CAT"))
defaultPrescriptionCAT$Espcodes[defaultPrescriptionCAT$Espcodes=="45.46"] ="45,46"
usethis::use_data(defaultPrescriptionCAT, overwrite = TRUE)
#Load default product destiny by species
defaultProductsCAT = as.data.frame(read_xlsx(paste0(dir_devel,"/ManagementData/Product_destination_by_species_CAT.xlsx"), sheet = "Destination"))
usethis::use_data(defaultProductsCAT, overwrite = TRUE)
#Creates a table with climate, topography and estratoIFN3 for plot data to be used in simulations
CCWD_year = readRDS(paste0(dir_devel,"/Rdata/CCWD_year_simulacion_FCC10_xy.rds"))
CSWD_year = readRDS(paste0(dir_devel,"/Rdata/CSWD_year_simulacion_FCC10_xy.rds"))
Rad_year = readRDS(paste0(dir_devel,"/Rdata/Rad_year_simulacion_FCC10_xy.rds"))
Temp_year = readRDS(paste0(dir_devel,"/Rdata/Temp_year_simulacion_FCC10_xy.rds"))
PET_year = readRDS(paste0(dir_devel,"/Rdata/PET_year_simulacion_FCC10_xy.rds"))
Prec_year = readRDS(paste0(dir_devel,"/Rdata/Prec_year_simulacion_FCC10_xy.rds"))
coords_topo_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/coords_topo_simulacion_FCC10_xy.rds"))
ID_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/ID_simulacion_FCC10_xy.rds"))
nplots = length(ID_simulacion_FCC10_xy)
plotDataIFN = data.frame(row.names=ID_simulacion_FCC10_xy)
plotDataIFN$Province = substr(ID_simulacion_FCC10_xy,1,nchar(ID_simulacion_FCC10_xy)-4)
plotDataIFN$CCWD = NA
plotDataIFN$CSWD = NA
plotDataIFN$Rad = NA
plotDataIFN$Temp = NA
plotDataIFN$Prec = NA
plotDataIFN$PET = NA
plotDataIFN$SWHC = NA
plotDataIFN$elevation = coords_topo_simulacion_FCC10_xy$elevation
plotDataIFN$slope = coords_topo_simulacion_FCC10_xy$slope
pb=txtProgressBar(1, length(ID_simulacion_FCC10_xy), style=3)
for(i in 1:length(ID_simulacion_FCC10_xy)) {
setTxtProgressBar(pb, i)
id = ID_simulacion_FCC10_xy[i]
plotDataIFN[i, "CCWD"] = mean(CCWD_year[id,])
plotDataIFN[i, "CSWD"] = mean(CSWD_year[id,])
plotDataIFN[i, "Temp"] = mean(Temp_year[id,])
plotDataIFN[i, "PET"] = mean(PET_year[id,])
plotDataIFN[i, "Rad"] = mean(Rad_year[id,])
plotDataIFN[i, "Prec"] = mean(Prec_year[id,])
}
rm(CCWD_year)
rm(CSWD_year)
rm(Temp_year)
rm(PET_year)
rm(Prec_year)
rm(Rad_year)
# Add SWHC to plot data
swhc_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/swhc_simulacion_FCC10_xy.rds"))
plotDataIFN[row.names(swhc_simulacion_FCC10_xy), "SWHC"] = swhc_simulacion_FCC10_xy$L13
# Add FCC to subset simulation plots
plotDataIFN3_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/plotDataIFN3_simulacion_FCC10_xy.rds"))
plotDataIFN$FCCARB = plotDataIFN3_simulacion_FCC10_xy$FCCARB
plotDataIFN$FCCTOT = plotDataIFN3_simulacion_FCC10_xy$FCCTOT
usethis::use_data(plotDataIFN, overwrite = TRUE)
# Read estrato from IFN3 sig (needed to scale output variables from plot to surface)
plotDataIFN$Estrato = NA
library(RODBC)
for(i in 1:50) {
cat(".")
ich = as.character(i)
if(i<10) ich = paste0("0",ich)
ch<-RODBC::odbcConnectAccess2007(paste0("D:/Datasets/IFN/Sources/IFN3/Sig/Sig_", ich,".accdb"))
parcPoly<-RODBC::sqlFetch(ch,"ParcPoly", as.is = TRUE)
parcPoly$IDs = paste0(i,parcPoly$Parcela)
sel = parcPoly$IDs %in% ID_simulacion_FCC10_xy
parcPolySel = parcPoly[sel, ]
plotDataIFN[parcPolySel$IDs, "Estrato"] = parcPolySel$Estrato
close(ch)
}
#NAs
sum(is.na(plotDataIFN$Estrato)) # Missing for 384 plots
usethis::use_data(plotDataIFN, overwrite = TRUE)
## Example neigbour data
lonlat<-coords_topo_simulacion_FCC10_xy[,c("lon","lat")]
lonlat<-as.matrix(lonlat)
examplePlotCoords = lonlat[plotCodes,]
usethis::use_data(examplePlotCoords, overwrite = TRUE)
#Identify protected areas
library(sf)
sp_latlon = SpatialPoints(cbind(coords_topo_simulacion_FCC10_xy$lon, coords_topo_simulacion_FCC10_xy$lat),
CRS("+proj=longlat +datum=WGS84"))
sf_latlon = sf::st_as_sf(sp_latlon)
nacionales_pb = sf::st_read("D:/Datasets/ProtectedAreas/Spain/ParquesNacionales/ParquesNacionales_P_B.shp")
nacionales_c = sf::st_read("D:/Datasets/ProtectedAreas/Spain/ParquesNacionales/ParquesNacionales_c.shp")
a=over(as(st_transform(sf_latlon,st_crs(nacionales_pb)), "Spatial"), as(nacionales_pb, "Spatial"))
b=over(as(st_transform(sf_latlon,st_crs(nacionales_c)), "Spatial"), as(nacionales_c, "Spatial"))
plotDataIFN$NationalPark = (!is.na(a$Shape_Area)) | (!is.na(b$Shape_Area))
enp = sf::st_read("D:/Datasets/ProtectedAreas/Spain/EspaciosNaturalesProtegidos/ENP_ES.shp")
c=over(as(st_transform(sf_latlon,st_crs(enp)), "Spatial"), as(enp, "Spatial"))
plotDataIFN$ProtectedArea = !is.na(c$area_ha)
table(plotDataIFN$ProtectedArea, plotDataIFN$NationalPark)
plotDataIFN$Managed = !plotDataIFN$NationalPark
usethis::use_data(plotDataIFN, overwrite = TRUE)
### Add regiones de Procedencia
library(sp)
library(sf)
coords_topo_simulacion_FCC10_xy = readRDS(paste0(dir_devel,"/Rdata/coords_topo_simulacion_FCC10_xy.rds"))
sp_latlon = SpatialPoints(cbind(coords_topo_simulacion_FCC10_xy$lon, coords_topo_simulacion_FCC10_xy$lat),
CRS("+proj=longlat +datum=WGS84"))
plotRPIFN = getRegionesProcedencia(sp_latlon, RPfile = paste0(dir_devel,"/Rdata/RegionesProcedencia/All_species.rds"))
plotRPIFN = plotRPIFN[,-1] # remove codigo INE municipio
row.names(plotRPIFN) = row.names(coords_topo_simulacion_FCC10_xy)
usethis::use_data(plotRPIFN, overwrite = TRUE)
# provincias_pol = readShapePoly("/media/miquel/Elements/UNITATD/Recerca/Datasets/Limits/Sources/Spain/Provincias_ETRS89_30N/Provincias_ETRS89_30N.shp")
# provincias_pol@data = provincias_pol@data[,1,drop=FALSE]
# provincias_pol@data$Codigo = as.character(as.numeric(provincias_pol@data$Codigo))
# usethis::use_data(provincias_pol, overwrite = TRUE)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.