2

I have a polygon that defines the boundaries of a study area. The area is quadrilateral. However, when I convert it from Lambert Conformal Conic projection to WGS84, one connection between two vertices is now drawn from west to east instead of from east to west, making the polygon into an hourglass shape, and no longer outlining the correct study area. How can I prevent this from happening?

library(terra)

#define projections

LCC <- "+init=EPSG:3347"
WGS84 <- "+init=EPSG:4326"

#create polygon
poly_LCC <- rbind(c(3847903, 1983584 ),  
              c(3847903, 5801864),  
              c(8301883, 5801864),  
              c(8301883, 1983584 ))



poly_LCC <- vect(poly_LCC, "polygons", crs = LCC)
plot(poly_LCC)

#project polygon
poly_WGS <- terra::project(poly_LCC, WGS84)
plot(poly_WGS)
canderson156
  • 1,045
  • 10
  • 24
  • I'm not going to leave an answer since I'm not familiar with the terra package, but in sf I'd just put a convex hull around a simple polygon if I wasn't sure the points were in the right order. https://rdrr.io/cran/terra/man/convhull.html – Brian Sep 20 '22 at 14:56
  • 1
    or perhaps `sf::st_make_valid`. Maybe i'm old school, but it used to be that poly(s) started and ended at the same point...so a quadrilateral would have 5? – Chris Sep 20 '22 at 15:33

1 Answers1

4

Your example data

library(terra)
poly_LCC <- rbind(c(3847903, 1983584 ), c(3847903, 5801864), c(8301883, 5801864), c(8301883, 1983584 ))
poly_LCC <- vect(poly_LCC, "polygons", crs="+init=EPSG:3347")
WGS84 <- "+init=EPSG:4326"

This shows what is going on

d <- densify(poly_LCC, 100000)
p <- terra::project(d, WGS84)
plot(p)

# add points to explain what you were seeing:
points(terra::project(poly_LCC, WGS84))

enter image description here

As you can see, some of the study area is in the Eastern hemisphere. I am guessing because the area includes all of North America (and the Aleutian Islands in Alaska cross the -180 to 180 longitude-line).

Also note that, normally, if you want to meaningfully project a line or polygon from or to longitude/latitude, the distance between the nodes should not be very large. You only used the four points on the image above. These are the extremes in the Lambert Conformal Conical projection, and these are not sufficient to capture the same extent in longitude/latitude coordinates. Hence the use of densify.

A quick and dirty workaround for this case could be

g <- crds(p)
i <- g[,1] > 50    
g[i, 1] <- g[i, 1] - 360
v <- vect(g, "polygons", crs=WGS84)

enter image description here

Now longitude goes beyond -180

v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 1, 0  (geometries, attributes)
# extent      : -184.2079, -0.416633, 48.35129, 88.09466  (xmin, xmax, ymin, ymax)
# coord. ref. : lon/lat WGS 84

And you could correct for that with normalize.longitude:

n <- normalize.longitude(v)
ext(n)
#SpatExtent : -180, 180, 48.3512922956366, 88.0946646700729 (xmin, xmax, ymin, ymax)
plot(n, col="blue", border="blue")

enter image description here

In terra 1.6-21 (currently the development version) you can do the above with rotate

d <- densify(poly_LCC, 100000)
p <- terra::project(d, WGS84)

x <- rotate(p, 50)
# or
y <- rotate(p, 50, normalize=TRUE)

That makes it easier to correct for this problem in some cases. But perhaps there could be a general solution in project to avoid this problem in the first place.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • Prior to this demonstration I had the impression that `normalize.longitude` would operate like `center_around(bisect(v, -90.0)` (an imagined combination of unrealized functions). – Chris Sep 20 '22 at 21:52
  • I have added a method `rotate` that (perhaps) is more like that. – Robert Hijmans Sep 20 '22 at 23:05