File:Pueblo culture area summer pdsi 1956 1-1.png

1. 1. nada pdsi, chelsa cruts
3. 21.01.2024 v 0000.0000

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

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

inprname1<-"./indata5/prec_07.tif" inprname2<-"./indata5/prec_08.tif" inprname3<-"./indata5/tmean_07.tif"

insmallr1<-raster("./resu1/minipdsi.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)

}

get_base_dem_1<-function(cutext1) { # url1="https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_30s_elev.zip" # download.file(url = url1,destfile = './indata3/elev.zip') # unzip("./indata3/elev.zip")

dir.create("./indata6") url1<-"https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/orog/CHELSA_TraCE21k_dem_0_V1.0.tif" destfile1<-"./indata6/CHELSA_TraCE21k_dem_0_V1.0.tif"

1. system("gdal_translate ")

download.file(url = url1,destfile =destfile1)

r1<-raster(destfile1)

1. plot(r1)

rout1<-raster::crop(r1, cutext1) crs(rout1)<-"lonlat"

writefile1<-"./indata6/dem.tif"

writeRaster(rout1, writefile1, format="GTiff", overwrite=TRUE)

r1<-raster(writefile1) plot(r1)

#sabluna1 <- raster(ncol=1200, nrow=960) #extent(sabluna1)<-extent(r1) ##crs(sabluna1)<-crs(r1) return(r1) }

convert_cruts_rasters_to_celsius_mm_and_create_tmean<-function(sourcepath1, destpath1, year1) {

1. quit("yes")

prefix1="CHELSAcruts" svar1="tmin" suffix1="V.1.0.tif" sytamp1=as.character(year1)

dvar1=svar1 #month1=1

coef1=0.1 ofset1=-273.15

for (month1 in 1:12) { smonth1=str_pad(month1, 2, pad = "0") smonth2=as.character(month1)

filename1<-paste0(prefix1,"_",svar1,"_", smonth2,"_",sytamp1, "_",suffix1 ) filename2<-paste0(prefix1,"_",svar1,"_",sytamp1, "_", smonth1,"_",suffix1 )

sourcename1<-paste0(sourcepath1, filename1) destname1<-paste0(destpath1,dvar1, "_",smonth1,".tif") print(sourcename1) print(destname1) rin1<-raster(sourcename1) #plot(rin1) rout1<-(rin1*coef1)+ofset1 #plot(rout1) raster::writeRaster(rout1,destname1, overwrite=TRUE)

}

prefix1="CHELSAcruts" svar1="tmax" suffix1="V.1.0.tif" sytamp1=as.character(year1)

dvar1=svar1 #month1=1

coef1=0.1 ofset1=-273.15 for (month1 in 1:12) { smonth1=str_pad(month1, 2, pad = "0") smonth2=as.character(month1)

filename1<-paste0(prefix1,"_",svar1,"_", smonth2,"_",sytamp1, "_",suffix1 ) filename2<-paste0(prefix1,"_",svar1,"_",sytamp1, "_", smonth1,"_",suffix1 )

sourcename1<-paste0(sourcepath1, filename1) destname1<-paste0(destpath1,dvar1, "_",smonth1,".tif") print(sourcename1) print(destname1) rin1<-raster(sourcename1) #plot(rin1) rout1<-(rin1*coef1)+ofset1 #plot(rout1) raster::writeRaster(rout1,destname1, overwrite=TRUE)

}

prefix1="CHELSAcruts" svar1="prec" suffix1="V.1.0.tif" sytamp1=as.character(year1)

dvar1=svar1 #month1=1

coef1=0.1 ofset1=0 for (month1 in 1:12) { smonth1=str_pad(month1, 2, pad = "0") smonth2=as.character(month1)

filename1<-paste0(prefix1,"_",svar1,"_", smonth2,"_",sytamp1, "_",suffix1 ) filename2<-paste0(prefix1,"_",svar1,"_",sytamp1, "_", smonth1,"_",suffix1 )

sourcename1<-paste0(sourcepath1, filename1) destname1<-paste0(destpath1,dvar1, "_",smonth1,".tif") print(sourcename1) print(destname1) rin1<-raster(sourcename1) #plot(rin1) rout1<-(rin1*coef1)+ofset1 #plot(rout1) raster::writeRaster(rout1,destname1, overwrite=TRUE)

}

prefix1="CHELSAcruts" suffix1="V.1.0.tif" sytamp1=as.character(year1)

svar1="tmin" svar2="tmax"

for (month1 in 1:12) { smonth1=str_pad(month1, 2, pad = "0") smonth2=as.character(month1)

sname1<-paste0(destpath1,svar1, "_",smonth1,".tif") sname2<-paste0(destpath1,svar2, "_",smonth1,".tif") destname1<- paste0(destpath1,"tmean", "_",smonth1,".tif") #print(sname1) print(destname1) rin1<-raster(sname1) rin2<-raster(sname2) rout1<-(rin1+rin2)/2 raster::writeRaster(rout1,destname1, overwrite=TRUE)

}

}

get_f_chelsa_cruts<-function(down1,destpath1, destpath2,year1, svar1,svar2, num1, ext1) {

1. https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_cruts/prec/CHELSAcruts_prec_1_1980_V.1.0.tif

urlbase1<-"https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_cruts"

prefix1="CHELSAcruts"

suffix1="V.1.0.tif" syear1=as.character(year1)

#print("k") print(svar1) #print(".") #stop(-1) urlbase2<-paste0(urlbase1,"/", svar1, "/")

#print(urlbase2)

#stop(-1)

#month1=1 for (month1 in 1:num1) { #smonth1=str_pad(month1, 2, pad = "0") smonth2=as.character(month1)

filename2<-paste0(prefix1,"_",svar2,"_", smonth2,"_",syear1, "_",suffix1 )

url1<-paste0(urlbase2, filename2) print(url1) #stop(-1)

destname1<-paste0(destpath1, filename2) destname2<-paste0(destpath2, filename2)

print(destname1)

if(down1==1) { download.file(url = url1,destfile =destname1) }

r1<-raster(destname1)

# plot(r1)

rout1<-crop(r1, ext1) crs(rout1)<-"lonlat"

#writefile1<-"./indata13/npp.tif"

writeRaster(rout1, destname2, format="GTiff", overwrite=TRUE)

r1<-raster(destname2) #plot(r1)

}

}

