-
Notifications
You must be signed in to change notification settings - Fork 1
/
Portuguese_Elv.R
107 lines (83 loc) · 4.41 KB
/
Portuguese_Elv.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
Elv_map <- function(iso3, country){
# Country...
if(country == 'Guinee'){
alt <- raster::raster('//dapadfs/workspace_cluster_13/WFP_ClimateRiskPr/1.Data/shps/guinee/GIN_alt/GIN_alt.gri')
}else{
alt <- raster::getData('alt', country = iso3, path = paste0(root,'/1.Data/shps/',country))
}
ext.files <- ls()
if(sum(ext.files %in% c('adm_lvl2')) < 2){
adm_lvl2 <- raster::shapefile('//dapadfs/workspace_cluster_13/WFP_ClimateRiskPr/1.Data/shps/guinea-bissau/request/GNB_adm2.shp')
adm_lvl2 <- adm_lvl2 %>% sf::st_as_sf()
adm_lvl2 <<- adm_lvl2
}
if(sum(ext.files %in% c('shp', 'shp_sf')) < 2){
# Read GDAM's shp at administrative level 1.
shp <<- raster::shapefile(paste0(root , "/1.Data/shps/", tolower(country), "/",tolower(iso3),"_gadm/",country,"_GADM1.shp"))
shp_sf <- shp %>% sf::st_as_sf() %>%
dplyr::group_by(NAME_0) %>% dplyr::summarise() %>% sf::as_Spatial()
shp_sf <<- shp_sf
}
if(sum(ext.files %in% c('regions_all')) == 0){
# Regions shp
regions_all <- raster::shapefile(paste0(root , "/1.Data/shps/", tolower(country), "/",tolower(iso3),"_regions/",tolower(iso3),"_regions.shp"))
regions_all <- regions_all %>% sf::st_as_sf() %>%
dplyr::group_by(region) %>% dplyr::summarise() %>% sf::as_Spatial()
regions_all <<- regions_all
}
if(sum(ext.files %in% c('glwd1', 'glwd2')) < 2){
glwd1 <- raster::shapefile('//dapadfs/workspace_cluster_13/WFP_ClimateRiskPr/1.Data/shps/GLWD/glwd_1.shp' )
crs(glwd1) <- crs(shp)
glwd2 <- raster::shapefile('//dapadfs/workspace_cluster_13/WFP_ClimateRiskPr/1.Data/shps/GLWD/glwd_2.shp' )
crs(glwd2) <- crs(shp)
if(!(iso3 %in% c('NPL', 'PAK', 'NER')) ){
ext.sp <- raster::crop(glwd1, raster::extent(shp))
glwd1 <- rgeos::gSimplify(ext.sp, tol = 0.05, topologyPreserve = TRUE) %>%
sf::st_as_sf()
ext.sp2 <- raster::crop(glwd2, raster::extent(shp))
glwd2 <- rgeos::gSimplify(ext.sp2, tol = 0.05, topologyPreserve = TRUE) %>%
sf::st_as_sf()
}else{
glwd1 <- rgeos::gSimplify(glwd1, tol = 0.05, topologyPreserve = TRUE) %>%
sf::st_as_sf()
glwd2 <- rgeos::gSimplify(glwd2, tol = 0.05, topologyPreserve = TRUE) %>%
sf::st_as_sf()
}
glwd1 <<- glwd1
glwd2 <<- glwd2
}
# =----
getAltitude <- function(iso3 = 'HTI', country = 'Haiti', Zone = 'all'){
adm_c <- regions_all
alt_c <- alt %>% raster::crop(., adm_c) %>% raster::mask(., adm_c)
alt_c <- alt_c %>% raster::rasterToPoints() %>% as_tibble()
xlims <- sf::st_bbox(adm_c)[c(1, 3)]
ylims <- sf::st_bbox(adm_c)[c(2, 4)]
adm_c <- sf::st_as_sf(adm_c)
test <- tibble::as_tibble(st_centroid(adm_c) %>% st_coordinates()) %>%
dplyr::mutate(name = adm_c$region)
# =----------------------------
pp <- ggplot() +
geom_tile(data = alt_c %>% setNames(c('x', 'y', 'alt')), aes(x = x, y = y, fill = alt)) +
geom_sf(data = adm_c, fill = NA, color = gray(.2)) +
geom_sf(data = adm_lvl2, fill = NA, color = gray(.2)) +
geom_sf(data = zone, fill = NA, color = 'black', size = 0.8) +
geom_sf(data = glwd1, fill = 'lightblue', color = 'lightblue') +
geom_sf(data = glwd2, fill = 'lightblue', color = 'lightblue') +
coord_sf(xlim = round(xlims, 2), ylim = round(ylims, 2)) +
scale_fill_gradientn(colours = terrain.colors(20),
guide = guide_colourbar(barheight = 12 ,
barwidth = 2, label.theme = element_text(size = 35))) +
scale_y_continuous(breaks = round(ylims, 2), n.breaks = 3) +
scale_x_continuous(breaks = round(xlims, 2), n.breaks = 3) +
labs(fill = glue::glue('Elevação (m)'), x = NULL, y = NULL) + # title = county ,
theme_bw() + theme(text = element_text(size=35),
legend.title=element_text(size=35),
legend.spacing = unit(5, units = 'cm'),
legend.spacing.x = unit(1.0, 'cm'), plot.title = element_text(hjust = 0.5))
path <- glue::glue('//dapadfs/workspace_cluster_13/WFP_ClimateRiskPr/7.Results/{country}/results/')
dir.create(path,recursive = TRUE)
ggsave(glue::glue('{path}/Elevation.png'), width = 12, height = 12, dpi = 300)
return('Map saved successfully\n')}
getAltitude(iso3 = iso3, country = country, Zone = 'all')
}