File:Chaco canyon pdsi drought index 500-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 drought atlas pdsi data extracting and viewing
3. "R" code
5. ## 14.01.2024 0000.0007

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

}

download_data=0

yeara=500 yearb=1600

year1=960

1. sitename="Mesa Verde"

1. 1. mesa verde
2. sitee_lon =-108.488611
3. sitee_lat =37.183889

1. sitename="Zuni"
2. sitee_lat = 35.069444
3. sitee_lon = -108.846667

1. sitename="Kewa"
2. sitee_lat =35.514444
3. sitee_lon =-106.363333

1. sitename="Acoma"
2. sitee_lat =34.896389
3. sitee_lon =-107.581944

1. sitename="Hopi reservation"
2. sitee_lat = 35.911667
3. sitee_lon = -110.615556

1. sitename="Taos Pueblo"
2. sitee_lat = 36.43917
3. sitee_lon = -105.54559

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

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

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. los angeles
2. sitename="Los Angeles"
3. sitee_lon <- -118.25
4. sitee_lat <- 34.05

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

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

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

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) }

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

1. print(t)

1. stop(-1)

nc_close(ncin1)

1. dim(pdsi0)

numu1=which(t==year1)

print(numu1)

print (dim(pdsi0))

1. stop(-1)

pdsi1 <- pdsi0[numu1,,]

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

1. r2<-resample(r1, s2)

1. writeRaster(r1 , filename="pdsi.nc", bandorder='BSQ',format="NetCDF", overwrite=TRUE)

1. quit(-1)

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

1. quit(-1)

print( dim(pdsi0))

1. pdsix0=as.matrix(pdsi0)

1. print( dim( t(pdsix0)))

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

1. str(r_brick)

1. quit(-1)

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)

1. quit(-1)

print("pazka")

selyears2=seq(from=yeara, to=yearb, by=1)

1. selitems2<-2005-selyears2

selitems2=selyears2 pdsis2=t(pdsi_series[yeara:yearb])

1. sitee_lat

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

x=selyears2 y=pdsis2

1. dex_lower_minus5<-lapply(y,function(y)which(y < -5))

1. print(dex_lower_minus5)

1. xminus5<-as.vector(x[dex_lower_minus5])
2. yminus5<-as.vector(y[dex_lower_minus5])

1. print(xminus5)

library(purrr)

dryval1<- -5.5 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]

1. print(xm5)

dd5<-as.vector(rep(dryval1, length(xd5) ) ) dm5<-as.vector(rep(moistval1, length(xm5) ) )

1. print(ym5)

1. stop(-1)

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")

1. y_fit1=movingAverage(y, n=10, centered=TRUE)

y_fit1=movingAverage(y, n=5, centered=FALSE)

print(y_fit1)

1. plot(x, y_fit1)

1. print(y-y_fit1)

1. stop(-1)

1. y_fit2= <- smooth.spline(x, y)

1. 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 = 16, height = 8, colormodel = "rgb")

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="#ff7f7f", 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)        

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

 lines(x, y_fit1, col = "#5f0000", lwd = 5)  
 

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

1. polygon(x=c(min(x), x, max(x) ) , c(0, y_dry1,0), col="red")
2. 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=2)
 points(xm5, dm5+2, col="#00007f", bg= "#00007f", pch=25, cex=2) 
 

dev.off()

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

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