library(terra) library(raster) library(ncdf4)

library(pals)

nada_pdsi_cut<-function(iname1, cutext1, beginyr1, endyr1) { ncin1 <- nc_open(iname1) lon1 <- ncvar_get(ncin1,"lon") nlon1 <- dim(lon1) lat1 <- ncvar_get(ncin1,"lat") nlat1 <- dim(lat1) time1 <- ncvar_get(ncin1,"time") ntime1 <- dim(time1) pdsi1 <- ncvar_get(ncin1,"pdsi") #npdsi1 <- dim(pdsi1)

minlon1<-min(lon1) maxlon1<-max(lon1) minlat1<-min(lat1) maxlat1<-max(lat1)

ext1<-c(minlon1, maxlon1, minlat1, maxlat1) crs1<-"lonlat"

#pdsi2<-pdsi1[beginyr1:endyr1,,] #slice2<-apply(pdsi2,c(2,3),mean)

slice2<-pdsi1[beginyr1,,] r1<-raster(slice2)

extent(r1)<-ext1 crs(r1) <-crs1

r1<-flip(r1)

r1[r1<-10]<- NA

   minir1<-crop(r1, cutext1)
   return(minir1)

}

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3

get_base_dem=1 get_chelsa_cruts=1 convert_and_create_tmeans=1 cut_nada_pdsi=1 downscale_mlr3=1

year1=1956

cutminlon1=-111 cutmaxlon1=-105 cutminlat1=34 cutmaxlat1=38

cutext1<-c(cutminlon1, cutmaxlon1, cutminlat1, cutmaxlat1 )

dir.create("./resu1") dir.create("./indata6") dir.create("./indata4") dir.create("./indata3") dir.create("./indata5")

if(get_base_dem==1) { get_base_dem_1(raster::extent(cutext1)) }

if(get_chelsa_cruts==1) { get_f_chelsa_cruts(1,"./indata3/", "./indata4/",year1, "prec", "prec", 12, cutext1) get_f_chelsa_cruts(1,"./indata3/", "./indata4/",year1, "tmin", "tmin", 12, cutext1) get_f_chelsa_cruts(1,"./indata3/", "./indata4/",year1, "tmax", "tmax", 12, cutext1) }

if(convert_and_create_tmeans==1) { convert_cruts_rasters_to_celsius_mm_and_create_tmean("./indata4/", "./indata5/", year1) }

if(cut_nada_pdsi==1) { r1<-nada_pdsi_cut("./indata1/nada_hd2_cl.nc", cutext1, year1,year1) writeRaster(r1, "./resu1/minipdsi.tif", overwrite=T) }

if(downscale_mlr3==1) { dsr1<-mlr3_downscale_pdsi(cutext1, 256,256) plot(dsr1) writeRaster(dsr1, "./resu1/dspdsi1.tif", overwrite=T) }

sabluna1<-raster(ncol=256, nrow=256)

extent(sabluna1)<-cutext1 crs(sabluna1)<-"lonlat"

plotr1<-resample(r1, sabluna1)

1. plot(cr1, col=rev(parula(256)))
2. contour(cr1, add=T)

plot(plotr1, col=rev(parula(256))) contour(plotr1, add=T)

2. 1. plot pueblo area PDSI
3. 1. linux gmt6 bash shell script
5. uses gmt6, gdal, cdo
6. 1. 11.01.2024 0000.0000

1. GMT set/configuration

