Étiquette : 30DayMapChallenge

Intersecting an area of interest with a layer at opening time

Days 3 and 4 of 30DayMapChallenge : « polygons » and « green » (previously).

The CORINE Landcover dataset is distributed as a geopackage weighting more than 8 Go. To limit the memory used when we only work on a subset, we can clip it at opening time. Here we will map the Cyprus Island :

library(dplyr)
library(ggplot2)
library(stringr)
library(sf)
library(rnaturalearth)
library(glue)

# Using the contour of Cyprus (enlarged) from naturalearth to clip
bb <- ne_countries(scale = "medium", country = "cyprus", returnclass = "sf") %>% 
  st_transform("EPSG:3035") %>% 
  st_buffer(90000) %>% 
  pull(geometry) %>% 
  st_as_text()

# Corine Land Cover 
# download from https://land.copernicus.eu/pan-european/corine-land-cover/clc2018
# (registration required)
# passing the bounding area
cyprus_clc <- read_sf("data/U2018_CLC2018_V2020_20u1.gpkg", query = glue("
  SELECT * 
  FROM U2018_CLC2018_V2020_20u1
  WHERE st_intersects(Shape, st_polygonfromtext('{bb}'))"))

legend_colors <- list("TRUE"  = "mediumaquamarine",
                      "FALSE" = "grey90")

legend_labs <- list("TRUE"  = "forest",
                    "FALSE" = "other")

# exclude sea (code 523)
# and classify on forest (codes 3xx)
cyprus_clc %>% 
  filter(Code_18 != "523") %>% 
  ggplot() +
  geom_sf(aes(fill = str_detect(Code_18, "^3"),
              color = str_detect(Code_18, "^3"))) +
  scale_fill_manual(name= "type",
                    values = legend_colors,
                    labels = legend_labs) +
  scale_color_manual(name = "type",
                     values = legend_colors,
                     labels = legend_labs) +
  labs(title = "Cyprus",
       subtitle = "Landcover",
       caption = glue("data: Copernicus CLC 2018
                      projection LAEA
                      r.iresmi.net {Sys.Date()}")) +
  theme_minimal() +
  theme(legend.position = "bottom",
               plot.caption = element_text(size = 7))

ggsave("cyprus.png", width = 20, height = 12.36, units = "cm", scale = 1.1, type = "cairo")

I shouldn’t have…

library(tidyverse)
library(sf)
library(glue)
library(rnaturalearth)
library(units)

# grams of carbon dioxide-equivalents per passenger kilometer
# https://en.wikipedia.org/wiki/Fuel_economy_in_aircraft
co2_eq <- set_units(88, g/km)

# countries map from Naturalearth
countries <- ne_countries(scale = "small", returnclass = "sf")

# airport code and coordinates to geolocate itineraries
airport <- read_csv("https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat",
                    col_names = c("airport",
                                  "name",
                                  "city",
                                  "country",
                                  "iata",
                                  "icao",
                                  "latitude",
                                  "longitude",
                                  "altitude",
                                  "timezone",
                                  "dst",
                                  "tz",
                                  "type",
                                  "source")) %>% 
  # Add Kai Tak, missig from the airport data
  add_row(iata = "HKGX",
          name = "Kai Tak", 
          city = "Hong Kong",
          latitude = 22.328611,
          longitude = 114.194167)

# itineraries
flight <- read_delim("from-to
LYS-LHR
LHR-LYS
LYS-BOD
LYS-BOD
LYS-BOD
LYS-BOD
BOD-LYS
BOD-LYS
BOD-LYS
LYS-BOD
BOD-LYS
BOD-LGW
LHR-JNB
CPT-JNB
JNB-LHR
LHR-ORY
BOD-ORY
CDG-HKGX
HKGX-PER
SYD-HKGX
HKGX-CDG
ORY-CAY
CAY-BEL
BEL-BSB
BSB-MAO
MAO-VVI
VVI-LPB
LPB-MAO
MAO-BEL
BEL-CAY
CAY-XAU
XAU-CAY
CAY-XAU
XAU-CAY
CAY-XAU
XAU-CAY
CAY-ORY
NCE-MXP
MXP-NCE
CDG-CAY
CAY-MPY
MPY-CAY
CAY-CDG
CDG-HKG
HKG-SYD
SYD-HKG
HKG-SYD
TLN-ORY
CDG-CPH
CPH-ORY
ORY-TLN
CDG-YYZ
YYZ-SFO
SFO-YYZ
YYZ-CDG
ORY-TLN
TLN-ORY
LYS-AMS
AMS-SHJ
SHJ-KTM
KTM-SHJ
SHJ-AMS
AMS-LYS
CDG-AUH
AUH-MCT
MCT-KTM
KTM-PKR
PKR-KTM
KTM-MCT
MCT-AUH
AUH-CDG
GVA-FCO
FCO-GVA
CDG-RUN
RUN-CDG
GVA-KEF
KEF-GVA
CDG-ARN
ARN-KRN
KRN-ARN
ARN-CDG
CDG-RUN
RUN-CDG
CDG-RUN
RUN-CDG
CDG-TLS
", delim = "-")

# geolocate
flight_geo <- flight %>% 
  left_join(airport, by = c("from" = "iata")) %>% 
  left_join(airport, by = c("to" = "iata"), suffix = c("_from", "_to"))

# create lines
flight_lines <- flight_geo %>% 
  mutate(line = glue("LINESTRING ({longitude_from} {latitude_from}, {longitude_to} {latitude_to})")) %>% 
  st_as_sf(wkt = "line", crs = "EPSG:4326")

# create great circles and compute costs
flight_geo_gc <- flight_lines %>% 
  st_segmentize(set_units(100, km)) %>% 
  mutate(distance = set_units(st_length(line), km),
         co2 = set_units(distance * co2_eq, t))

# totals
total_flight <- flight_geo_gc %>% 
  st_drop_geometry() %>% 
  summarise(total_distance = sum(distance, na.rm = TRUE),
            total_co2 = sum(co2, na.rm = TRUE))

# map
ggplot() +
  geom_sf(data = countries, fill = "lightgrey", color = "lightgrey") +
  geom_sf(data = flight_geo_gc, color = "red") + 
  # geom_sf(data = flight_lines, color = "blue") + 
  coord_sf(crs = "+proj=eqearth") +
  # coord_sf(crs = "+proj=robin") +
  # coord_sf(crs = "+proj=fouc") +
  # coord_sf(crs = "+proj=eck1") +
  # coord_sf(crs = "+proj=moll") +
  # coord_sf(crs = "+proj=bonne +lat_1=10") +
  # coord_sf(crs = "+proj=laea") +
  labs(title = "My air travel carbon footprint 1993-2023",
       subtitle = glue("{round(total_flight$total_distance, -2)} km - {round(total_flight$total_co2, 1)} teqCO₂")) +
  theme_minimal()