File:Lgm habitat europe ccsm4 1.svg

From Wikimedia Commons, the free media repository
Jump to navigation Jump to search

Original file ‎(SVG file, nominally 1,650 × 1,275 pixels, file size: 757 KB)

Captions

Captions

Habitat of humanin Europe during LGM

Summary

[edit]
Description
English: Human habitat in Europe during Last Glacial Maximum. CCSM4 multimethod.
Date
Source Own work
Author Merikanto
SVG development
InfoField
 
The SVG code is valid.
 
This map was created with Adobe Illustrator by Merikanto.
This map is not a real SVG file and is just a wrapper containing one or more raster graphics without vector coding.
Using SVG as just a wrapper is undesirable.
Note: This template should be used if the SVG file contains only raster graphics. If the SVG file contains both vector and raster graphics, then use {{BadSVG}}.

Source of data:

CCSM4 LGM Worldclim GMT rivers fine data, rasterized with QGIS

Visualization is from mean of seven methods with Panoply.

Sites data is from

"Demographic estimates of hunter–gatherers during the Last Glacial Maximum in Europe against the background of palaeoenvironmental data"

plus additional intuitive points,

https://www.sciencedirect.com/science/article/abs/pii/S1040618216300490

library(raster) library(dismo) library(maptools) library(kernlab) require(rJava) library(ecospat) library(randomForest)

data(wrld_simpl)

kerros1<- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

createstack<-function(rext1) { ptopet0<-raster("d:/data_worldclim/ccsm4_lgm_nc/world_lgm_ccsm4_ptopet_25m.nc") gdd00<-raster("d:/data_worldclim/ccsm4_lgm_base/lgm_ccsm4_2-5arcmin_growingDegDays0.tif") 

   tk70<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn7.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx7.tif"))/2
   tk10<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn1.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx1.tif"))/2
   tka0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi1.tif")
   pra0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi12.tif")    
 
   rivers0<-raster("./data_elev/distance_to_rivers.nc")  
  
   tk1<-crop(tk10, rext1)  
   tk7<-crop(tk70, rext1) 
   tka<-crop(tka0, rext1) 
   pra<-crop(pra0, rext1) 
     
   ptopet1<-crop(ptopet0, rext1)  
   gdd01<-crop(gdd00, rext1)  
  
  rivers1<-resample(rivers0,tk7)
     
   rivers <- focal(rivers1, w=matrix(1, 15, 15), mean)   
     
   biome<-tk7
   biome=(biome*0.0)+1.0
    names(biome)<-"biome"

names(tk1)<-"tk1" names(tk7)<-"tk7" names(tka)<-"tka" names(pra)<-"pra"

names(ptopet1)<-"ptopet" names(gdd01)<-"gdd0"

names(rivers)<-"riverz"

  1. tk1[is.na(tk1[])] <- 0
  2. tk7[is.na(tk7[])] <- 0
  3. tka[is.na(tka[])] <- 0
  1. ptopet1<-crop(ptopet0, rext1)

##gdd01<-crop(gdd00, rext1) #pazka

  1. names(ptopet1)<-"bio2"
  2. names(gdd01)<-"biome"
  3. predictors <- stack(gdd01)

#predictors<-stack(tk1,tk7,tka) ##predictors<-stack(tk7,tk1, tka) 

   predictors<-stack(biome,tk7,rivers)
     
   crs(predictors)<-CRS(kerros1)
   
   return(predictors)

}

read_table1<-function(faili, sar1, sar2) {

sitez0<-read.table(faili,sep=';',header=TRUE) occ<-sitez0[,sar1:sar2] occ<-na.omit(occ) #print (head(occ)) #print (predictors) colnames(occ) <- c("Lon","Lat") #print (head(occ)) #okkam1<-as.SpatialPoints(occ) #ok1<-as.matrix(occ) #okkam1<-SpatialPoints(ok1, "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

  1. okkam1<-SpatialPoints(ok1)

#crs(okkam1)<-CRS(kerros1)

  1. return(okkam1)
   return(occ)

}

writeout<-function(oras, outn, varnamex, varunitx, longnamex) {

crs(oras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(oras, filename=outn, overwrite=TRUE, format="CDF", varname=varnamex, varunit=varunitx, longname=longnamex, xname="lon", yname="lat")

}

      2. MAIN ####
  1. predictors <- stack(list.files(file.path(system.file(package="dismo"), 'ex'), pattern='grd$', full.names=TRUE ))
  2. file <- file.path(system.file(package="dismo"), "ex/bradypus.csv")
  3. bradypus <- read.table(file, header=TRUE, sep=',')
  4. bradypus <- bradypus[,-1]

rext1<-c(-15,70,30,65)

  1. rext1<-c(120,180,40,80)

predictors<-createstack(rext1)

  1. occ00<-read.table("./wanha1/lgmas1.txt",sep=';',header=TRUE)
  1. faili1="./wanha1/lgmas1.txt"
  2. occ<-read_table1(faili1, 2, 3)

faili1<-"./csvs/paikatx3.csv"

  1. faili1="./wanha1/lgmas1.txt"
  1. faili1<-"./csvs/solutre3.txt"

bradypus<-read_table1(faili1,1,2)

presvals <- extract(predictors, bradypus) set.seed(0) print ("p1") backgr <- randomPoints(predictors, 500) absvals <- extract(predictors, backgr)

print ("p1")

pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals))) sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals))) sdmdata[,'biome'] <- as.factor(sdmdata[,'biome'])

pred_nf <- dropLayer(predictors, 'biome')