gmt set MAP_FRAME_TYPE plain gmt set PS_MEDIA 500x500

1. gdal_translate ./indata13/srtm.tif ./indata13/srtm.nc
2. gdal_translate ./indata13/CHELSA_npp_1981-2010_V.2.1.tif ./indata13/npp.nc

1. cdo remapbil,./indata13/srtm.nc ./indata13/npp.nc ./indata13/npp_srtm.nc

1. exit

1. gmt grdgradient ./indata13/srtm.tif -Ghillshade-grad.nc -A345 -Ne0.6 -V
2. gmt grdhisteq hillshade-grad.nc -Ghillshade-hist.nc -N -V
3. gmt grdmath hillshade-hist.nc 2 DIV = shade.nc

1. gmt makecpt -Cbatlow -Iz -T0/10000/2000 -Z > my.cpt

gmt makecpt -Cbatlow -Iz -D -T-5/5/1 -Z > my.cpt

gmt begin southwest pdf

1. gmt basemap -JM12 -R-111/-105/34/38 -Bafg +w100 -B+t"Pueblo region potential NPP gC/m2/yr "

       gmt basemap -JM12 -R-111/-105/34/38 -Bafg   -B+t"Pueblo region summer PDSI 1956 "
 

#gmt grdimage ./indata13/CHELSA_npp_1981-2010_V.2.1.tif -Cmy.cpt -Ishade.nc #-Ei ## -nc gmt grdimage ./resu1/dspdsi1.tif -Cmy.cpt -nc -Ei # -Ishade.nc #-Ei ## -nc

   gmt coast -LjBL+c0+w200k+f+l+o1/2 -Df  -I1/1p,blue  -N2/0.5p,red,-. -N1/0.25p,red,-. -I2/0.25p,blue -I3/0.15p,blue -I4/0.12p,blue -I5/0.12p,blue  -W0.25p,blue -Slightblue
  gmt grdcontour ./resu1/dspdsi1.tif -A1+ap+f20p,0,black -C0.2 -t50
  
  # gmt psscale -Cmy.cpt -DjRM+o0.6i/0+jLM+w2.88i/0.4i+mc -I0.3 -Bx2+lNPP -By+200

 # gmt psscale -Cmy.cpt -B50+l"NPP"
 gmt psscale   -Cmy.cpt -I1 -By+lPDSI_summer
  
  echo -108.472778  37.166667  | gmt psxy -St0.25 -Wblack -Gred  
  echo -108.472778  37.166667  "Cliff Palace"| gmt pstext -F+f8p,1,white+a0+jTL
  echo 107.9617 36.0607 | gmt psxy -St0.25 -Wblack -Gred  
  echo 107.9617 36.0607  "Chaco Canyon"| gmt pstext -F+f8p,1,white+a0+jTL  
  echo -110.53448 36.68415 | gmt psxy -St0.25 -Wblack -Gred 
  echo -110.53448 36.68415 "Betatakin" | gmt pstext -F+f8p,1,white+a0+jTL

echo -107.998124 36.835837 | gmt psxy -St0.25 -Wblack -Gred echo -107.998124 36.835837 "Pueblo Aztec" | gmt pstext -F+f8p,1,white+a0+jTL

  echo  -106.805947 36.23169 |  gmt psxy -St0.25 -Wblack -Gred 
   echo  -106.805947 36.23169 "Gallina" | gmt pstext -F+f8p,1,white+a0+jTL
   

echo -105.54559 36.43917 | gmt psxy -Sc0.25 -Wblack -Gred

	echo -105.54559 36.43917 "Taos" | gmt pstext -F+f8p,1,white+a0+jTL
  echo -108.846667 35.069444 |  gmt psxy -Sc0.25 -Wblack -Gred 
  echo -108.846667 35.069444 "Zuni" | gmt pstext -F+f8p,1,white+a0+jTL
  echo -110.378333 35.835   |  gmt psxy -Sc0.25 -Wblack -Gred 
  echo -110.378333 35.835  "First Mesa" | gmt pstext -F+f8p,1,white+a0+jTL
  echo -106.7275 35.610556 |  gmt psxy -Sc0.25 -Wblack -Gred 
  echo -106.7275 35.610556   "Jemez" | gmt pstext -F+f8p,1,white+a0+jTL

  echo -107.384444 35.018889 |  gmt psxy -Sc0.25 -Wblack -Gred 
  echo -107.384444 35.018889  "Laguna" |  gmt pstext -F+f8p,1,white+a0+jTL 
  echo -106.089167 35.971111 |   gmt psxy -Sc0.25 -Wblack -Gred 
  echo -106.089167 35.971111 "Santa Clara" |   gmt pstext -F+f8p,1,white+a0+jTL 

gmt end

1. pdf2svg southwest.pdf pueblo_culture_area_summer_pdsi__1.svg

pdftoppm -png -rx 1000 -ry 1000 southwest.pdf pueblo_culture_area_summer_pdsi_1956_1