File:Mesa verde drought index pdsi 900 1500 ad 1.svg
Original file (SVG file, nominally 1,620 × 720 pixels, file size: 114 KB)
Captions
Summary
[edit]DescriptionMesa verde drought index pdsi 900 1500 ad 1.svg |
English: Mesa Verde drought index PDSI between 900 - 1500 AD |
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 to download and plot data
- 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=950
yearb=1600
year1=960
- mesa verde
sitename="Mesa Verde"
sitee_lon = -108.488611
sitee_lat = 37.183889
- gila cliff
- sitename="Gila Cliff"
- sitee_lat =-108.272222
- sitee_lon = 33.227222
- casa grande
- 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="Salt Lake City"
- sitee_lon = 40.760833
- sitee_lat = -111.891111
- pharo vilage (fremont culture)
- sitename="Pharo Village"
- sitee_lon =-112.104722
- sitee_lat =39.2475
- 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)
dryval1<- -5.0
moistval1<- 5.0
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 | 12:52, 21 January 2024 | 1,620 × 720 (114 KB) | Merikanto (talk | contribs) | Correction of flaws | |
14:19, 13 January 2024 | 1,440 × 720 (115 KB) | Merikanto (talk | contribs) | Update | ||
14:05, 13 January 2024 | 1,440 × 720 (96 KB) | Merikanto (talk | contribs) | Update | ||
13:30, 10 January 2024 | 1,440 × 720 (85 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 |