File:Salt lake city summer pdsi at 100 1700 ad 1.svg
Original file (SVG file, nominally 1,620 × 720 pixels, file size: 217 KB)
Captions
Summary
[edit]DescriptionSalt lake city summer pdsi at 100 1700 ad 1.svg |
English: Salt Lake City summer PDSI at 100 - 1700 AD |
Date | |
Source | Own work |
Author | Merikanto |
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=1700
- yeara=950
- yearb=1600
- dryval1<- -5.0
- moistval1<- 5.0
dryval1<- -4.5
moistval1<- 4.5
year1=960
- gila cliff
- sitename="Gila Cliff"
- sitee_lon =-108.272222
- sitee_lat = 33.227222
sitename="Salt Lake City"
sitee_lat = 40.760833
sitee_lon = -111.891111
- 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
- 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 | 14:05, 21 January 2024 | 1,620 × 720 (217 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 |