rspatial / dismo

R dismo package
GNU General Public License v3.0
25 stars 11 forks source link

circles() fails again (still?) at -180 lon #31

Open plantarum opened 2 years ago

plantarum commented 2 years ago

I think this is the same issue as #4. The example from that issue works with the current code, but the following doesn't:

library(dismo)

tmpPts <- data.frame(decimalLongitude = -179.916667,
                     decimalLatitude = -29.25)

circles(tmpPts, 10000, lonlat = TRUE)

Error: TopologyException: side location conflict at -179.9201118234478 -29.160219827926703. This can occur if the input geometry is invalid.

plantarum commented 1 year ago

I took another look, and I the issue seems to be with .generateCircles, particularly this clause in the 'unwrapping' code:

            } else {
                r1[c(1,nrow(r1)),1] <- -180
                r2[c(1,nrow(r2)),1] <- 180
                pols <- c(pols, Polygons(list(Polygon(r1), Polygon(r2)), i))
            }

As I understand it, you break the polygon up into the negative and positive halves, and in this section, make sure the first and last lines of each sub-polygon are 180 or -180. This works if the sub-polygons are contiguous points, but fails when there are two discontinuous sections in the positive or negative portion. In those cases, you are inserting a 180/-180 coordinate at a random location along the circumference.

I made a correction that addresses this:

        if (diff(range(res[,1])) > 350) {  # basic test for wrapping
            res <- cbind(res, 1:nrow(res)) ## add an index column
            r1 <- res[res[,1] < 0, ,drop=FALSE]
            r2 <- res[res[,1] >= 0, ,drop=FALSE]
            if (nrow(r1) < 3) {
                res[res[,1] < 0, 1] <- res[res[,1] < 0, 1] + 360
                pols <- c(pols, Polygons(list(Polygon(
                                                  res[, 1:2] )), i))
            } else if (nrow(r2) < 3) {
                res[res[,1] >= 0, 1] <- res[res[,1] >= 0, 1] - 360
                pols <- c(pols, Polygons(list(Polygon(
                                                  res[, 1:2] )), i))
            } else {
                                if(sum(diff(r1[, 3]) != 1) == 0){
                                  ## continous arc, use the first and last:
                                  r1[c(1,nrow(r1)),1] <- -180
                                } else {
                                  ## discontinuous arc, find the endpoint:
                                  endpoint <- which(diff(r1[, 3]) != 1)
                                  r1[c(endpoint, endpoint + 1), 1] <- -180
                                }
                                if(sum(diff(r2[, 3]) != 1) == 0){
                                  r2[c(1,nrow(r2)),1] <- 180
                                } else {
                                  endpoint <- which(diff(r2[, 3]) != 1)
                                  r2[c(endpoint, endpoint + 1), 1] <- 180
                                }
                pols <- c(pols,
                                          Polygons(list(Polygon(r1[, 1:2]),
                                                        Polygon(r2[, 1:2])), i))
            }
        } else {
            pols <- c(pols, Polygons(list(Polygon( res )), i))
        }
    }

Not the most elegant, but what I've done is add an index column to res, and use that to determine whether each sub-polygon is a continuous series of points or not. If it isn't, then I used the index to locate the actual breakpoint where the 180/-180 should be inserted. The index column is excluded when the result is passed to Polygon.

If this looks ok to you I can make a pull request.

Best,

Tyler