Search code examples
rmapsgeospatialspatialgeo

Percentage of ocean out of total area and the complexity of landmasses at each latitude

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!

Solution

  • 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