Her ønsker vi å finne fjellarealet per region og andel fjell per region.

nor <- readRDS('../data/norway_outline.RDS')%>%
  st_as_sf()
reg <- st_read("/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp")
## Reading layer `regNorway_wgs84 - MERGED' from data source 
##   `/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/FINAL/Raw_data/Geografisk_oppdeling/regioner_2010/regNorway_wgs84 - MERGED.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.087524 ymin: 57.75901 xmax: 31.76159 ymax: 71.38488
## Geodetic CRS:  WGS 84
reg <-   st_transform(reg, crs = crs(nor))
reg <- st_intersection(reg, nor)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
myColours <- c("blue","red","green","yellow","brown")
reg$region[reg$region=="Ã\u0098stlandet"] <- "Østlandet"
reg$region[reg$region=="Sørlandet"] <- "Sørlandet"
regionNames <- reg$region
plot(nor$geometry, axes=T)
  plot(reg$geometry, add=T, border = "black", 
       col = scales::alpha(myColours, .2))
  legend("bottomright",   
      legend = regionNames, 
      fill = myColours)

Samme kartet, bare med tmap

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(reg) + 
  tm_polygons(col="region", 
              border.col = "white",
              title = "")+
  tm_layout(title="Regioner i Norge",
            legend.text.size = 1)+
  tm_shape(nor)+
  tm_polygons(alpha = 0,border.col = "black")

file <- "/data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/M/Fjell/Mountain ecosystem Norge 50m.tif"
fjell <- raster::raster(file, proxy=F)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
fjell
## class      : RasterLayer 
## dimensions : 32021, 32776, 1049520296  (nrow, ncol, ncell)
## resolution : 50, 50  (x, y)
## extent     : -147675, 1491125, 6401860, 8002910  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs 
## source     : /data/P-Prosjekter/41201042_okologisk_tilstand_fastlandsnorge_2020_dataanaly/M/Fjell/Mountain ecosystem Norge 50m.tif 
## names      : Mountain_ecosystem_Norge_50m 
## values     : 0, 65535  (min, max)

Vi trenger ikke så god oppløsning for dette. Reduserer fra 50x50m til 1x1km ved å ta gjennomsnittet. Jeg finner ikke noen løsning for å redusere oppløsning til stars-objekter så jeg går tilbake til raster. Mulig tarra-pakken kunne klart det like fint.

#fjell[fjell[]>0] <- 1
#fjell_low <- aggregate(fjell, fact=20)
#saveRDS(fjell_low, "../output/fjell_1km.rds")
# tar noen minutter...

fjell_low <- readRDS("../output/fjell_1km.rds")

fjell_low_star <- st_as_stars(fjell_low)
fjell_low_star[fjell_low_star[]==0] <- NA # for plotting
summary(fjell_low)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (3.81% of all cells)
##         layer
## Min.        0
## 1st Qu.     0
## Median      0
## 3rd Qu.     0
## Max.        1
## NA's        0
# reg$arealFjell <- raster::extract(fjell_low, reg, fun = sum)
#saveRDS(reg, "../data/reg_helperfile.rds")
reg_inn <- readRDS("../data/reg_helperfile.rds")

Lager et datafil som jeg lagrer på disk slik at jeg kan bruken den i andre sammenhenger

(fjellareal <- data.frame(Region = reg_inn$region,
                         Fjellareal = reg_inn$arealFjell))
##       Region Fjellareal
## 1 Nord-Norge   59822.74
## 2 Midt-Norge   18294.64
## 3  Østlandet   19004.18
## 4 Vestlandet   20329.63
## 5  Sørlandet    7085.62
#saveRDS(fjellareal, "../data/fjellareal.rds")

Total antall km2 fjell:

sum(fjellareal$Fjellareal) # i km2
## [1] 124536.8
barplot(fjellareal$Fjellareal, names.arg = fjellareal$Region,
        ylab = "Fjellareal (km2)")

reg_inn$arealFjell_round <- round(reg_inn$arealFjell)
tmap_mode("view")
tm_shape(reg_inn) + 
  tm_polygons(col="region", 
              border.col = "white",
              title = "Regioner i Norge",
              alpha=.5)+
  tm_text("arealFjell_round",
          size = 1.5)+
  tm_layout(title="Fjellareal km2",
            legend.text.size = 1)+
  tm_shape(nor$geometry)+
  tm_polygons(alpha = 0,border.col = "black")+
  tm_shape(fjell_low_star)+
    tm_raster(title = "Andel fjellareal")

```