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

Loop 0 is not valid: Non-empty, non-full loops must have at least 3 vertices #2388

Open
TripleEmma opened this issue May 6, 2024 · 7 comments
Labels
reprex needs a minimal reproducible example

Comments

@TripleEmma
Copy link

I have a sf object created like the following (I need the polygons to be equal size), and in order to accommodate with other downloaded data, I convert these polygons back to 'EPSG:4326'. Both objects have the same number of polygons, but different coordinates. Specifically, the latter sf has some polygons with missing vertices. I noticed this when applying st_join() on the latter polygons (grid4326). There is error: Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, Loop 0 is not valid: Non-empty, non-full loops must have at least 3 vertices'.

worldSingleCell <- st_polygon(list(
    rbind(c(-180, -90),
          c(180, -90),
          c(180, 90),
          c(-180, 90),
          c(-180, -90)
    )
))
# create sfc object with CRS as "EPSG:4326"
worldSingleCell_sfc <- st_sfc(worldSingleCell, crs = "EPSG:4326")

# Function to create grids
gridsList <- function(edgeLen, worldSingleCell_sfc){
    edge <- units::as_units(edgeLen, "km")
    grid <- worldSingleCell_sfc %>%  
        st_transform("ESRI:54017") %>% 
        st_make_grid(cellsize = c(edge, edge)) %>% 
        st_sf() %>% 
        st_cast()
    grid$gridNo <- paste0('grid', seq(nrow(grid)))
    return(grid)
}
grid <- gridsList(50, worldSingleCell_sfc)

# covert the grid back to "EPSG:4326" to be consistent with other downloaded data
grid4326 <- st_transform(grid, 'EPSG:4326')

print(grid)
print(grid4326)

# however, the coordinates is different
grid_coordinate <- grid %>%
    st_coordinates() %>%
    as_tibble()

grid4326_coordinate <- grid4326 %>%
    st_coordinates() %>%
    as_tibble()

nrow(grid_coordinate)
nrow(grid4326_coordinate)

>print(grid)
Simple feature collection with 204330 features and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -17367530 ymin: -7342230 xmax: 17382470 ymax: 7357770
Projected CRS: World_Behrmann
First 10 features:
geometry gridNo
1 POLYGON ((-17367530 -734223... grid1
2 POLYGON ((-17317530 -734223... grid2
3 POLYGON ((-17267530 -734223... grid3
4 POLYGON ((-17217530 -734223... grid4
5 POLYGON ((-17167530 -734223... grid5
6 POLYGON ((-17117530 -734223... grid6
7 POLYGON ((-17067530 -734223... grid7
8 POLYGON ((-17017530 -734223... grid8
9 POLYGON ((-16967530 -734223... grid9
10 POLYGON ((-16917530 -734223... grid10

>print(grid4326)
Simple feature collection with 204330 features and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -90 xmax: 179.6366 ymax: 84.47133
Geodetic CRS: WGS 84
First 10 features:
geometry gridNo
1 POLYGON ((-180 -90, -179.48... grid1
2 POLYGON ((-179.4818 -90, -1... grid2
3 POLYGON ((-178.9636 -90, -1... grid3
4 POLYGON ((-178.4454 -90, -1... grid4
5 POLYGON ((-177.9272 -90, -1... grid5
6 POLYGON ((-177.409 -90, -17... grid6
7 POLYGON ((-176.8907 -90, -1... grid7
8 POLYGON ((-176.3725 -90, -1... grid8
9 POLYGON ((-175.8543 -90, -1... grid9
10 POLYGON ((-175.3361 -90, -1... grid10

# however, the coordinates is different
>nrow(grid_coordinate)
[1] 1021650
>nrow(grid4326_coordinate)
[1] 1020260

@edzer
Copy link
Member

edzer commented May 6, 2024

Yes, but the number of polygons returned is the same. Maybe some grid cells touching the pole become triangles in geographical coordinates, and loose one node, hence the smaller number of coordinates?

@TripleEmma
Copy link
Author

Sorry to ask the very basic question but how to fix this then? Thank you!

@edzer
Copy link
Member

edzer commented May 7, 2024

Fix what?

@TripleEmma
Copy link
Author

TripleEmma commented May 7, 2024

There is error when applying st_join() on grid4326.

geneCoor <- geneMetaUsed %>% dplyr::select(ID, longitude, latitude)
geneMetaUsed.sf <- st_as_sf(geneCoor, coords=c("longitude", "latitude"))
st_crs(geneMetaUsed.sf) <- 4326
gridGenes <- st_join(grid4326, geneMetaUsed.sf)  # grid4326 resulted from grid4326 <- st_transform(grid, 'EPSG:4326')

Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, Loop 0 is not valid: Non-empty, non-full loops must have at least 3 vertices

@edzer
Copy link
Member

edzer commented May 7, 2024

Please modify your reproducible example such that it exposes the error; geneMetaUsed is not known.

@TripleEmma
Copy link
Author

TripleEmma commented May 7, 2024

geneMetaUsed is a data frame with several columns. Here I just selected three of them, including 'longitude' and 'latitude'.

geneCoor <- geneMetaUsed %>% dplyr::select(ID, longitude, latitude)
head(geneCoor)

> head(geneCoor)
# A tibble: 6 × 3
ID longitude latitude

1 MK262583_2309838992 172. -40.4
2 MK262588_2308522646 172. -43.2
3 MK262601_2311115935 172. -40.4
4 GQ481466_2309171845 50.4 58.4
5 GQ481467_2310943848 178. 64.2
6 GQ481468_2308837004 34.9 69.0

I once did it as following, and it works fine, no error.

geneCoor <- geneMetaUsed %>% dplyr::select(ID, longitude, latitude)
geneMetaUsed.sf <- st_as_sf(geneCoor, coords=c("longitude", "latitude"))
st_crs(geneMetaUsed.sf) <- 4326
geneMetaUsed.sf <- st_transform(geneMetaUsed.sf, "ESRI:54017")
gridGenes <- st_join(grid, geneMetaUsed.sf)  # CRS of grid is "ESRI:54017"

@edzer edzer added the reprex needs a minimal reproducible example label May 7, 2024
@TripleEmma TripleEmma reopened this May 7, 2024
@TripleEmma
Copy link
Author

@edzer I am sorry if I misunderstood your idea. I had a mistake in the original code, the final geneMetaUsed.sf should be derived from geneCoor (fixed in previous posts). geneCoo has three columns, including "longitude" and "latitude".

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
reprex needs a minimal reproducible example
Projects
None yet
Development

No branches or pull requests

2 participants