set.seed(0) group <- kfold(bradypus, 5) pres_train <- bradypus[group != 1, ] pres_test <- bradypus[group == 1, ]

ext<-rext1

set.seed(10) backg <- randomPoints(pred_nf, n=1000, ext=ext, extf = 1.25) colnames(backg) = c('Lon', 'Lat') group <- kfold(backg, 5) backg_train <- backg[group != 1, ] backg_test <- backg[group == 1, ]

train <- rbind(pres_train, backg_train) pb_train <- c(rep(1, nrow(pres_train)), rep(0, nrow(backg_train))) envtrain <- extract(predictors, train) envtrain <- data.frame( cbind(pa=pb_train, envtrain) )

envtrain[,'biome'] = factor(envtrain[,'biome'], levels=1:3) head(envtrain)

testpres <- data.frame( extract(predictors, pres_test) ) testbackg <- data.frame( extract(predictors, backg_test) ) testpres[ ,'biome'] = factor(testpres[ ,'biome'], levels=1:3) testbackg[ ,'biome'] = factor(testbackg[ ,'biome'], levels=1:3)

    1. svm

svm <- ksvm(pa ~ tk7, data=envtrain)

esv <- evaluate(testpres, testbackg, svm) esv

ps <- predict(predictors, svm, ext=ext) par(mfrow=c(1,2)) plot(ps, main='Support Vector Machine') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(esv, 'spec_sens') plot(ps > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25)

writeout(ps, "svm.nc", "svm", "u", "svm")

pr1=ps>tr

plot(pr1)

writeout(ps, "svm_pres.nc", "svm_p", "u", "svm_p")

    1. randomForest 4.6-14
    2. Type rfNews() to see new features/changes/bug fixes.

model <- pa ~ tk7 rf1 <- randomForest(model, data=envtrain)

    1. Warning in randomForest.default(m, y, ...): The response has five or fewer
    2. unique values. Are you sure you want to do regression?

model <- factor(pa) ~ tk7 rf2 <- randomForest(model, data=envtrain) rf3 <- randomForest(envtrain[,1:3], factor(pb_train)) erf <- evaluate(testpres, testbackg, rf1) erf

pr <- predict(predictors, rf1, ext=ext) par(mfrow=c(1,2)) plot(pr, main='Random Forest, regression') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(erf, 'spec_sens') plot(pr > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25)

writeout(pr, "rf.nc", "rf", "u", "rf")

pr1=pr>tr

plot(pr1) writeout(pr1, "rf_pres.nc", "rf_p", "u", "rf_p")

    1. glm

gm2 <- glm(pa ~ tk7, family = gaussian(link = "identity"), data=envtrain) ge2 <- evaluate(testpres, testbackg, gm2)

