File:Chaco canyon summer plus winter precipitation index spi 100 1600 ad 1.svg
![File:Chaco canyon summer plus winter precipitation index spi 100 1600 ad 1.svg](https://upload.wikimedia.org/wikipedia/commons/thumb/7/77/Chaco_canyon_summer_plus_winter_precipitation_index_spi_100_1600_ad_1.svg/800px-Chaco_canyon_summer_plus_winter_precipitation_index_spi_100_1600_ad_1.svg.png?20240115100838)
Original file (SVG file, nominally 1,800 × 720 pixels, file size: 226 KB)
Captions
Captions
Summary
[edit]DescriptionChaco canyon summer plus winter precipitation index spi 100 1600 ad 1.svg |
English: Chaco Canyon summer plus winter precipitation index SPI 100 - 1600 AD |
Date | |
Source | Own work |
Author | Merikanto |
This plot is based on North American Seasonal Precipitation Atlas (NASPA)
David W. Stahle, Edward R. Cook, Dorian J. Burnette, Max C. A. Torbenson, Ian M. Howard, Daniel Griffin, Jose Villanueva Diaz, Benjamin I. Cook, A. Park Williams, Emma Watson, David J. Sauchyn, Neil Pederson, Connie A. Woodhouse, Gregory T. Pederson, David Meko, Bethany Coulthard, Christopher J. Crawford. 2020. Dynamics, Variability, and Change in Seasonal Precipitation Reconstructions for North America. Journal of Climate, 33, 3173-3195. doi: 10.1175/JCLI-D-19-0270.1
https://www.ncei.noaa.gov/access/paleo-search/study/29372
Code to generate this plot, some edit needed is under this page
- 3
-
- north american paleo precipitation atlas naspa data extracting and viewing
- "R" code
- 15.1.2024 0000.0000
-
library(raster)
library(terra)
library(ncdf4)
library(ggplot2)
library(pals)
library(zoo)
library(purrr)
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
}
get_raster_from_data<-function(iname1,variable1, year1)
{
ncin1<- nc_open(iname1)
lon <- ncvar_get(ncin1, "lon")
lat <- ncvar_get(ncin1, "lat")
t <- ncvar_get(ncin1, "time")
pdsi0 <- ncvar_get(ncin1, variable1)
#print(t)
#stop(-1)
nc_close(ncin1)
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=200, ncols=200)
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)))
#quit(-1)
return(r1)
}
get_data_rain<-function(iname1,variable1, sitee_lat, sitee_lon, yeara, yearb )
{
ncin1<- nc_open(iname1)
lon <- ncvar_get(ncin1, "lon")
lat <- ncvar_get(ncin1, "lat")
t <- ncvar_get(ncin1, "time")
pdsi0 <- ncvar_get(ncin1, variable1)
#print(t)
#stop(-1)
nc_close(ncin1)
#dim(pdsi0)
#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])
x=selyears2
y=pdsis2
# print(x)
# print(y)
# print( length(x))
# print( length(y))
#sitee_lat
df1<-data.frame(matrix(c(x, y), ncol = 2))
names(df1)<-c("x", "y")
return(df1)
}
- main program
download_data=0
yeara=100
yearb=1600
year1=1175
- Cahokia
- sitename="Cahokia"
- sitee_lat=38.654722
- sitee_lon=-90.059444
- chichen itza
- sitename="Chichén Itzá"
- sitee_lat=20.684167
- sitee_lon =-88.567778
- copan NOK
- sitename="Copán"
- sitee_lat = 14.838139
- sitee_lon = -89.142222
- chaco
sitename="Chaco Canyon"
sitee_lat=36.058333
sitee_lon =-107.958889
- sitename="Mexico City"
- sitee_lat=19.433333
- sitee_lon=-99.133333
- mesa verde
- sitename="Mesa Verde"
- sitee_lon =-108.488611
- sitee_lat =37.183889
- los angeles
- sitename="Los Angeles"
- sitee_lon <- -118.25
- sitee_lat <- 34.05
- sitee_lon <- -80
- sitee_lat <- 40
iname1<-"NASPA_WARM_SPI.nc"
iname2<-"NASPA_COOL_SPI.nc"
- summer
- winter
if(download_data==1)
{
download.file(url = url1,destfile = iname1)
download.file(url = url2,destfile = iname2)
}
- rs1<-get_raster_from_data(iname1,"SPI", 1175)
- plot(rs1)
- rs2<-get_raster_from_data(iname2,"SPI", 1175)
- plot(rs2)
df1<-get_data_rain(iname1,"SPI", sitee_lat, sitee_lon, yeara, yearb )
df2<-get_data_rain(iname2,"SPI", sitee_lat, sitee_lon, yeara, yearb )
print (head(df1, 10))
x1<-as.vector(df1$x)
y1<-as.vector(df1$y)
x2<-as.vector(df2$x)
y2<-as.vector(df2$y)
- print (x1)
- print(y1)
x<-x1
y<-y1+y2
y_moving3=movingAverage(y, n=3, centered=FALSE)
y_mean12=mean(y)
- print (y)
sdy12<-sd(y-y_mean12)
- print(sdy12)
- quit(-1)
dev1=sdy12*2
dryval1<- y_mean12-dev1
moistval1<- y_mean12+dev1
idf1<-which(y < dryval1, arr.ind = TRUE) %>% as.data.frame()
imf1<-which(y > moistval1, arr.ind = TRUE) %>% as.data.frame()
- print(head(imf1))
- ixd5<-as.vector(idf1$col)
- ixm5<-as.vector(imf1$col)
ixd5<-unlist(idf1)
ixm5<-unlist(imf1)
- print(ixm5)
xd5<-(x[ixd5])
xm5<-(x[ixm5])
yd5<-(y[ixd5])
ym5<-(y[ixm5])
- print("moist")
- print(xm5)
- stop(-1)
dd5<-as.vector(rep(dryval1, length(xd5) ) )
dm5<-as.vector(rep(moistval1, length(xm5) ) )
y_moist1<-y_moving3
y_dry1<-y_moving3
y_moist2<-y_moving3
y_dry2<-y_moving3
y_dry1[y_dry1>0]<-0
y_moist1[y_moist1<0]<-0
y_dry2[y_dry2>y_mean12]<-y_mean12
y_moist2[y_moist2<y_mean12]<-y_mean12
- plot(x1,y1)
- quit(-1)
- print(selyears2)
- print(selitems2)
- print(pdsis2)
- fit1 <- smooth.spline(x, y, all.knots=F, nknots=501)
- valu2= predict(loess(y ~ x))
- valu2<-filter(df2, filter = rep(1/3, 3), method = 'convolution', sides = 1)
- valu2<-rollmean(df2, 3,fill = list(NA, NULL, NA))
- valu2<-rollmean(zoo(y,x), 3,fill = list(NA, NULL, NA))
- print(head(valu2))
- str(valu2)
- stop(-1)
- dev.new(width = 1200, height = 600, unit = "px")
title1=paste0("Summer+Winter precipitation index SPI at ", sitename)
- title2=paste0("Winter precipitation index SPI at ", sitename)
pdf(file = paste0("out.pdf"), width = 20, height = 8, colormodel = "rgb")
plot(x, y, type="l", lwd=2, col="red", lty=1,
main=title1,
xlab="Year AD",
ylab="SPI",
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, # Grid line type
col = "gray", # Grid line color
lwd = 1)
lines(x,y , col="#ff7f7f", lwd=2,add=T)
lines(x, y_moving3, col = "#5f0000", lwd = 5)
# abline(h = 2, col="blue", lwd=2, lty=2)
abline(h = y_mean12, col="green", lwd=2, lty=2)
# abline(h = -2, col="red", lwd=2, lty=2)
polygon(x=c(min(x), x, max(x) ) , c(y_mean12, y_dry2,y_mean12), col="brown")
polygon(x=c(min(x), x, max(x) ) , c(y_mean12, y_moist2,y_mean12), col="#007f00")
points(xd5, dd5, col="#7f0000", bg= "#7f0000", pch=24, cex=2)
points(xm5, dm5, col="#00007f", bg= "#00007f", pch=25, cex=2)
- Add the smooth curve to the plot
- lines(fit1, col = "blue", lwd = 2)
- lines(valu2, col = "blue", lwd = 2 )
- axis(1, xaxp=c(700, 1800, 19), las=2)
- lines(x,y,col = "red",lwd = 4, add=T)
dev.off()
system("pdf2svg out.pdf out.svg")
print(".")
quit("yes")
Licensing
[edit]![w:en:Creative Commons](https://upload.wikimedia.org/wikipedia/commons/thumb/7/79/CC_some_rights_reserved.svg/90px-CC_some_rights_reserved.svg.png)
![attribution](https://upload.wikimedia.org/wikipedia/commons/thumb/1/11/Cc-by_new_white.svg/24px-Cc-by_new_white.svg.png)
![share alike](https://upload.wikimedia.org/wikipedia/commons/thumb/d/df/Cc-sa_white.svg/24px-Cc-sa_white.svg.png)
- 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 | 10:08, 15 January 2024 | ![]() | 1,800 × 720 (226 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 | 1440pt |
---|---|
Height | 576pt |