File:Chaco canyon summer plus winter precipitation index spi 100 1600 ad 1.svg

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. north american paleo precipitation atlas naspa data extracting and viewing
3. "R" code
5. 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)

1. 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')

1. print(pdsi_series)

#quit(-1)

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)

}

1. 2. main program

download_data=0

yeara=100 yearb=1600

year1=1175

1. 1. Cahokia
2. sitename="Cahokia"
3. sitee_lat=38.654722
4. sitee_lon=-90.059444

1. 1. chichen itza
2. sitename="Chichén Itzá"
3. sitee_lat=20.684167
4. sitee_lon =-88.567778

1. 1. copan NOK
2. sitename="Copán"
3. sitee_lat = 14.838139
4. sitee_lon = -89.142222

1. 1. chaco

sitename="Chaco Canyon" sitee_lat=36.058333 sitee_lon =-107.958889

1. sitename="Mexico City"
2. sitee_lat=19.433333
3. sitee_lon=-99.133333

1. 1. mesa verde
2. sitename="Mesa Verde"
3. sitee_lon =-108.488611
4. sitee_lat =37.183889

1. 1. los angeles
2. sitename="Los Angeles"
3. sitee_lon <- -118.25
4. sitee_lat <- 34.05

1. sitee_lon <- -80
2. sitee_lat <- 40

1. https://www.ncei.noaa.gov/access/paleo-search/study/19119

iname1<-"NASPA_WARM_SPI.nc"

iname2<-"NASPA_COOL_SPI.nc"

1. 1. summer

url1<-"https://www.ncei.noaa.gov/pub/data/paleo/treering/reconstructions/northamerica/NASPA/NASPA_WARM_SPI.nc"

1. 1. winter

url2<-"https://www.ncei.noaa.gov/pub/data/paleo/treering/reconstructions/northamerica/NASPA/NASPA_COOL_SPI.nc"

if(download_data==1) { download.file(url = url1,destfile = iname1) download.file(url = url2,destfile = iname2) }

1. rs1<-get_raster_from_data(iname1,"SPI", 1175)
2. plot(rs1)
3. rs2<-get_raster_from_data(iname2,"SPI", 1175)
4. 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)

1. print (x1)

1. print(y1)

x<-x1 y<-y1+y2

y_moving3=movingAverage(y, n=3, centered=FALSE)

y_mean12=mean(y)

1. print (y)

sdy12<-sd(y-y_mean12)

1. print(sdy12)

1. 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()

1. print(head(imf1))

1. ixd5<-as.vector(idf1$col)
2. ixm5<-as.vector(imf1$col)

ixd5<-unlist(idf1) ixm5<-unlist(imf1)

1. print(ixm5)

xd5<-(x[ixd5]) xm5<-(x[ixm5]) yd5<-(y[ixd5]) ym5<-(y[ixm5])

1. print("moist")
2. print(xm5)

1. 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

1. plot(x1,y1)

1. quit(-1)

1. print(selyears2)
2. print(selitems2)
3. print(pdsis2)

1. fit1 <- smooth.spline(x, y, all.knots=F, nknots=501)
2. valu2= predict(loess(y ~ x))
3. valu2<-filter(df2, filter = rep(1/3, 3), method = 'convolution', sides = 1)

1. valu2<-rollmean(df2, 3,fill = list(NA, NULL, NA))
2. valu2<-rollmean(zoo(y,x), 3,fill = list(NA, NULL, NA))

1. print(head(valu2))

1. str(valu2)
2. stop(-1)

1. dev.new(width = 1200, height = 600, unit = "px")

title1=paste0("Summer+Winter precipitation index SPI at ", sitename)

1. 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) 

1. Add the smooth curve to the plot
2. lines(fit1, col = "blue", lwd = 2)
3. lines(valu2, col = "blue", lwd = 2 )

1. axis(1, xaxp=c(700, 1800, 19), las=2)

1. lines(x,y,col = "red",lwd = 4, add=T)

dev.off()

system("pdf2svg out.pdf out.svg")

print(".") quit("yes")