pg <- predict(predictors, gm2, ext=ext) par(mfrow=c(1,2)) plot(pg, main='GLM/gaussian, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(ge2, 'spec_sens') plot(pg > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25)

writeout(pg, "glm.nc", "glm", "u", "glm") pg1=pg>tr plot(pg1) writeout(pg1, "glm_pres.nc", "glm_p", "u", "glm_p")

    1. maha

mm <- mahal(pred_nf, pres_train) e <- evaluate(pres_test, backg_test, mm, pred_nf) e pm = predict(pred_nf, mm, ext=ext, progress=) par(mfrow=c(1,2)) pm[pm < -10] <- -10 plot(pm, main='Mahalanobis distance') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(e, 'spec_sens') plot(pm > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+')

writeout(pm, "maha.nc", "maha", "u", "maha") pm1=pm>tr plot(pm1) writeout(pm1, "maha_pres.nc", "maha_p", "u", "maha_p")

    1. bioclim

bc <- bioclim(pred_nf, pres_train) e <- evaluate(pres_test, backg_test, bc, pred_nf) tr <- threshold(e, 'spec_sens') pb <- predict(pred_nf, bc, ext=ext, progress=)

writeout(pb, "bio.nc", "bio", "u", "bio") pb1=pb>tr plot(pb1) writeout(pb1, "bio_pres.nc", "bio_p", "u", "bio_p")

    1. domain

dm <- domain(pred_nf, pres_train) e <- evaluate(pres_test, backg_test, dm, pred_nf) e

pd = predict(pred_nf, dm, ext=ext, progress=) par(mfrow=c(1,2)) plot(pd, main='Domain, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(e, 'spec_sens') plot(pd > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+')

writeout(pd, "domain.nc", "domain", "u", "domain") pd1=pd>tr plot(pd1) writeout(pd1, "domain_pres.nc", "domain_p", "u", "domain_p")

models <- stack(pb, pd, pg) names(models) <- c("bioclim", "domain", "glm") plot(models)

m <- mean(models) plot(m, main='average score')

auc <- sapply(list(ge2, erf, esv), function(x) x@auc) w <- (auc-0.5)^2 m2 <- weighted.mean( modelsc("glm", "bioclim", "domain"), w) plot(m2, main='weighted mean of three models')

writeout(m, "m.nc", "m", "u", "m")

  1. pb1=pb>tr
  2. plot(pb1)

writeout(m2, "m2.nc", "m2", "u", "m2")

Sites data

Lon;Lat 0.248959292267969;48.2382232173701 -0.325470462716793;48.0785621679185 6.02317505852104;49.2951557335836 8.34482865158445;50.153867184884 11.8153417546174;49.0212369201497 13.2873180017658;48.9191096574132 7.59088959816695;47.709455035682 4.82046275902167;46.3720630297676 4.08447463544747;46.0656957116877 4.12636013841508;46.6356511470423 3.9288999101391;45.1447743615906 4.10840920857183;45.2880863575831 3.59381588639798;45.1363329962802 3.51004488046271;45.2206904095848 3.46217573421398;45.2712448690432 2.85184411954267;46.643867684342 2.91168055235358;46.758768203571 2.349218083931;46.9061381888109 4.3357876532533;47.4673148103602 11.5640287368116;45.5568703449381 11.1810755668217;45.4562263075943 3.79127611467397;47.515832480265 3.80922704451727;47.6772343903265 3.51004488046271;47.5481526935534 3.42627387452743;48.086557022435 2.74413854048303;48.3337813981293 19.528257943944;50.6724363128828 19.7915382483121;50.2610955825263 19.9351456870582;50.138529172223 20.1744914183019;49.5287383327041 18.2836601414771;49.9386872085095 17.9366088311738;48.6984324336702 17.2903753568159;48.5718949601122 17.4579173686865;49.5132007330877 17.0869314852588;49.4120859706974 16.6920110287068;49.2483057946455 16.7159456018311;48.9348351606992 14.4660957281408;48.4609142575879 16.249221425906;48.8482836003954 15.6508570977969;48.7221228817887 15.6628243843591;48.5243615709371 16.0338102677867;48.5322869030909 15.8662682559162;48.3814934517334 14.2387172834594;45.866056797224 16.0936467005976;46.3555463342765 18.0802162699199;46.9143135621808 19.3008794992626;46.1486662645181 20.1984259914263;46.3472861139455 20.5574445882917;48.1345000573429 18.7743188905265;47.5077493145738 18.9059590427105;47.8542016173879 18.5708750189694;47.830105173461 29.8919281067941;49.9463888222076 28.312246280586;48.2860246562238 27.1514194840543;48.6510180725924 26.7684663140645;48.7142273050695 26.4812514365721;48.5956449128672 27.259125063114;48.1584548029475 27.0317466184325;48.3337813981293 26.8881391796863;48.2382232173702 26.6607607350048;47.4673148103602 26.3376439978259;47.0041603721743 25.9786254009604;47.142717821753 13.0300213406789;43.5358327352401 15.5670860918616;41.7936282292632 15.9261046887271;40.0760469764719 18.2238237086661;40.1126661462058 18.0802162699199;40.2224053996106 28.0848678359045;44.4313429884533 27.1753540571787;44.7041661480803 24.0399249778868;43.2575742860116 22.7235234560468;43.7436878467773 18.2716928549149;45.1785273250628 18.1280854161687;44.8400978428107 20.7309702434434;41.1299871143924 20.1685077750208;39.727182852891 20.9104795418761;39.727182852891 21.7122877415423;39.7823847172648 20.8506431090652;39.395043345843 23.196231275253;37.5684449734027 11.7794398949308;43.6138626668477 10.8938606893293;43.074272911103 11.8273090411796;42.6532337643238 11.923047333677;42.459298177568 12.4017387961643;42.459298177568 -9.13937701576427;39.422782372328 -8.32560152953586;39.7547893151501 -8.68462012640133;37.2261822618582 -1.48031361596744;37.8524682245273 0.195106502738117;38.8286264581541 -0.187846667251724;39.0891863125127 2.92364783891574;42.1228995812507 3.11512442391066;43.4229443505996 5.46071259009844;43.3359636207116 1.34396601270764;41.7668560048073 0.003629917743197;42.2824702401768 -2.22228538282276;43.4403255224348 -5.74066763210443;43.6485101806209 -8.75642384577443;39.7363862357604 -8.56494726077951;40.0668891054192 -3.61049062403593;40.5596503876726 -4.08918208652323;38.0223591605297 -2.19835080969839;38.5109457081537 -5.76460220522879;36.460019467409 -5.71673305898006;36.9206385746047 -4.40033153713998;36.7865764545038 -3.80196720903085;36.9397711067791 1.89446119456804;47.9504754788391 2.18167607206042;48.6035590838006 4.76660996949185;44.8485829507466 4.55119881137257;44.0283142843797 2.5406946689259;43.3881670322776 0.745601684598517;43.3359636207116 -0.618668983490293;43.8387144332808 -0.044239228505532;44.8316114850993 1.77478832894622;44.148653977359 0.673797965225421;44.2173095848955 1.5833117439513;43.7350415743018 -1.83933221283291;43.3881670322776 -3.82590178215522;43.5098002904194 -5.04656501149784;43.5098002904194 1.65511546332439;44.6445954277273 1.55937717082693;45.2206904095848 1.91839576769241;45.0349391292476 1.67905003644876;46.5204886010659 1.94233034081677;47.3052654655651 1.24822772021018;47.3377151911363 6.70531039256542;43.5965314198953 7.56695502504257;43.8732320936943 8.04564648752987;44.148653977359 6.1069460644563;43.8041767925782 4.50332966512384;44.4740559570813 4.21611478763146;44.508203839511 -0.403257825371008;45.506570819052 0.649863392101057;44.6275639745576 -0.427192398495373;44.148653977359 1.24822772021018;43.1267008563638 -1.0255567266045;43.7696191725782 -4.71148098775672;43.5271565017702 -4.30459324464252;43.3707708861284 -2.62917312593696;43.4403255224348 -3.94557464777704;43.3881670322776 -5.81247135147752;43.4229443505996 -2.58130397968823;43.2488581222097 -0.451126971619738;45.0856585589659 0.458386807106135;45.1532144771222 1.15248942771272;44.6105275226309 1.17642400083709;45.0856585589659 0.601994245852326;45.4226383073099 0.937078269593437;46.9306605679274 1.24822772021018;46.6849316551859 0.314779368359945;45.6573344459973 1.32003143958328;45.2375468950702 1.20035857396145;44.8994673529734 -5.52525647398514;36.7098640154208 -5.2141070233684;36.8440604383833 -4.92689214587601;36.8823591061342 -4.06524751339887;37.0353617099913 -4.08918208652324;37.5304930021437 -3.51475233153847;37.0926583926782 -3.25147202717046;37.2261822618582 -1.83933221283292;37.3594701352906 -2.0547433709522;37.3404434707763 -2.07867794407657;37.9846395600665 -1.48031361596744;38.2106659393478 -0.738341849112125;38.3797265811846 -0.427192398495378;38.9031698858952 -0.044239228505537;38.9217935285666 -0.379323252246648;39.0891863125127 -8.08625579829221;38.6419261330417 -7.48789147018308;42.6180175802716 -6.91346171519832;42.9517664724742 -2.15048166344966;43.0392960078135 -1.59998648158927;43.1267008563638 -1.48031361596744;43.5791951797614 -0.355388679122282;43.7523328706596 -0.283584959749187;44.0799184183039 1.03281656209089;43.213980992628 1.41576973208073;43.2314220522468 -3.41901403904101;43.3707708861284 -3.6344251971603;43.405558187158 -3.65835977028467;43.561853946595 -6.02788250959681;43.561853946595 -8.85216213827189;39.4782273161934 -9.35478817388356;38.7726676068199 -9.35478817388356;39.1820110748149 -9.33085360075919;39.4597505740167 -9.21118073513737;38.9962392216343 -9.23511530826173;39.2376471801766 -8.58888183390387;39.6074270921491 -8.68462012640133;39.9018369392061 -8.44527439515768;40.194987219913 -2.31802367532022;43.1965349435417 -4.06524751339887;43.5271565017702 -4.49606982963744;43.4229443505996 -5.02263043837348;43.3533697488819 -5.6209947664826;43.405558187158 0.841339977095973;44.9164188221005 1.67905003644875;46.7505691326905 0.817405403971606;46.504016841466 0.793470830847241;46.156956451219 -0.283584959749187;45.8076932469028 0.482321380230495;45.791006703284 0.4105176608574;45.5736268570226 0.027564490867559;45.0180226543298 1.87052662144367;44.8655494172844 1.70298460957312;45.0349391292476 1.46363887832946;44.9333652920693 1.32003143958327;44.7806670936732 0.147237356489383;44.8655494172844 0.990931059123252;45.6154962703381 1.40978608879964;45.4394348067478 1.63716453348111;45.2880863575831 1.46962252161055;45.2375468950702 1.01486563224762;45.060304468171 1.08666935162071;44.6701332347303 -0.265634029905915;44.8909897437284 -0.205797597095003;44.5934860719867 -0.325470462716828;45.1363329962802 0.524206883198133;45.2038289250683 0.715683468193053;45.2712448690432 0.787487187566148;45.4394348067478 1.03880020537198;45.4730128098952 1.12257121130726;45.3469922011284 0.907160053187974;45.2965052272537 0.141253713208292;45.6322352890653 0.907160053187974;46.6109940485112 0.835356333814879;47.0286378647561 1.06273477849635;46.7833579316585 0.895192766625792;44.508203839511 0.703716181630871;45.0856585589659 0.691748895068688;44.9757095956275 1.0986366381829;45.1532144771222 40.0;54.5 50.0;60.0 60.0;60.0 31.249444;30.064722 23.716667;37.966667 30.523333;50.45 26.255556;47.651389 19.945;50.064722 39.052;51.386 5.37;43.2964 12.5;41.883333 16.366667;48.2 19.051389;47.4925

Additional habitat scripts.

library(raster) library(sp) library(dismo) library(maptools) library(rJava) library(ENMeval) library(viridis) library(RColorBrewer)

rext1<-c(-15,70,30,70) kerros1<- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

create_rstack1<-function(rext1) {

  1. bio1<-raster("d:/data_worldclim/ccsm4_lgm_base/cclgmbi1.tif")
  2. bio10<-raster("d:/data_worldclim/ccsm4_lgm_base/cclgmbi10.tif")
    1. europe

smallxsiz1=275 smallysiz1=200

    1. france
  2. rext1<-c(-5,10,40,50)
    1. eurasia
    1. rext1<-c(-15,180,30,80)
  2. all
  3. rext1<-c(-180,180,-90,90)
   tk70<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn7.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx7.tif"))/2
   tk10<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn1.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx1.tif"))/2

ptopet0<-raster("d:/data_worldclim/ccsm4_lgm_nc/world_lgm_ccsm4_ptopet_25m.nc") annprecip0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi12.tif") anntemp0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi1.tif") warmprecip0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi18.tif") warmtemp0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi10.tif") coldprecip0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi19.tif") coldtemp0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi11.tif")

