Catégories
R

Opening a spatial subset with {sf}

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")
Cyprus landcover
Catégories
R

My air travel carbon footprint

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", 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-2020",
       subtitle = glue("{round(total_flight$total_distance, -2)} km - {round(total_flight$total_co2, 1)} teqCO₂")) +
  theme_minimal()
map of flights
Flygskam (equal-earth projection)

Catégories
R

Extract POIs from a Suunto watch

The Suunto watches (Spartan, Suunto 9,…) can record waypoints (or POIs) but although they can be visualized in the Suunto app (or on the watch), they cannot be exported to be used with other tools. It used to be possible to access them from the Movescount website but it was discontinued a few months ago.

So we have to use some black magic tricks to extract them :

  • Make sure Suunto app has storage permission.
  • Go to settings and tap many times the version number. Logging will start.
  • Go to home , pull to refresh the feed
  • Wait a bit
  • Go again to settings tap many times the version it will stop logging
  • On your Android internal storage there will be a folder called stt

Actually, with the current app version, we get an SQLite database in :

android > data > com.stt.android.suunto > files > stt-dump > amer_app.db

Now it’s just a matter of copying the file (via USB, bluetooth, etc.) to a PC, connecting to the database and finding where the POIs are stored (in a table called… pois) ; we can use QGIS or R for that…

With R, we can select this year POIs (dates are stored as UNIX timestamp) and export them as GeoJSON :

library(RSQLite)
library(dplyr)
library(lubridate)
library(sf)
library(leaflet)

cnx <- dbConnect(SQLite(), "~/amer_app.db")

poi <- dbGetQuery(cnx, "SELECT * FROM pois") %>% 
  filter(creation >= format(ymd("2022-01-01"), "%s")) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326")

poi %>% 
  st_write("~/pois.geojson")

Or we can easily map them :

poi %>% 
  leaflet() %>% 
  addCircleMarkers() %>% 
  addTiles()