This is what a million miles looks like

I don’t have any more flights planned in 2024, so I pulled up my R script from last year to generate a map of the flights I made this year. Then I wondered what the flights would look like for the past decade, from 2015 through 2024. So I made this map.

That’s a little over a million miles in those ten years. I had the most takeoffs and landings at BKK–206 of them. HND was second with 100, and NRT was third with 95. In the USA, I was in and out of SFO 41 times, IAD 19 times, and DEN 18. In Europe, it was close between FRA 16 times, CPH 15 times, and MUC 14 times.

The route I flew the most was Trang to Bangkok (Don Mueang)—I took that route 44 times in the past decade. And I flew Don Mueang to Trang 39 times. I try to take that route by train when I have time. I’ve also driven it more times than I can recall. Again, I try to take the train when I can.

When I take the train I have time to read, think, look out the window, maybe even do a bit of work on my computer. And that brings me to the rturf section of this post.

I hadn’t opened my flight map script for almost a year. I am using a different computer than I used a year ago, so some of the packages in the script needed to be installed. One of those packages was rgdal, and when trying to install it I got this message:

Warning in install.packages :
packagergdalis not available for this version of R

I installed an archived package from CRAN, but I also got messages such as

> Please note that rgdal will be retired during October 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.

I knew I needed to update the code for my flight maps. I thought maybe I could quickly update my code using ChatGPT. There are a lot of problems that I can get help with through ChatGPT. I tried to make the necessary updates using ChatGPT, and I can report here that for this particular issue, the solutions to my problems were primarily from other sources.

I wanted to use the sf library and the Winkel tripel projection. Claus Wilke’s vignette here was extremely helpful to get me started with the base map.

Then I needed to get my great circle flight routes and the airport locations into the sf format. This was the part where I got a working solution from ChatGPT, which I ended up modifying, but the flow goes something like this, taking a dataframe with lat and lon coordinates for airport cities, for example, and converting to an sf object that can be plotted:

## make the airport locations a sf object
cities_sf <- st_as_sf(cities2, coords = c("lon", "lat"), crs = 4326)

## make great circle routes for mapping
gc_routes <- gcIntermediate(flights_unique[c("longitude.y", "latitude.y")],
                            flights_unique[c("longitude.x", "latitude.x")],
                            n = 360, addStartEnd = TRUE, sp = TRUE, 
                            breakAtDateLine = TRUE)

# Convert to sf object
gc_routes_sf <- st_as_sf(gc_routes)

At that point I had a map and cities and routes. But the standard map I had was centered on the prime meridian, at 0° longitude. A lot of my travel is in Asia, and I wanted to make a map that was centered at a different longitude. I got no help (that worked) from ChatGPT with this problem.

But I found the answers I needed from Stack Overflow, especially in answers to these two questions.

And I found in Position Dodge Doesn’t Dodge geom_sf_label a link to ggsflabel which I used for labeling.

I wanted to share these links because they were so useful to me today.

This section was key. I can get the map, state the projection I want to use, with a center on longitude 150 and latitude 13 degrees north, and then the st_break_antimeridian(lon_0 = 150) statement fixes the map polygons so they break properly based on how I’ve set the projection.

## get the world map
world_sf <- giscoR::gisco_get_countries()

## now try to transform this map to center on Asia
crs_wintri_150 <- "+proj=wintri +y_0=0 +lon_0=150 +lat_0=13 +no_defs +datum=WGS84"

world2 <- world_sf |>
  dplyr::filter(NAME_ENGL != "Antarctica") |>
  st_break_antimeridian(lon_0 = 150) |>
  st_transform(crs = crs_wintri_150)

I love maps, but I’m not very good at making them. Whenever I make this kind of progress with my map-making skills, I’m ecstatic.

Related Posts

Previous