topowet0<-raster("d:/data_worldclim/ccsm4_lgm_base/lgm_2-5arcmin_topoWet.tif") gdd00<-raster("d:/data_worldclim/ccsm4_lgm_base/lgm_ccsm4_2-5arcmin_growingDegDays0.tif") 


    rivers0<-raster("./data_elev/distance_to_rivers.nc")  
 	 
 	 ## bio 18 warmest precip
 	 ## bio 10 warmest temp
  1. small1<- raster(xy2)
  2. extent(small1) <- rext1
  3. crs(small1) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
 	 # 		xy2 <- matrix(rnorm(smallxsiz1*smallysiz1 ), smallxsiz1,smallysiz1)  
 	  	
 	ptopet1<-crop(ptopet0, rext1) 	
 	annprecip1<-crop(annprecip0, rext1)
   anntemp1<-crop(anntemp0, rext1) 
   warmprecip1<-crop(warmprecip0, rext1)
   warmtemp1<-crop(warmtemp0, rext1)  	
   coldprecip1<-crop(coldprecip0, rext1)
   coldtemp1<-crop(coldtemp0, rext1)  	 
 	topowet1<-crop(topowet0, rext1)	
 	gdd01<-crop(gdd00, rext1)
    tk1<-crop(tk10, rext1)  
   tk7<-crop(tk70, rext1) 

     rivers<-resample(rivers0,anntemp1)

   names(ptopet1)<-"PTOPET"
   names(annprecip1)<-"PrecipAnn"    
   names(anntemp1)<-"TempAnn"
   names(warmprecip1)<-"PrecipWarm"    
   names(warmtemp1)<-"TempWarm"
   names(coldprecip1)<-"PrecipCold"    
   names(coldtemp1)<-"TempCold"
     
   names(topowet1)<-"Topowet" 
   names(gdd01)<-"GDD0"

