File:Mesa verde winter precipitation index 800 1500 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

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

download_data=0

yeara=800 yearb=1500

year1=1200

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
2. sitename="Chaco Canyon"
3. sitee_lat=36.058333
4. sitee_lon =-107.958889

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

1. 1. mesa verde

sitename="Mesa Verde"

sitee_lon =-108.488611 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) }

ncin1<- nc_open(iname1)

lon <- ncvar_get(ncin1, "lon") lat <- ncvar_get(ncin1, "lat") t <- ncvar_get(ncin1, "time") pdsi0 <- ncvar_get(ncin1, "SPI")

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=200, ncols=200)

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

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

df2<-data.frame(cbind(x,y)) names(df2)<-c("x", "y")

fit1 <- smooth.spline(x, y, all.knots=F, nknots=501)

1. valu2= predict(loess(y ~ x))
2. valu2<-filter(df2, filter = rep(1/3, 3), method = 'convolution', sides = 1)

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

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

print(head(valu2))

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

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

title1=paste0("Summer 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="#3f0000", lwd=5,add=T)

1. Add the smooth curve to the plot

lines(fit1, col = "blue", lwd = 2)

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