rgisspatialr-sfr-sp

density curve of spatial object area in R

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!

Solution

  • I think this will do the trick:

    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:

    enter image description here

    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.