density curve of spatial object area in R

Question

I'm trying to view the latitudinal and longitudinal distributions of different shapes in R. I'm trying to find the area of an sf object across bands of latitude, basically flattening the shape sideways to see where most of its mass is distributed. How would I go about this?

I thought about using st_cast() to reproject the multipolygon into multipoints, then making a density curve of those points across latitude, but that only makes points around the perimeter of the shape, not the inside. Maybe I can find the distance between multipoints at the same latitudes, then plot those distances as a pseudo-histogram?

I think maybe it will be easier if I convert to raster and count the number of cells in each latitudinal row? Any thoughts appreciated.

Code to get a small shapefile below to start with:

# try with the country of brazil
brazil <- rnaturalearth::ne_countries(continent = "South America", 
                                  returnclass = "sf") %>%
  filter(name == "Brazil")

brazil %>%
  ggplot() +
  geom_sf() +
  theme_minimal()

Desired output: density plot or histogram showing latitude on the Y-axis and area of the country on the x axis.

Thanks!

Answer 1

@Spacedman's answer is very neat. I've included mine in case it's of interest. It calculates an area by dividing the shape into a grid of squares and determines how much area from the shape is in each square. The resulting data frame is then summarised using dplyr to find the area for each latitude (use of dplyr allows summary by longitude as well if desired).

library(rnaturalearth)
library(ggplot2)
library(sf)
library(dplyr)

brazil <- 
    rnaturalearth::ne_countries(continent = "South America", returnclass = "sf") %>%
    filter(name == "Brazil") 

# make a grid of squares (actually, c(1,20 will work as well))
brazil_grid <- sf::st_make_grid(brazil, n=c(20,20))
# find the centres of each square and make a data frame of their coordinates
brazil_grid_centroids = st_centroid(brazil_grid)
# round the latitude and longitude so that later we can group precisely
brazil_grid_details <- 
    st_coordinates(brazil_grid_centroids) %>% as_tibble() %>% 
    mutate(Longitude=round(X,4), Latitude=round(Y,4)) %>%
    select(-X, -Y)

# intersect each square with the shape
brazil_intersection <- st_intersection(brazil_grid, st_geometry(brazil))
# make a true/false flag to allow selection of the rows where there is an intersection
brazil_intersects <- st_intersects(brazil_grid, st_geometry(brazil), sparse=F)
# find the area of each intersection
brazil_areas <- st_area(brazil_intersection)
# create an area column in the data frame of centroids
brazil_grid_details$area <- 0
# divide by 10e6 to make the area square km
brazil_grid_details$area[brazil_intersects] <- brazil_areas/1000000

# sum the area for each Latitude
summary <- brazil_grid_details %>% group_by(Latitude) %>% summarize(area=sum(area))
ggplot(data=summary) + geom_col(mapping=aes(x=Latitude,y=area))

As the number of Latitude points increases, the width of each bar will reduce which causes the area of each to decrease. The graph below shows 20 Latitude points.

Answer 2

I think this will do the trick:

• Put a horizontal line going through each of the vertex points of the polygon, extending past the bounds of the polygon.

• Crop those lines with the polygon.

• Get the line lengths, and sort by latitude.

• Interpolate linearly with latitude.

I think this will work because the length of cross-section changes linearly except at vertex points, so you only need to find the cross-section length and vertex points and then linear interpolation gets you the rest. Can all be done with sf functions, sf_intersection and sf_length might be all you need except for geometry constructors to make the lines.

12 hours later, here's this:

flatten_EW <- function(polygon){
    ## get all the vertex coords
    bp = st_coordinates(polygon)

    ## get the limits in longitude over the polygon
    xr = range(bp[,"X"])

    ## build a data frame of big lines across the polygon
    xxyy = data.frame(t(xr), bp[,"Y"], bp[,"Y"])
    xylines = st_sfc(apply(xxyy, 1, function(v){st_linestring(matrix(v,2,2))}, simplify=FALSE), crs=st_crs(polygon))

    ## make into a spatial data frame
    xylines = st_as_sf(data.frame(y=bp[,"Y"], geometry=xylines))

    ## crop to the polygon
    xylines = st_intersection(xylines, st_geometry(polygon))

    ## order by increasing latitude
    xylines = xylines[order(xylines$y),]

    ## return data frame with line length
    data.frame(y = xylines$y, L = st_length(xylines))
}

Which you use like this:

flBr = flatten_EW(brazil)
plot(flBr$y, flBr$L, type="l")

to produce this:

Interpreting this as an area is tricky - the line lengths are in metres and the latitude is in degrees. And Brazil is large enough that the ellipsoidal earth shape will probably have an effect. For a planar polygon in cartesian coordinates this will be a shape that has the same area. I think.