diff --git a/leer_climas_paleoclim.R b/leer_climas_paleoclim.R index 32d51ec..7d8527f 100644 --- a/leer_climas_paleoclim.R +++ b/leer_climas_paleoclim.R @@ -8,6 +8,7 @@ setwd("C:\\Users\\sara.varela\\Documents\\CIENCIAS\\mariana_araucaria") costa<- readShapePoly ("ne_10m_land.shp") # load paleoclim layers +setwd ("C:\\Users\\sara.varela\\Documents\\CIENCIAS\\mariana_araucaria\\paleoclim_org") carpetas<- list.files () i<- 8 for (i in 1:8){ @@ -56,7 +57,7 @@ datos<- datos [complete.cases(datos), ] summary (datos) # set percentile to drop out -x1<- 0.07 +x1<- 0.01 min_t<- quantile (datos$bio_1, x1) max_t<- quantile (datos$bio_1, 1-x1) @@ -102,48 +103,56 @@ eval(parse(text=foo)) # build raw figure -tiff ("lala.tif", width = 500, height = 1000, units = "px") +tiff ("PaleoClim_envelope.tif", width = 500, height = 1000, units = "px") par (mfrow=c(2,4)) plot (map_1, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="present", xlim=c(-55, -30), ylim=c(-45, 0)) + main="present", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) plot (araucaria, add=T) points (coord2) plot (map_2, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="0.3-4.2 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="0.3-4.2 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X0.3_4.2 ==1,2:3 ]) plot (map_3, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="4.2-8.3 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="4.2-8.3 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X4.2_8.3 ==1, 2:3 ]) plot (map_4, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="8.3-11.7 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="8.3-11.7 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X8.3_11.7 ==1,2:3 ]) plot (map_5, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="11.7-12.9 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="11.7-12.9 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X11.7_12.9 ==1, 2:3 ]) plot (map_6, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="12.9-14.7 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="12.9-14.7 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X12.9_14.7 ==1, 2:3 ]) plot (map_7, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="14.7-17 kybp", xlim=c(-55, -30), ylim=c(-45, 0)) + main="14.7-17 kybp", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$X14.7_17 ==1, 2:3 ]) plot (map_8, col=c("#00000000","#00006090"), legend=F, axes=F, box=F, - main="LGM", xlim=c(-55, -30), ylim=c(-45, 0)) + main="LGM", xlim=c(-55, -30), ylim=c(-45, -20)) plot (wrld_simpl, add=T) points (polen [polen$lgm ==1, 2:3 ]) dev.off() +stabi<- (map_1 + map_2 + map_3 + map_4 + map_5 + map_6 + map_7 + map_8) + +pdf ("stability_paleoclim.pdf") +plot (stabi) +plot (araucaria, add=T) +points (coord2) +plot (wrld_simpl, add=T) +dev.off() \ No newline at end of file diff --git a/leer_ecoclimate_data.R b/leer_ecoclimate_data.R index f7fd179..140a39d 100644 --- a/leer_ecoclimate_data.R +++ b/leer_ecoclimate_data.R @@ -88,7 +88,7 @@ map_holocene<- CCSMpres[[1]]$bio.1*0 map_lgmaximum<- CCSMpres[[1]]$bio.1*0 sta<- 0 areas_res<- NULL - +res_maps<- list () x<- 1 for (i in 1:7){ nicho_pres<- as.data.frame (extract (lala[[x]], coord)) @@ -103,7 +103,7 @@ for (i in 1:7){ datos<- datos [complete.cases(datos), ] # set percentile to be drop out - x1<- 0.02 + x1<- 0.01 min_t<- quantile (datos$bio.1, x1) max_t<- quantile (datos$bio.1, 1-x1) @@ -135,11 +135,9 @@ for (i in 1:7){ map_pres_precipmin<- reclassify (lala[[x]]$bio.16, c(-Inf, min_pmin, 0, min_pmin, max_pmin, 1, max_pmin, +Inf, 0)) map_ts<- reclassify (lala[[x]]$bio.4, c(-Inf, min_ts, 0, min_ts, max_ts, 1, max_ts, +Inf, 0)) - map1<- map_pres * map_pres_tmin * map_pres_tmax * map_ts * - map_pres_precip * map_pres_precipmax * map_pres_precipmin + map1<- map_pres * map_pres_tmin * map_pres_tmax * map_ts * map_pres_precip * map_pres_precipmax * map_pres_precipmin map1<- mask (map1, costa) - map_pres<- reclassify (lala[[x+1]]$bio.1, c(-Inf, min_t, 0, min_t, max_t, 1, max_t, +Inf, 0)) map_pres_tmin<- reclassify (lala[[x+1]]$bio.11, c(-Inf, min_tmin, 0, min_tmin, max_tmin, 1, max_tmin, +Inf, 0)) map_pres_tmax<- reclassify (lala[[x+1]]$bio.10, c(-Inf, min_tmax, 0, min_tmax, max_tmax, 1, max_tmax, +Inf, 0)) @@ -148,13 +146,11 @@ for (i in 1:7){ map_pres_precipmin<- reclassify (lala[[x+1]]$bio.16, c(-Inf, min_pmin, 0, min_pmin, max_pmin, 1, max_pmin, +Inf, 0)) map_ts<- reclassify (lala[[x]]$bio.4, c(-Inf, min_ts, 0, min_ts, max_ts, 1, max_ts, +Inf, 0)) - map2<- map_pres * map_pres_tmin * map_pres_tmax * map_ts * map_pres_precip * map_pres_precipmax * map_pres_precipmin map2<- mask (map2, costa) - map_pres<- reclassify (lala[[x+2]]$bio.1, c(-Inf, min_t, 0, min_t, max_t, 1, max_t, +Inf, 0)) map_pres_tmin<- reclassify (lala[[x+2]]$bio.11, c(-Inf, min_tmin, 0, min_tmin, max_tmin, 1, max_tmin, +Inf, 0)) map_pres_tmax<- reclassify (lala[[x+2]]$bio.10, c(-Inf, min_tmax, 0, min_tmax, max_tmax, 1, max_tmax, +Inf, 0)) @@ -168,7 +164,8 @@ for (i in 1:7){ map3<- mask (map3, costa) stab<- map1*map2*map3 - + stab2<- map1+map2+map3 + res_maps [[i]]<- stab2 todo_1_0<- extract (stab, coord) matri_1_0<- extract (stab, coord2) big_area_1_0<- extract (stab, araucaria) [[1]] @@ -245,8 +242,17 @@ for (i in 1:7){ } - - +## the other figures +pdf ("stab_ecoclimate.pdf") +par (mfrow =c(2,4)) +for (i in 1:7){ +plot (res_maps [[1]], ylim=c(-45, -10), axes=F, box=F, + legend=F, main=modelos [i]) +plot (araucaria, add=T) +points (coord2) +plot (wrld_simpl, add=T) +} +dev.off() alt2<- resample (alt, map1, method="bilinear") alt_pres<- data.frame (alt= extract (alt2, coord), time= 0) @@ -262,9 +268,6 @@ boxplot (alt_datos$alt ~alt_datos$time, axis (1, labels=c("present", "Mid-Holocene", "LGM"), at=c(1,2,3)) axis(2) -summary (alt_pres) -summary (alt_glm) - pdf("alt_tiempo.pdf") par (mfrow=c(1,3)) plot (alt2, zlim=c(0, 1500), xlim=c(-65, 30), @@ -288,28 +291,30 @@ axis(2) dev.off() - +setwd ("C:\\Users\\sara.varela\\Documents\\CIENCIAS\\mariana_araucaria") # ensemble -pdf("alt_tiempo.pdf") + +pdf("ensemble.pdf") par (mfrow=c(1,3), mar=c(0,0,5,0)) plot (map_lgmaximum, xlim=c(-65, 30), ylim=c(-40, 5), axes=F, box=F, - main="LGM distribution", legend=F, zlim=c(4,7)) + main="LGM distribution", legend=F, zlim=c(2,7)) points (coord_glm, col="#00006090", pch=16) plot (wrld_simpl, add=T) plot (map_holocene, xlim=c(-65, 30), ylim=c(-40, 5), axes=F, box=F, - main="Holocene distribution", legend=F, zlim=c(4,7)) + main="Holocene distribution", legend=F, zlim=c(2,7)) points (coord_hol, col="#00006090", pch=16) plot (wrld_simpl, add=T) plot (map_presente, xlim=c(-65, 30), ylim=c(-40, 5), axes=F, box=F, - main="present distribution", zlim=c(4,7)) + main="present distribution", zlim=c(2,7)) plot (araucaria, add=T) points (coord2) plot (wrld_simpl, add=T) dev.off() + pdf("stability.pdf") palette (colorRampPalette(c("red", "blue"))(7)) plot (wrld_simpl, xlim=c(-60, -30), @@ -317,14 +322,21 @@ plot (wrld_simpl, xlim=c(-60, -30), points (coord, col=sta, pch=16) dev.off() +nombres<- rep (modelos, 3) + library (ggplot2) -a<- data.frame (area=as.vector (t(areas_res)), time=c(23,6,0)) +variable<- nombres [order (rep (modelos, 3))] +length (variable) +dim (a) +a<- data.frame (area=as.vector (t(areas_res)), time=c(23,6,0), variable) + + areas_plot<- ggplot(a, aes(x=time,y=area)) + geom_smooth(color="#00006090") + geom_point (color="#00006090") pdf("area.pdf") plot (areas_plot) dev.off() -sta<- sta + coord3$stable -areas_res<- rbind (areas, areas_res) - +pdf("area_modelos.pdf") +ggplot(a, aes(x=time,y=area)) + geom_line(aes(colour=variable)) +dev.off()