File:Chaco canyon pdsi 100-1700 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

library(raster) library(terra) library(ncdf4) library(ggplot2) library(pals)

download_data=0

yeara=100 yearb=1700

year1=901

1. sitename="Mesa Verde"

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

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

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

1. 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. sitee_lat=38.654722
3. 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"))

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

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

fit1 <- smooth.spline(x, y)

1. valu2=values <- loess(y ~ x)

1. y2=predict(valu2)
2. dev.new(width = 1200, height = 600, unit = "px")

title1=paste0("Drought index PDSI at ", sitename)

pdf(file = paste0("out.pdf"), width = 16, height = 8, colormodel = "rgb")

plot(x, y, type="l", lwd=2, col="red", 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="#3f0000", lwd=5,add=T)

# lines(fit1, col = "blue", lwd = 2)  
    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)

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

dev.off()

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

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