File:Pueblo area summer spi at 1174 ad 1 1 1 1.png

2. 1. north american paleo precipitation atlas naspa data extracting and viewing
3. "R" code
5. 26.1.2024 0000.0003a1

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

library(stringr) library(mlr3) library(mlr3learners)

mlr3_downscale_spi<-function (ext1, xdim1, ydim1) {

inprname1<- "./indata6/gdd0.tif" inprname2<- "./indata6/bio18.tif" inprname3<- "./indata6/annual_precip.tif"

insmallr1<-raster("./outdata1/spi_summer.tif") accupr10<-raster(inprname1) accupr20<-raster(inprname2) accupr30<-raster(inprname3)

plot(accupr10)

#stop(-1)

covo0<-insmallr1

accupr1 <- crop(x = accupr10, y = ext1) accupr2 <- crop(x = accupr20, y = ext1) accupr3 <- crop(x = accupr30, y = ext1)

#xdim1=200 #ydim1=200

sabluna1<-raster(nrow=xdim1,ncol=ydim1)

extent(sabluna1)<-extent(accupr2)

accupr1[is.na(accupr1)] <- -0 accupr2[is.na(accupr2)] <- -0 accupr3[is.na(accupr3)] <- -0

#covo0[is.na(covo0)] <- -0

covo1<-raster::resample(covo0, sabluna1, method="bilinear")

ak1<-raster::resample(accupr1, sabluna1, method="bilinear") ak2<-raster::resample(accupr2, sabluna1, method="bilinear") ak3<-raster::resample(accupr3, sabluna1, method="bilinear")

ck1<-raster::resample(accupr1, covo0, method="bilinear") ck2<-raster::resample(accupr2, covo0, method="bilinear") ck3<-raster::resample(accupr3, covo0, method="bilinear")

#plot(ak4)

#stop(-1)

adf1<-data.frame(cbind(values(ak1), values(ak2), values(ak3) )) cdf1<-data.frame(cbind(values(covo0),values(ck1), values(ck2), values(ck3) ))

names(adf1)<-c("a","b", "c") names(cdf1)<-c("re","a", "b", "c")

print(head(cdf1))

re<-cdf1

#tsk_peba = tsk("cdf1")

tsk_peba = as_task_regr(cdf1, target = "re", id = "peba")

#lrn_rpart = lrn("regr.rpart") #mlr_learners$get("regr.lm") # lrn_rpart = lrn("regr.rpart")

	#mlr_learners$get("regr.ranger")

#lrn_rpart = lrn("regr.ranger")

mlr_learners$get("regr.svm") lrn_rpart = lrn("regr.svm")

# mlr_learners$get("regr.glmnet") # lrn_rpart = lrn("regr.glmnet")

lrn_rpart$train(tsk_peba)

lrn_rpart$model

splits = partition(tsk_peba) splits

prediction = lrn_rpart$predict_newdata(adf1)

#prediction

#prediction$response

red1<-as.numeric(prediction$response)

#print(red1)

red2<-t(matrix(red1, nrow=xdim1, ncol=ydim1)) #print(red2)

#image(red2)

#red2[red2<0]=0

rout1<-raster(red2) crs(rout1)<-crs(sabluna1) extent(rout1)<-extent(sabluna1)

1. plot(rout1, col=rev(parula(100)))

1. contour(rout1, add=T)

1. plot(covo1)
2. contour(covo1, add=T)

  return(rout1)

}

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) #plot(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)

#print("..") 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

1. 1. in first run must set this or obtain data from elsewhere ...

download_data=1 download_chelsa=1

yeara=1170 yearb=1180

1. year1=1156

year1=1174

outpngname1=paste0("pueblo_area_summer_spi_at_",as.character(year1),"_ad_1_1_1_1.png")

ext1<-extent(c(-111,-104, 34, 39))

1. 1. mesa verde

sitename="Mesa Verde" sitee_lon =-108.488611 sitee_lat =37.183889

1. 1. chaco
2. sitename="Chaco Canyon"
3. sitee_lat=36.058333
4. sitee_lon =-107.958889

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. 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. sitee_lon <- -80
2. sitee_lat <- 40

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

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) { iname1<-"NASPA_WARM_SPI.nc" iname2<-"NASPA_COOL_SPI.nc"

dir.create("./indata1") dir.create("./indata2") iname1=paste0("./indata1/", iname1) iname2=paste0("./indata1/", iname2) download.file(url = url1,destfile = iname1) download.file(url = url2,destfile = iname2) }

