5

I have two polygons, where the smaller is completed contained within the larger. These polygons cover a large portion of North America. Here I create them using a bounding box in 4326.

#create two polygons
b1 <- sf::st_bbox(c(xmin = -170, ymin = 23, xmax = -50, ymax = 80), crs = 4326) %>% 
  sf::st_as_sfc()
b2 <- sf::st_bbox(c(xmin = -160, ymin = 26, xmax = -80, ymax = 60), crs = 4326) %>% 
  sf::st_as_sfc()

#plot polygons (4236)
plot(b1, col = 'red')
plot(b2, col = 'yellow', add = TRUE)

enter image description here

We can also visualize in Google Earth.

enter image description here

When I reproject these polygons (sinusoidal projection), however, the smaller is no longer completely contained within the larger.

Why is this?

I assume this has to do with the properties of the projection, but don't fully understand. I would have thought that a polygon contained within a larger polygon would need to stay contained within that polygon regardless of any transformations that are applied.

Apparently not though?

Note I also checked this in QGIS to make sure this isn't some sort of coding/software issue.

#repoject to MODIS sinusoidal projection
PROJ <- '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs'
b1_tr <- sf::st_transform(b1, crs = PROJ)
b2_tr <- sf::st_transform(b2, crs = PROJ)

#plot reprojected polygons
plot(b1_tr, col = 'red')
plot(b2_tr, col = 'yellow', add = TRUE)

Note how the yellow polygon is no longer completely inside the red polygon.

enter image description here

No such weirdness applies to rasters though (that have the same extent as the polygons) that are reprojected. The yellow raster is completely contained within the red raster. See below for rasters in 4326 and Sinusoidal.

#create rasters with same extent as polygons
r1 <- terra::rast(matrix(1, nrow = 100, ncol = 100))
terra::crs(r1) <- "epsg:4326"
terra::ext(r1) <- c(-170, -50, 23, 80)
r2 <- terra::rast(matrix(2, nrow = 100, ncol = 100))
terra::crs(r2) <- "epsg:4326"
terra::ext(r2) <- c(-160, -80, 26, 60)

#plot rasters (4326)
plot(r1, col = 'red')
plot(r2, col = 'yellow', add = TRUE, axes = FALSE)

enter image description here

#repoject to MODIS sinusoidal projection
r1_tr <- terra::project(r1, PROJ)
r2_tr <- terra::project(r2, PROJ)

#plot reprojected polygons
plot(r1_tr, col = 'red')
plot(r2_tr, col = 'yellow', add = TRUE, axes = FALSE)

enter image description here

My take away from this is not to reproject vector files that cover a large spatial area, which makes me somewhat uneasy!

4
  • The sinusoidal projection is a pseudocylindrical equal-area projection displaying all parallels and the central meridian at true scale. see the deformation en.wikipedia.org/wiki/Sinusoidal_projection
    – Mapperz
    Jun 7, 2023 at 20:07
  • Thanks. I get that it's going to be distorted, but I don't understand why the two polygons would be distorted differently. If they were distorted the same, the yellow would still be within the red, right? Particularly since the rasters are distorted the same. Does this have to do with some area preservation?
    – Caseyy
    Jun 7, 2023 at 20:12
  • 2
    Densify your polygons before reprojecting. If only the corner coordinates are reprojected and connected with straight lines the result is not correct.
    – user30184
    Jun 7, 2023 at 20:36
  • 3
    Ah, right. The straight lines drawn between the vertices are missing the "bulge" that occurs between the vertices because of the transformation. Thank you!
    – Caseyy
    Jun 7, 2023 at 20:40

1 Answer 1

5

As user30184 pointed out, this behavior is due to the fact that only the vertices were transformed. A straight line was then drawn between the transformed vertices when in reality the line should not be straight. By densifying the line (i.e., adding extra vertices), the curvature of the polygon is apparent (as in the raster example) and the yellow polygon is contained within the red polygon. This can be done using the smoothr package in R (probably among other methods).

b1_tr_den <- smoothr::densify(b1, n = 1000) %>%
  sf::st_transform(crs = PROJ)
b2_tr_den <- smoothr::densify(b2, n = 1000) %>%
  sf::st_transform(crs = PROJ)

#plot reprojected polygons after densifying
plot(b1_tr_den, col = 'red')
plot(b2_tr_den, col = 'yellow', add = TRUE)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.