File:Pueblo culture area summer pdsi 1051 1-1.png
Original file (5,481 × 5,075 pixels, file size: 10.25 MB, MIME type: image/png)
Captions
Summary
[edit]DescriptionPueblo culture area summer pdsi 1051 1-1.png |
English: Pueblo culture area summer PDSI 1051. PDSI is Palmers Drought Index. |
Date | |
Source | Own work |
Author | Merikanto |
This image is based on chalsa cruos, chelsa trace21ka dem and north american drought atlas version 2.
North American drought atlas NADA
Cook, E.R., Seager, R., Heim, R.R., Vose, R.S., Herweijer, C., and Woodhouse, C. 2010. Megadroughts in North America: Placing IPCC projections of hydroclimatic change in a long-term paleoclimate context. Journal of Quaternary Science, 25(1), 48-61. doi: 10.1002/jqs.1303
NOAA Study Page:
https://www.ncei.noaa.gov/access/paleo-search/study/19119
JSON Metadata:
https://www.ncei.noaa.gov/access/paleo-search/study/search.json?xmlId=16785
DIF Metadata:
http://www1.ncdc.noaa.gov/pub/data/metadata/published/paleo/dif/xml/noaa-recon-19119.xml
ISO Metadata:
http://www1.ncdc.noaa.gov/pub/data/metadata/published/paleo/iso/xml/noaa-recon-19119.xml
DOI:
https://doi.org/10.25921/xyj1-3836
Chelsa cruts v 1.0
Karger, D.N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E., Linder, H.P. & Kessler, M. (2017): Climatologies at high resolution for the earth’s land surface areas. Scientific Data. 4, 170122.
https://chelsa-climate.org/chelsacruts/
Chelsa trace21k
Karger, D.N., Nobis, M.P., Normand, S., Graham, C.H., Zimmermann, N. (2023) CHELSA-TraCE21k – High resolution (1 km) downscaled transient temperature and precipitation data since the Last Glacial Maximum. Climate of the Past. https://doi.org/10.5194/cp-2021-30
https://chelsa-climate.org/chelsa-trace21k/
"R" code
-
- nada pdsi, chelsa cruts
- 21.01.2024 v 0000.0002
library(stringr)
library(raster)
library(terra)
library(ncdf4)
library(pals)
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)
- plot(rout1, col=rev(parula(100)))
- contour(rout1, add=T)
- plot(covo1)
- 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"
- system("gdal_translate ")
download.file(url = url1,destfile =destfile1)
r1<-raster(destfile1)
- 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)
{
- 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)
{
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)
}
}
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)
}
- 3
-
- 3
get_base_dem=1
get_chelsa_cruts=1
convert_and_create_tmeans=1
cut_nada_pdsi=1
downscale_mlr3=1
year2=1151 ## "target year" in NADA
year1=1955 ## accurate map year from CHELSA
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, year2,year2)
writeRaster(r1, "./resu1/minipdsi.tif", overwrite=T)
}
if(downscale_mlr3==1)
{
dsr1<-mlr3_downscale_pdsi(cutext1, 256,256)
writeRaster(dsr1, "./resu1/dspdsi1.tif", overwrite=T)
#plot(dsr1)
plot(dsr1, col=rev(parula(256)) )
contour(dsr1, add=T)
}
- sabluna1<-raster(ncol=256, nrow=256)
- extent(sabluna1)<-cutext1
- crs(sabluna1)<-"lonlat"
- r1<- "./resu1/minipdsi.tif"
- plotr1<-resample(r1, sabluna1)
- plot(cr1, col=rev(parula(256)))
- contour(cr1, add=T)
- plot(plotr1, col=rev(parula(256)) )
- contour(plotr1, add=T)
GMT 6.0 code
import numpy as np
import math
import matplotlib.pyplot as plt
mstar=1
a=0.06
hr=0.05 ## disk aspect ratioo in early disk stahes in inner sys o.05, but later stage s0.02 ... 0.03
- hr=0.025
f=10 # separation of protoplanets hill radius
sigmab=10 ## pebbles? !! kg m-2
sigmad=8 ## g cm-2! at 1 au 8 10 etc
sigmagas=math.pow(a, -0.5)*610
- sigmap=5*1e-3*1e-4*math.pow(a, -1) #https://iopscience.iop.org/article/10.3847/1538-4357/ab9eb2/pdf
stokes=0.01 ## 0.01 ... 0.5 ?
alpha=1e-3
mplanet=1e-3*1/330
rhill1=math.pow((mplanet/mstar),1/3)*a
- snow line
- https://arxiv.org/pdf/1505.03516.pdf
a_snow=2.1*math.pow(mstar,1/3 )
spla1=0
aexp1=1.5*(2-spla1)
miso_planetesimal=3.3e-3*math.pow(f, 1.5)*math.pow((sigmad/10), 1.5)*math.pow(a,3)*math.pow(mstar, -0.5)
sidmada=10
- a=np.linspace(0.01,10,999)
a=np.linspace(0.1,6,999)
sigmad=np.power(a, -1.0)*10.0
s2=np.copy(a)
s2=np.where(s2>a_snow,1.0,0.0 )
sidmad=sigmad*s2*4
plt.plot(a,sigmad)
plt.show()
quit(-1)
- print(s2)
- quit(-1)
- mstar=0.1
miso_pt=3.3e-3*math.pow(f, 1.5)*np.power((sigmad/10), 1.5)*np.power(a,3)*math.pow(mstar, -0.5)
- plt.xlim(-0.2, 6)
plt.plot(a,sigmad)
- plt.plot(a,miso_pt)
plt.show()
quit(-1)
- dead zones alpha=1e-6
sigmagas=1e3 ## OK if 8e-8 mo/yr
q=-1/2 # -3/4
p=-1 # -3/4
hr=0.035 ## 1au
eta=0.01
- print(aexp1)
quit(-1)
- miso1_me=3.3e-3*math.pow(f, 1.5)*math.pow((sigmad), 1.5)*math.pow(a,3)
miso1_me=3.3e-3*math.pow(f, 1.5)*math.pow((sigmad/10), 1.5)*math.pow(a,3)
- https://ui.adsabs.harvard.edu/abs/2003DPS....35.2507W/abstract
- miso1_me_2=2.1e-3*math.pow((sigmad), 1.5)*math.pow(a,3)
- --- > mmsn , sigmap=8 g cm-2 , miso 0.05 me
miso1_me2=0.1*math.pow((sigmad/5), 1.5)*math.pow(a,3)*math.pow(mstar, -1/2)
spla1=1 ## gradient of planetesimals
b=10*rhill1 ## number of hilll radiis between planets
miso1_me3=0.16*math.pow((b/(10*rhill1)), 3/2)*math.pow((sigmad/10),3/2)*math.pow(a, 1.5*(2-spla1))*math.pow(mstar, -0.5)
Licensing
[edit]This file is made available under the Creative Commons CC0 1.0 Universal Public Domain Dedication. | |
The person who associated a work with this deed has dedicated the work to the public domain by waiving all of their rights to the work worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. You can copy, modify, distribute and perform the work, even for commercial purposes, all without asking permission.
http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse |
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 09:09, 22 January 2024 | 5,481 × 5,075 (10.25 MB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Horizontal resolution | 393.7 dpc |
---|---|
Vertical resolution | 393.7 dpc |