if(download_chelsa==1) { dir.create("./indata2") dir.create("./indata3") dir.create("./indata4") dir.create("./indata5") dir.create("./indata6")

url1="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio12_1981-2010_V.2.1.tif" url2="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_gdd0_1981-2010_V.2.1.tif" url3="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/orog/CHELSA_TraCE21k_dem_-1_V1.0.tif" url4="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio18_1981-2010_V.2.1.tif" #https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_30s_elev.zip oname1="annual_precip.tif" oname2="gdd0.tif" oname3="dem.tif" oname4="bio18.tif" oname1=paste0("./indata2/", oname1) oname2=paste0("./indata2/", oname2) oname3=paste0("./indata2/", oname3) oname4=paste0("./indata2/", oname4) download.file(url = url1,destfile = oname1) download.file(url = url2,destfile = oname2) download.file(url = url3,destfile = oname3) download.file(url = url4,destfile = oname4)

r1<-raster(oname4) rout1=crop(r1, ext1) raster::writeRaster(rout1, "./indata6/bio18.tif", overwrite=TRUE)

r1<-raster(oname2) rout1=crop(r1, ext1) raster::writeRaster(rout1, "./indata6/gdd0.tif", overwrite=TRUE) r1<-raster(oname1) rout1=crop(r1, ext1) raster::writeRaster(rout1, "./indata6/annual_precip.tif", overwrite=TRUE)

r2<-raster::raster(oname3)

rout1=crop(r2, ext1) sabluna1<-raster::raster("./indata6/bio18.tif") plot(sabluna1) plot(rout1)

rout2<-raster::resample(rout1,sabluna1)

plot(rout2)

raster::writeRaster(rout2, "./indata6/dem.tif", overwrite=TRUE )

## you need gdal: qgis etc!

system("gdaldem hillshade ./indata6/dem.tif ./indata6/hillshade.tif")

}

1. quit(-1)

print("...")

r0<-get_raster_from_data("./indata1/NASPA_WARM_SPI.nc","SPI", year1)

1. quit(-1)

r1<-crop(r0, ext1) rcoarse1<-r1

dir.create("./outdata1/")

raster::writeRaster(r1, "./outdata1/spi_summer.tif", overwrite=TRUE )

1. raccu1<-mlr3_downscale_spi(ext1, 14*2, 10*2)

1. raccu1<-mlr3_downscale_spi(ext1, 840,600)

1. raccu1<-mlr3_downscale_spi(ext1, 120,120)

dem1<-raster("./indata6/dem.tif") hillshade1<-raster("./indata6/hillshade.tif")

1. plot(hillshade1, col=grey(1:100/100))
2. plot(dem1,col=rainbow(100), alpha=0.4, add=T,legend=F)

1. stop(-1)

raccu1<-mlr3_downscale_spi(ext1, 512,512)

r1<-raccu1

raster::writeRaster(r1, "./outdata1/downscaled_spi_summer.tif", overwrite=TRUE )

pale1=(rev(coolwarm(256))) colo1="#001f00"

hillshade2=hillshade1*0.25+0.75

png( outpngname1, width = 1600, height = 1600,

   units = "px",
   pointsize = 36, bg = "white"
   #, res = NA,
   #restoreConsole = TRUE
   )
   
   
   

plot(hillshade2, col=grey(1:100/100), legend=F, main=paste0("Pueblo area summer SPI at ",as.character(year1)), xlab="lon", ylab="lat" )

1. plot(dem1,col=rainbow(100), alpha=0.4, add=T,legend=F)

plot (r1, col=pale1, cex=1, alpha=0.7, add=T)

contour(r1, add=T, lwd=1) points( c(-107.958889),c( 36.058333) , col = colo1, pch=2, cex = 3, lwd=15) text ( c(-107.958889),c( 36.058333+0.15) , "Chaco Canyon" , cex=2,col=colo1 ) points( c(-108.488611),( 37.183889), col = colo1, pch=2,cex = 3, lwd=15) text ( c(-108.488611),( 37.183889+0.15), "Mesa Verde" , cex=2, col=colo1 )

dev.off()

r1<-rcoarse1

plot (r1, col=pale1, main=paste0("Pueblo area summer SPI",as.character(year1) ), xlab="lon", ylab="lat", cex=1 )

contour(r1, add=T, lwd=0.25) points( c(-107.958889),c( 36.058333) , col = colo1, cex = 1.5, lwd=5) text ( c(-107.958889),c( 36.058333+0.15) , "Chaco Canyon" , cex=2,col=colo1 ) points( c(-108.488611),( 37.183889), col = colo1, cex = 1.5, lwd=5) text ( c(-108.488611),( 37.183889+0.15), "Mesa Verde" , cex=2, col=colo1 )

r1<-raccu1

1. r1<-rcoarse1

1. plot (r1, col=pale1,
2. main="Pueblo area summer SPI",
3. xlab="lon", ylab="lat", cex=1)

plot(hillshade2, col=grey(1:100/100), legend=F, main=paste0("Pueblo area summer SPI at ",as.character(year1)), xlab="lon", ylab="lat" )

1. plot(dem1,col=rainbow(100), alpha=0.4, add=T,legend=F)

plot (r1, col=pale1, cex=1, alpha=0.7, add=T)

contour(r1, add=T, lwd=0.52) points( c(-107.958889),c( 36.058333) , col = colo1, cex = 1.5, lwd=5) text ( c(-107.958889),c( 36.058333+0.15) , "Chaco Canyon" , cex=2,col=colo1 ) points( c(-108.488611),( 37.183889), col = colo1, cex = 1.5, lwd=5) text ( c(-108.488611),( 37.183889+0.15), "Mesa Verde" , cex=2, col=colo1 )

1. plot(rs1)
2. contour(rs1, add=T)

1. rs2<-get_raster_from_data(iname2,"SPI", 1175)
2. plot(rs2)

quit("yes")

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

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=3, lw=52
 points(xm5, dm5, col="#00007f", bg= "#00007f", pch=25, cex=3, lw=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")