File:Pharo village summer pdsi 100 2000 ad 1.svg
Original file (SVG file, nominally 1,620 × 720 pixels, file size: 275 KB)
Captions
Summary
[edit]DescriptionPharo village summer pdsi 100 2000 ad 1.svg |
English: Pharo Village summer grought index PDSI at 100 - 2000 AD. Fremont culture. |
Date | |
Source | Own work |
Author | Merikanto |
This image is based on Living Blended Drought Atlas LBDA.
Source of PDSI data
The 'Living Blended Drought Atlas (LBDA)' North American Drought Reconstruction for the last 2000 years
Cook, E.R., Seager, R., Heim, R.R., Vose, R.S., Herweijer, C., and Woodhouse, C. 2010. Megadroughts in North America: Placing IPCC projections of hydroclimatic change in a long-term paleoclimate context. Journal of Quaternary Science, 25(1), 48-61. doi: 10.1002/jqs.1303
NOAA Study Page:
https://www.ncei.noaa.gov/access/paleo-search/study/19119
R code at
https://commons.wikimedia.org/wiki/File:Mesa_verde_drought_index_pdsi_900_1500_ad_1.svg
R code
- 3
-
- north american drought atlas pdsi data extracting and viewing
- "R" code
- ## 21.01.2024 0000.0008
-
library(raster)
library(terra)
library(ncdf4)
library(ggplot2)
library(pals)
library(stats)
movingAverage <- function(x, n=1, centered=FALSE) {
if (centered) {
before <- floor ((n-1)/2)
after <- ceiling((n-1)/2)
} else {
before <- n-1
after <- 0
}
s <- rep(0, length(x))
count <- rep(0, length(x))
new <- x
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- 1
while (i <= before) {
new <- c(rep(NA, i), x[1:(length(x)-i)])
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- i+1
}
i <- 1
while (i <= after) {
new <- c(x[(i+1):length(x)], rep(NA, i))
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- i+1
}
s/count
}
- main program
download_data=0
yeara=100
yearb=2000
- yeara=950
- yearb=1600
- dryval1<- -5.0
- moistval1<- 5.0
dryval1<- -4.5
moistval1<- 4.5
year1=960
- gila cliff
- sitename="Gila Cliff"
- sitee_lat =-108.272222
- sitee_lon = 33.227222
- mexicali
- sitename="Mexicali"
- sitee_lon = -115.467778
- sitee_lat = 32.663333
- pharo village (fremont culture)
sitename="Pharo Village"
sitee_lon =-112.104722
sitee_lat =39.2475
- mesa verde
- sitename="Mesa Verde"
- sitee_lon = -108.488611
- sitee_lat = 37.183889
- sitename="Salt Lake City"
- sitee_lon = 40.760833
- sitee_lat = -111.891111
- casa grande hohokam
- sitename="Casa Grande"
- sitee_lat = 32.997005
- sitee_lon = -111.532069
- casa grandes, paquime (mogollon culture)
- sitename="Paquime"
- sitee_lon = -107.9475
- sitee_lat = 30.366389
- salt lake city (fremont culture)
- sitename="Zuni"
- sitee_lat = 35.069444
- sitee_lon = -108.846667
- sitename="Kewa"
- sitee_lat =35.514444
- sitee_lon =-106.363333
- sitename="Acoma"
- sitee_lat =34.896389
- sitee_lon =-107.581944
- sitename="Hopi reservation"
- sitee_lat = 35.911667
- sitee_lon = -110.615556
- sitename="Taos Pueblo"
- sitee_lat = 36.43917
- sitee_lon = -105.54559
- copan NOK
- sitename="Copán"
- sitee_lat = 14.838139
- sitee_lon = -89.142222
- chichen itza
- sitename="Chichén Itzá"
- sitee_lat=20.684167
- sitee_lon =-88.567778
- chaco
- sitename="Chaco Canyon"
- sitee_lat=36.058333
- sitee_lon =-107.958889
- sitename="Mexico City"
- sitee_lat=19.433333
- sitee_lon=-99.133333
- los angeles
- sitename="Los Angeles"
- sitee_lon <- -118.25
- sitee_lat <- 34.05
- Cahokia
- sitename="Cahokia"
- sitee_lat=38.654722
- sitee_lon=-90.059444
- sitee_lon <- -80
- sitee_lat <- 40
iname1<-"nada_hd2_cl.nc"
url1<-"https://www.ncei.noaa.gov/pub/data/paleo/drought/LBDA2010/nada_hd2_cl.nc"
if(download_data==1)
{
download.file(url = url1,destfile = iname1)
}
- iname1<-"./northdata1/mex/NADAv2-2008.nc"
ncin1<- nc_open(iname1)
lon <- ncvar_get(ncin1, "lon")
lat <- ncvar_get(ncin1, "lat")
t <- ncvar_get(ncin1, "time")
pdsi0 <- ncvar_get(ncin1, "pdsi")
- print(t)
- stop(-1)
nc_close(ncin1)
- dim(pdsi0)
numu1=which(t==year1)
print(numu1)
print (dim(pdsi0))
- stop(-1)
pdsi1 <- pdsi0[numu1,,]
- image(pdsi1)
r1 <- raster(pdsi1, xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
print("s")
r1<-flip(r1)
print("s2")
s2<-raster(nrows=1024, ncols=1024)
crs(s2)<-crs(r1)
extent(s2)<-extent(r1)
- r2<-resample(r1, s2)
- writeRaster(r1 , filename="pdsi.nc", bandorder='BSQ',format="NetCDF", overwrite=TRUE)
- quit(-1)
- r1 <- flip(r1, direction='y')
plot(r1, col=rev(parula(64)))
ext1 <- extent(c(xmin = -96, xmax = -85,
ymin = 13, ymax = 22))
r2 <- crop(x = r1, y = ext1)
plot(r2, col=rev(parula(64)) )
- quit(-1)
print( dim(pdsi0))
- pdsix0=as.matrix(pdsi0)
- print( dim( t(pdsix0)))
- quit(-1)
pdsix0<- aperm(pdsi0, c(3,2,1))
print( dim(pdsix0))
r_brick <- brick(pdsix0, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
- str(r_brick)
- quit(-1)
- print( dim( t(pdsi0)))
r_brick <- flip(t(r_brick), direction='y')
pdsi_series <- extract(r_brick, SpatialPoints(cbind(sitee_lon,sitee_lat)), method='simple')
- print(pdsi_series)
- quit(-1)
print("pazka")
selyears2=seq(from=yeara, to=yearb, by=1)
- selitems2<-2005-selyears2
selitems2=selyears2
pdsis2=t(pdsi_series[yeara:yearb])
- sitee_lat
- print(selyears2)
- print(selitems2)
- print(pdsis2)
x=selyears2
y=pdsis2
- dex_lower_minus5<-lapply(y,function(y)which(y < -5))
- print(dex_lower_minus5)
- xminus5<-as.vector(x[dex_lower_minus5])
- yminus5<-as.vector(y[dex_lower_minus5])
- print(xminus5)
library(purrr)
idf1<-which(y < dryval1, arr.ind = TRUE) %>% as.data.frame()
imf1<-which(y > moistval1, arr.ind = TRUE) %>% as.data.frame()
ixd5<-as.vector(idf1$col)
ixm5<-as.vector(imf1$col)
xd5<-x[ixd5]
xm5<-x[ixm5]
yd5<-x[ixd5]
ym5<-x[ixm5]
- print(xm5)
dd5<-as.vector(rep(dryval1, length(xd5) ) )
dm5<-as.vector(rep(moistval1, length(xm5) ) )
- print(ym5)
- stop(-1)
- myts1 <- ts(y, start=c(min(x), 1), end=c(max(x), 1), frequency=1)
ts1 <- ts(y, start=min(x), frequency=1)
df3<-data.frame(x,y)
names(df3)<-c("x", "y")
- y_fit1=movingAverage(y, n=10, centered=TRUE)
y_fit1=movingAverage(y, n=5, centered=FALSE)
print(y_fit1)
- plot(x, y_fit1)
- print(y-y_fit1)
- stop(-1)
- y_fit2= <- smooth.spline(x, y)
- dev.new(width = 1200, height = 600, unit = "px")
title1=paste0("Drought index PDSI at ", sitename)
y_mean1<-mean(y)
y_moist1<-y_fit1
y_dry1<-y_fit1
y_moist2<-y_fit1
y_dry2<-y_fit1
y_dry1[y_dry1>0]<-0
y_moist1[y_moist1<0]<-0
y_dry2[y_dry2>y_mean1]<-y_mean1
y_moist2[y_moist2<y_mean1]<-y_mean1
pdf(file = paste0("out.pdf"), width = 18, height = 8, colormodel = "rgb")
- png(file = paste0("out.png"), width = 1600, height = 800)
par(mar = c(6, 6, 6, 6))
plot(x, y, type="l", lwd=2, col="#ffffff", lty=1,
main=title1,
xlab="Year AD",
ylab="PDSI",
cex.lab=2, cex.axis=1.5, cex.main=2, cex.sub=1.5
,xaxt="n"
)
axis(1, at = seq(yeara, yearb, by = 50), cex.axis=1.5)
grid(nx = NULL, ny = NULL,
lty = 2,
col = "gray",
lwd = 1)
lines(x,y , col="#af5f5f", lwd=2,add=T)
abline(h = 3, col="blue", lwd=2, lty=2)
abline(h = 0, col="green", lwd=2, lty=2)
abline(h = -3, col="red", lwd=2, lty=2)
- axis(1, xaxp=c(700, 1800, 19), las=2)
lines(x, y_fit1, col = "#5f0000", lwd = 5)
- lines(x,y,col = "red",lwd = 4, add=T)
- polygon(x=c(min(x), x, max(x) ) , c(0, y_dry1,0), col="red")
- polygon(x=c(min(x), x, max(x) ) , c(0, y_moist1,0), col="blue")
polygon(x=c(min(x), x, max(x) ) , c(y_mean1, y_dry2,y_mean1), col="red")
polygon(x=c(min(x), x, max(x) ) , c(y_mean1, y_moist2,y_mean1), col="blue")
points(xd5, dd5+1, col="#7f0000", bg= "#7f0000", pch=24, cex=1.5)
points(xm5, dm5+1, col="#00007f", bg= "#00007f", pch=25, cex=1.5)
dev.off()
system("pdf2svg out.pdf out.svg")
print(".")
quit("yes")
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 13:54, 21 January 2024 | 1,620 × 720 (275 KB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Width | 1296pt |
---|---|
Height | 576pt |