names(tk1)<-"tk1" names(tk7)<-"tk7" 

  elev0<-raster("./data_elev/area_glacial.nc")
  elev1<-resample(elev0, anntemp1)
  names(elev1)<-"Elevation"

names(rivers)<-"riverz"

#slope0 <- terrain(rout3, opt='slope') #aspect0 <- terrain(rout3, opt='aspect') #hill0 <- hillShade(slope, aspect, hillshade_slope, hillshade_aspect) #plot(hill0, col=grey(0:100/100), legend=FALSE, main='France')

#blurelev1 <- focal(elev1, w=matrix(1,nrow=7,ncol=7), fun=mean, pad=TRUE, na.rm = TRUE) 

   # elev2=elev1
    # elev21[elev2 <0] <- nan

# distrasa1 <- distance(elev2) #tpi5 <- tpiw(elev1, w=15) # rasa1[rasa1 <0] <- 1 


   # names(anntemp1)<-"Elevation"

    ##dstak1<-stack(anntemp1, annprecip1, ptopet1,topowet1)
 ##    dstak1<-stack(coldtemp1, warmtemp1,ptopet1,gdd01)
    # dstak1<-stack(gdd01)
     
     dstak1<-stack(tk7, topowet1, elev1)
     crs(dstak1)<-CRS(kerros1)   
      
      return(dstak1)

}

plottaus1<-function(zvarnaame1, outfilenaame1) {

  1. elevenaame0<-"./data_processing/area_neworog.nc"

elevenaame1<-"./data_elev/area_glacial.nc" 

   while (dev.cur()>1) dev.off()

grayscale_colors <- gray.colors(100,start = 0.0,end = 1.0,gamma = 2.2,alpha = NULL)

relev1<-raster(elevenaame1)

razteri1<-raster(zvarnaame1)

  1. plot(relev1)


#aage<<-fage


ram1<-relev1

#maintitle1<-paste0("Temperature of July, ",toString(aage), " BP") maintitle1<-"Niche of LGM human" subtitle1<-"CCSM4" xlab1<-"Lon" ylab1<-"Lat"


print (maintitle1)

limx1<- -15 limx2<-45 

   limy1<-35
   limy2<-65	
      
   limitx1<-c(limx1, limx2)
   limity1<-c(limy1, limy2)
   
  # templevs1<-c(-15,-10,-5,0,5,10,15,20,25,30,35)
       templevs1<-c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)
   ltemp1<-length(templevs1)

tempmin1<-templevs1[0] 

   tempmax1<-templevs1[ltemp1-1]

print(tempmin1) print(tempmax1)

x1 <- rasterToContour(razteri1, levels=templevs1)

    1. lev1 <- seq(tempmin1,tempmax1, by=5)
 lev1=seq(-15,35,5)

class(x1) ek1<-rasterToContour(ram1)

svg(filename="out.svg", width=12, height=10, pointsize=20)

plot( razteri1, main=maintitle1, sub=subtitle1, xlab=xlab1, ylab=ylab1, xlim=limitx1, ylim=limity1, breaks=templevs1, col=rev(rainbow(ltemp1)) ) 

#  plot(x1,xlim=c( -15, 50 ), ylim=c( 30, 60 ), alpha=0.5,add=TRUE)

plot(ram1, xlim=limitx1, ylim=limity1 , zlim=c(0,1), breaks=c(-1,0,1), col=viridis(2) , legend=FALSE, add=TRUE)


  1. plot(ek1, add=TRUE, legend=FALSE)

dev.off() }

