Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
plantarum opened this issue Nov 7, 2022 · 1 comment
Open

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

plantarum opened this issue Nov 7, 2022 · 1 comment

Comments

@plantarum
Copy link

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
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant