-
Notifications
You must be signed in to change notification settings - Fork 1
/
PalEON_arctic_project.R
116 lines (77 loc) · 3.77 KB
/
PalEON_arctic_project.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
library(neotoma)
<<<<<<< HEAD
# get arctic sites pollen datasets
# surface sample
arctic_suface_pollen_sites <- get_dataset(datasettype = "pollen surface sample", loc=c(-170, 66, -135, 72))
# pollen
arctic_pollen_sites <- get_dataset(datasettype = "pollen", loc=c(-170, 66, -135, 72))
print(arctic_suface_pollen_sites)
print(arctic_pollen_sites)
# download arctic pollen data
arctic_suface_pollen_data <- get_download(arctic_suface_pollen_sites)
arctic_pollen_data <- get_download(arctic_pollen_sites)
print(arctic_suface_pollen_data)
print(arctic_pollen_data)
#---------------------- read shrub map ----------------------#
library(raster)
# locate the dataset in your computer
# setwd("../../PalEON_Arctic_Project/data")
imported_image <- raster("AllShrubMap_masked.img")
projected_raster <- projectRaster(imported_image, crs = "+init=epsg:4326")
plot(projected_raster)
# lat-lon of pollen sites
xy <- surface_pollen[,2:3]
# extract the data from raster within the 10km radius
z <- extract(projected_raster, xy, buffer=20000)
# take the mean
prcnt <- sapply(z, mean, na.rm=TRUE)
#crop the raster a little bit to make it more managable
e <- extent(130000, 250000, 2000000, 2300000)
imported_image_cr <- crop(imported_image, e)
plot(imported_image_cr)
# according to the internets this is how you get lon-lat-data matrix
# Convert raster to SpatialPointsDataFrame
shrub_data <- rasterToPoints(imported_image_cr, spatial=TRUE)
# reproject sp object to lat-lon
shrub_data <- spTransform(shrub_data, CRS("+init=epsg:4326"))
# Assign coordinates to @data slot, display first 6 rows of data.frame
shrub_data@data <- data.frame(shrub_data@data, long=coordinates(shrub_data)[,1],
lat=coordinates(shrub_data)[,2])
head(shrub_data@data)
shrub_df <- shrub_data@data
##---------------------Map of Alaska---------------------------------------#
library(maps)
Alaska.map <- map("world", c("USA:Alaska"), xlim=c(-180,-135), ylim=c(50,72),interior = FALSE)
northslope.map <- map("world", c("USA:Alaska"), xlim=c(-170,-135), ylim=c(66,72),interior = FALSE)
#---------------------- reformat surface pollen data ----------------------#
surface_pollen_taxons <- pollen_taxons <- list()
surface_pollen_counts <- pollen_counts <- list()
#---------------------- reformat surface pollen data ----------------------#
# site.name - lat -lon age Alnus Betula Ericales Populus Salix Total
# site.id - lat -lon age Alnus Betula Ericales Populus Salix Total
names(arctic_surface_pollen_data[[1]])
head(arctic_surface_pollen_data[[1]]$taxon.list)
head(arctic_surface_pollen_data[[1]]$counts)
head(arctic_pollen_data[[1]]$dataset$site.data)
head(arctic_pollen_data[[1]]$ecological)
head(arctic_pollen_data[[1]]$taxon.list)
surface_pollen_taxons <- pollen_taxons <- list()
surface_pollen_counts <- pollen_counts <- list()
# site.name - lat -lon age Alnus Betula Ericales Populus Salix Total
sites.id<-names(arctic_surface_pollen_data)
site.lat<- rep(NA, 148)
site.long<-rep(NA, 148)
site.total<-rep(NA,148)
shrub.total <- rep(NA,148)
shrub.taxa <- c("Alnus", "Betula", "Salix" )
for(i in 1:length(arctic_surface_pollen_data)){
site.lat[i]<-arctic_surface_pollen_data[[i]]$dataset$site.data$lat
site.long[i]<-arctic_surface_pollen_data[[i]]$dataset$site.data$long
site.total[i]<-sum(arctic_surface_pollen_data[[i]]$count)
shrub.total[i] <- sum(arctic_surface_pollen_data[[i]]$count[which(colnames(arctic_surface_pollen_data[[i]]$counts) %in% shrub.taxa)])
}
surface_pollen <- data.frame(site.id = sites.id, site.long = site.long, site.lat = site.lat, site.total = site.total, shrub.total = shrub.total)
surface_pollen$prop<-(surface_pollen$shrub.total/surface_pollen$site.total)*100
library(fields)
dist.mat <- rdist(surface_pollen[,2:3], shrub_df[,2:3])
print("I love camp PalEON!")