writeout<-function(oras, outn, varnamex, varunitx, longnamex) {

crs(oras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(oras, filename=outn, overwrite=TRUE, format="CDF", varname=varnamex, varunit=varunitx, longname=longnamex, xname="lon", yname="lat")

}

read_table1<-function(faili, sar1, sar2) {

sitez0<-read.table(faili,sep=';',header=TRUE) occ<-sitez0[,sar1:sar2] occ<-na.omit(occ) #print (head(occ)) #print (predictors) colnames(occ) <- c("Lon","Lat") #print (head(occ)) #okkam1<-as.SpatialPoints(occ) #ok1<-as.matrix(occ) #okkam1<-SpatialPoints(ok1, "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

  1. okkam1<-SpatialPoints(ok1)

#crs(okkam1)<-CRS(kerros1)

  1. return(okkam1)
   return(occ)

}

faili1<-"./csvs/paikatx3.csv"

sitez2<-read_table1(faili1,1,2)

dstak1<-create_rstack1(rext1)

head (sitez2)

#BREIK 

    pef0 <- extract(dstak1, sitez2)

    pef1<-as.data.frame(pef0)
    print(pef1)
    
 #   pef2<-data.frame(sitez2$Lon, sitez2$Lat, pef1$TempCold, pef1$TempWarm, pef1$PTOPET, pef1$GDD0   )
     pef2<-data.frame(sitez2$Lon, sitez2$Lat, pef1$tk7)
  
   # pef2<-c(sitez2, pef1)
    
    names(pef2)<-c("lon", "lat", "tk7")
    
    print (pef2)

occ1 <- data.frame(lon=pef2$lon,lat=pef2$lat) occ1 <- na.omit(occ1) coordinates(occ1) <- c("lon", "lat")

crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") proj4string(occ1) <- crs.geo summary(occ1)

bg <- randomPoints(dstak11, n=30000) bg <- as.data.frame(bg)

  1. proj4string(bg) <- crs.geo

eval2 <- ENMevaluate(occ=occ1, env=dstak1, bg.coords=bg, method='checkerboard2', RMvalues=c(1,2), fc=c('L','LQ','LQP'), algorithm='maxent.jar')

  1. <- ENMevaluate(occs, envs, bg, method='checkerboard2', RMvalues=c(1,2), fc=c('L'), algorithm='maxent.jar')

plot(eval2@predictionswhich(eval2@results$delta.AICc==0), main="Relative human occurrence" ) 

  plot(eval2@predictions'L_2', ylim=c(30,70), xlim=c(-15,70), legend=F, main='L_2 prediction')

 plot(eval2@predictions'LQP_1', ylim=c(30,70), xlim=c(-15,70), legend=F, main='LQP_1 prediction')

  # plot(eval2@predictionswhich(eval2@results$avg.test.AUC==0), main="avg auc"  )
   	
   



   raic<-eval2@predictionswhich(eval2@results$delta.AICc==0) # ok

#  r1<-eval2@predictionswhich(eval2@results$delta.AICc==1)
  
  print(eval2@results)

#   r2<-eval2@predictionswhich(eval2@results$avg.test.AUC==0)



#  xm <- maxent(dstak1, occ1, path='MaxEnt_Default') 
#  r <- predict(xm, dstak1, progress='text') 


 	writeout(raic,"./maxent0.nc","LGM_niche", "", "Niche of LGM human")
   plot(raic)

    m1=maxValue(raic)
    m2=minValue(raic)
    m3=m1-m2
    print(m1,m2,m3)   
    rax1=(raic-m2)/abs(m3)
    
 	writeout(rax1,"./maxent1.nc","LGM_niche", "", "Niche of LGM human")
   plot(rax1)
  1. plottaus1("./maxent0.nc", "pazka.nc")


   rl2=eval2@predictions'L_2';
   rqp2=eval2@predictions'LQP_1';
   
    m1=maxValue(rqp2)
    m2=minValue(rqp2)
    m3=m1-m2
    rqp2s=(rqp2-m2)/abs(m3)
    print(m1,m2,m3)
   
 
 
	writeout(rl2,"./mexnet3.nc","l2", "", "l2")
  	writeout(rqp2,"./mexnet4.nc","qp2", "", "qp2")
  	
  	writeout(rqp2s,"./mexnet4s.nc","qp1s", "p", "qp2s")
  1. plot(r2)

Biomod

library(biomod2)

DataSpecies <- read.csv(system.file("external/species/mammals_table.csv", package="biomod2"), row.names = 1)

head(DataSpecies) myRespName <- 'GuloGulo'

myResp <- as.numeric(DataSpecies[,myRespName])

myRespXY <- DataSpecies[,c("X_WGS84","Y_WGS84")]

myExpl = stack( system.file( "external/bioclim/current/bio3.grd", 

                     package="biomod2"),
                system.file( "external/bioclim/current/bio4.grd",
                             package="biomod2"),
                system.file( "external/bioclim/current/bio7.grd",
                             package="biomod2"),
                system.file( "external/bioclim/current/bio11.grd",
                             package="biomod2"),
                system.file( "external/bioclim/current/bio12.grd",
                             package="biomod2"))

myBiomodData <- BIOMOD_FormatingData(resp.var = myResp, 

                                     expl.var = myExpl,
                                     resp.xy = myRespXY,
                                     resp.name = myRespName)

myBiomodOption <- BIOMOD_ModelingOptions()

myBiomodModelOut <- BIOMOD_Modeling( myBiomodData, 

                                       models = c('SRE','RF'),
                                       models.options = myBiomodOption,
                                       NbRunEval=1,
                                       DataSplit=70,
                                       models.eval.meth = c('TSS'),
                                       do.full.models = FALSE)
    1. environmental variables ( bio3 bio4 bio7 bio11 bio12 )

myBiomodProjection <- BIOMOD_Projection(modeling.output = myBiomodModelOut, 

                                          new.env = myExpl,
                                          proj.name = 'current',
                                          selected.models = 'all',
                                          binary.meth = 'TSS',
                                          compress = FALSE,
                                          build.clamping.mask = FALSE)

myExplFuture = stack(system.file("external/bioclim/future/bio3.grd",package="biomod2"), 

                     system.file("external/bioclim/future/bio4.grd",package="biomod2"),
                     system.file("external/bioclim/future/bio7.grd",package="biomod2"),
                     system.file("external/bioclim/future/bio11.grd",package="biomod2"),
                     system.file("external/bioclim/future/bio12.grd",package="biomod2"))

myBiomodProjectionFuture <- BIOMOD_Projection(modeling.output = myBiomodModelOut, 

                                              new.env = myExplFuture,
                                              proj.name = 'future',
                                              selected.models = 'all',
                                              binary.meth = 'TSS',
                                              compress = FALSE,
                                              build.clamping.mask = TRUE)

myBiomodProjectionFuture plot(myBiomodProjectionFuture)

    1. habitat calculator maha

library(raster) library(SDMTools) library(adehabitatHS) library(viridis) library(RColorBrewer)

  1. library(dismo)
  2. require(rJava)
  3. library(ecospat)

kerros1<- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

  1. rext1<-c(-5.0,5.0,42.0,48.0)

rext1<-c(-15.0,75.0,25.0,70.0)

rivers1<-function(rext1) { sabluna0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi10.tif") sabluna<-crop(sabluna0, rext1) 

   elev0<-raster("./data_elev/distance_to_rivers.nc")
    
   rw1<-resample(elev0,sabluna) 
   #rw1<-crop(tk10, rext1)  

   names(rw1)<-"Rivers"

## remove na rw1[is.na(rw1)] <- 0 wxcoords <- xFromCol(rw1) wycoords <- yFromRow(rw1)

wvalues <- values(rw1)

wNewRow <- rep(wycoords,each=NCOL(rw1)) wNewCol <- rep(wxcoords,NROW(rw1))

wdx1 <- round(wNewCol,digits=8) wdy1 <- round(wNewRow,digits=8) wdz1 <- round(wvalues,digits=0)

warm1<-data.frame(wdx1, wdy1, wdz1)

names(warm1)<-c("Lon", "Lat", "Rivers") #print (head(warm1)) return(warm1)


}

elev1<-function(rext1) { sabluna0<-raster("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi10.tif") sabluna<-crop(sabluna0, rext1) 

   elev0<-raster("./data_elev/area_glacial.nc")
    
   rw1<-resample(elev0,sabluna) 
   #rw1<-crop(tk10, rext1)  

   names(rw1)<-"Elevation"

## remove na rw1[is.na(rw1)] <- 0 wxcoords <- xFromCol(rw1) wycoords <- yFromRow(rw1)

wvalues <- values(rw1)

wNewRow <- rep(wycoords,each=NCOL(rw1)) wNewCol <- rep(wxcoords,NROW(rw1))

wdx1 <- round(wNewCol,digits=8) wdy1 <- round(wNewRow,digits=8) wdz1 <- round(wvalues,digits=0)

warm1<-data.frame(wdx1, wdy1, wdz1)

names(warm1)<-c("Lon", "Lat", "tk1") #print (head(warm1)) return(warm1)


}

wclim_temp_7<-function(rext1) { 

   tk10<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn7.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx7.tif"))/2
    
   rw1<-crop(tk10, rext1)  

   names(rw1)<-"tk1"

## remove na rw1[is.na(rw1)] <- 0 wxcoords <- xFromCol(rw1) wycoords <- yFromRow(rw1)

wvalues <- values(rw1)

wNewRow <- rep(wycoords,each=NCOL(rw1)) wNewCol <- rep(wxcoords,NROW(rw1))

wdx1 <- round(wNewCol,digits=8) wdy1 <- round(wNewRow,digits=8) wdz1 <- round(wvalues,digits=0)

warm1<-data.frame(wdx1, wdy1, wdz1)

names(warm1)<-c("Lon", "Lat", "tk1") #print (head(warm1)) return(warm1)


}

wclim_temp_1<-function(rext1) { 

   tk10<-(raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtn1.tif")+raster("d:/data_worldclim/ccsm4_lgm_base/cclgmtx1.tif"))/2
    
   rw1<-crop(tk10, rext1)  

   names(rw1)<-"tk1"

## remove na rw1[is.na(rw1)] <- 0 wxcoords <- xFromCol(rw1) wycoords <- yFromRow(rw1)

wvalues <- values(rw1)

wNewRow <- rep(wycoords,each=NCOL(rw1)) wNewCol <- rep(wxcoords,NROW(rw1))

wdx1 <- round(wNewCol,digits=8) wdy1 <- round(wNewRow,digits=8) wdz1 <- round(wvalues,digits=0)

warm1<-data.frame(wdx1, wdy1, wdz1)

names(warm1)<-c("Lon", "Lat", "tk1") #print (head(warm1)) return(warm1)


}

plottaus1<-function(zvarnaame1, outfilenaame1) {

elevenaame1<-"./data_elev/area_glacial.nc" 

   while (dev.cur()>1) dev.off()

grayscale_colors <- gray.colors(100,start = 0.0,end = 1.0,gamma = 2.2,alpha = NULL)

relev1<-raster(elevenaame1)

razteri1<-raster(zvarnaame1)

  1. plot(relev1)


#aage<<-fage


ram1<-relev1

#maintitle1<-paste0("Temperature of July, ",toString(aage), " BP") maintitle1<-"Niche of Human" subtitle1<-"CCSM4" xlab1<-"Lon" ylab1<-"Lat"


print (maintitle1)

limx1<- -15 limx2<-45 

   limy1<-35
   limy2<-65	
   
      
   limitx1<-c(limx1, limx2)
   limity1<-c(limy1, limy2)
   
  # templevs1<-c(-15,-10,-5,0,5,10,15,20,25,30,35)
       templevs1<-c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0)
   ltemp1<-length(templevs1)

tempmin1<-templevs1[0] 

   tempmax1<-templevs1[ltemp1-1]

print(tempmin1) print(tempmax1)

x1 <- rasterToContour(razteri1, levels=templevs1)

    1. lev1 <- seq(tempmin1,tempmax1, by=5)
 lev1=seq(-15,35,5)

class(x1) ek1<-rasterToContour(ram1)

svg(filename="out.svg", width=12, height=10, pointsize=20)

plot( razteri1, main=maintitle1, sub=subtitle1, xlab=xlab1, ylab=ylab1, xlim=limitx1, ylim=limity1, breaks=templevs1, col=rev(rainbow(ltemp1)) ) 

#  plot(x1,xlim=c( -15, 50 ), ylim=c( 30, 60 ), alpha=0.5,add=TRUE)

plot(ram1, xlim=limitx1, ylim=limity1 , zlim=c(0,1), breaks=c(-1,0,1), col=viridis(2) , legend=FALSE, add=TRUE)


  1. plot(ek1, add=TRUE, legend=FALSE)

dev.off() }

writeout<-function(oras, outn, varnamex, varunitx, longnamex) {

crs(oras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(oras, filename=outn, overwrite=TRUE, format="CDF", varname=varnamex, varunit=varunitx, longname=longnamex, xname="lon", yname="lat")

}

readat1<-function(faili) { rw0<-raster(faili) rw1<-crop(rw0, rext1) ## remove na rw1[is.na(rw1)] <- 0 wxcoords <- xFromCol(rw1) wycoords <- yFromRow(rw1)

wvalues <- values(rw1)

wNewRow <- rep(wycoords,each=NCOL(rw1)) wNewCol <- rep(wxcoords,NROW(rw1))

wdx1 <- round(wNewCol,digits=8) wdy1 <- round(wNewRow,digits=8) wdz1 <- round(wvalues,digits=0)

warm1<-data.frame(wdx1, wdy1, wdz1)

names(warm1)<-c("Lon", "Lat", "w") #print (head(warm1)) return(warm1) }

read_table1<-function(faili, sar1, sar2) {

sitez0<-read.table(faili,sep=';',header=TRUE) occ<-sitez0[,sar1:sar2] occ<-na.omit(occ) print (head(occ)) #print (predictors) colnames(occ) <- c("Lon","Lat") #print (head(occ)) #okkam1<-as.SpatialPoints(occ) ok1<-as.matrix(occ) #okkam1<-SpatialPoints(ok1, "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") okkam1<-SpatialPoints(ok1) crs(okkam1)<-CRS(kerros1) return(okkam1)

}

    1. main proggis
  1. faili1<-read.table("./wanha1/solutre3.txt"
  2. faili1<-"./wanha1/lgmas1.txt"

faili1<-"./csvs/paikatx3.csv"

  1. points1<-read_table1(faili1,2,3)

points1<-read_table1(faili1,1,2)

w1<-readat1("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi10.tif") w2<-readat1("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi11.tif") precip<-readat1("d:/data_worldclim/ccsm4_lgm_bio/cclgmbi12.tif") ptopet<-readat1("d:/data_worldclim/ccsm4_lgm_nc/world_lgm_ccsm4_ptopet_25m.nc") gdd0<-readat1("d:/data_worldclim/ccsm4_lgm_base/lgm_ccsm4_2-5arcmin_growingDegDays0.tif") elev<-elev1(rext1) rivers<-rivers1(rext1)

tk1<-wclim_temp_1(rext1)

tk7<-wclim_temp_7(rext1)

  1. df1<-data.frame(w1[1], w1[2],w1[3],w2[3],w3[3])
  2. names(df1)<-c("Lon","Lat","TW", "TC","PR")
  1. df1<-data.frame(gdd0, ptopet, precip, elev)
  2. names(df1)<-c("Lon","Lat","gdd0", "Lon2","Lat2", "ptopet", "Lon3","Lat3", "precip", "Lon4","Lat4", "Elev" )
  3. head(df1)
  4. selekt1<-c("Lon","Lat","gdd0", "ptopet", "precip", "Elev")

df1<-data.frame(tk7, tk1,rivers)

names(df1)<-c("Lon","Lat","tk7", "Lon2","Lat2", "tk1","Lon3","Lat3", "Rivers" ) head(df1) selekt1<-c("Lon","Lat","tk7", "tk1","Rivers")

keeps <- selekt1

df1<-df1[ , keeps, drop = FALSE]

  1. DF[keeps]

head(df1)

  1. pazka

map1 <- SpatialPixelsDataFrame( points=df1[,c("Lon","Lat")], 

                                    data        = df1,
                                    proj4string = CRS(kerros1))

head(map1)

  1. habit1 <- mahasuhab(map1, points1, type = "probability")

habit1 <- mahasuhab(map1, points1, type = "distance")

plot(habit1)

oras1<-raster(habit1)

names(oras1)<-"p"

writeout(oras1, "katja.nc", "p", "p", "p")

plottaus1("katja.nc", "mhabitat.svg")

Licensing

[edit]
I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current09:42, 1 December 2019Thumbnail for version as of 09:42, 1 December 20191,650 × 1,275 (757 KB)Merikanto (talk | contribs)User created page with UploadWizard

You cannot overwrite this file.

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.

Retrieved from "https://commons.wikimedia.org/w/index.php?title=File:Lgm_habitat_europe_ccsm4_1.svg&oldid=649716730"