How to avoid connecting polygon vertices in the wrong direction when projecting a polygon?

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)

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))

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)

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")

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.