0

I am interested in getting a line graph that shows the proportion of ocean (i.e. marine water bodies, that will include inland seas such as the Black sea) out of the total surface area at each latitude. As well as to quantify the complexity of land masses at each latitude. But am not sure where or how to start as I do not typically work with maps, and google has not been particularly insightful.

For the first part to get the proportion (or percentage) of ocean out of total surface area at each latitude, I am interested in recreating the km of land at latitude graph but as a proportion instead of absolute values: https://gis.stackexchange.com/questions/190753/using-gis-to-determine-average-land-distance-from-equator

That would look very similar to the first plot in: http://www.radicalcartography.net/index.html?histland

Not sure where to start on this, though a possibly way will be to calculate the area by specifying an entire polygon that spans each latitude.

For the second part, I imagine it'd be something along the lines of putting a line through each latitude and count the number of intersections of the line with the world map shape file? Such that at 70°N there would be less intersections than at the equator.

Projection wise, I am not particular, and would be happy with what is most accurate and will go with any recommendations.

I am hoping to be able to do the above in R, but would not mind doing the above in either Python or GIS as well if its more logical to do so.

Thank you for your tips and advice in advance!

  • considering the field and the things you referenced, it might be better posting this on GIS.stackexchange.com – Mark Aug 05 '23 at 04:08
  • Interestingly, 70°N is at the north edge of many parts of Canada and Russia where there are many islands, resulting in the highest number of coasts of any latitude, provided we look at it with a high enough resolution. – Jon Spring Aug 05 '23 at 05:47

1 Answers1

2

Borrowing from https://github.com/gkaramanis/xkcd-map to create a data frame with coordinates for the earth's landmasses:

library(rnaturalearth)
library(terra)
library(tidyverse)

world <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
  rmapshaper::ms_dissolve()
r <- terra::rast(world, nrows = 1000, ncols = 1000) # adj. for higher res. if wanted
rr <- rasterize(terra::vect(world), r)
rr_df <- terra::as.data.frame(rr, xy = TRUE)

We can plot this to confirm what it is: 343k (at this resolution) points on the globe where land is. These are unprojected points, so the appearance is quite distorted (e.g. the south pole point gets as much display width as the equator), but is fine for this purpose.

ggplot(rr_df, aes(x,y)) +
  geom_point(size = 0.1) +
  scale_y_continuous(breaks = scales::breaks_width(10))

enter image description here

At this resolution, the earth's surface is sampled 1000x at each of 1000 latitudes, so to get the share of land we can use:

rr_df %>%
  count(y) %>%
  ggplot(aes(n/1000, y)) +
  geom_area(orientation = "y") +
  scale_y_continuous(breaks = scales::breaks_width(10)) +
  scale_x_continuous(labels = scales::percent_format())

enter image description here

For the 2nd part, we could take the data frame and count how often each point is farther than a standard spacing from its prior neighbor on the left (0.36 degrees in my example since I sampled at 1000 pixels per latitude consisting of 360 degrees), indicating a western-facing beach. I think by necessity every latitude must have an eastern-facing beach for every western-facing one, so you could double the numbers to get the total number of coasts, if you prefer.

rr_df %>%
  arrange(y, x) %>%
  mutate(lat_grp = cumsum(x > lag(x, default = -179.82) + 0.37), .by = y) %>%
  summarize(coasts_per_lat = max(lat_grp), .by = y) %>%
  ggplot(aes(coasts_per_lat, y)) +
  geom_line(orientation = "y") +
  scale_y_continuous(breaks = scales::breaks_width(10)) 

Here I've composed this with the map above using the patchwork package:

enter image description here

Or we could combine both by mapping the # of coasts in the background with the land on top.

enter image description here

Jon Spring
  • 55,165
  • 4
  • 35
  • 53
  • Thanks @JonSpring, this works wonders! For the first part, I was unable to get `rmapshaper::ms_dissolve()` to work, there was this error: `Error: SyntaxError: Unexpected token .` but the code works well otherwise. I am trying to understand the 2nd part of the code, could you briefly elaborate on what the mutate and summarize line does? And why is the default value -179.82 - this I assume is the smallest longitude that the calculations start with? And why 0.37, not 0.36 since it was sampled at every 0.36 degrees? Surprised that there is more complexity in the north – Lilnet Cloud Aug 06 '23 at 01:49
  • 1
    Hmm, not sure. My first thought was maybe base pipe |> vs magrittr pipe %>%? For the second part, I wanted a default number that is within one tick from -180. And the 0.37 because floating point numbers meant that some points (which should all be exactly 0.36 apart) were evaluating as being > 0.36; picking something higher fixed that. See https://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal – Jon Spring Aug 06 '23 at 03:37
  • 1
    The mutate line is counting, for each latitude, the cumulative number of times that a point is more than one grid spacing from the prior point. The summarize line grabs the largest of these. One could shorten by using summarize(coasts = sum ()). – Jon Spring Aug 06 '23 at 03:40