r-spatial / sf

Simple Features for R
https://r-spatial.github.io/sf/
Other
1.34k stars 298 forks source link

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

Open TripleEmma opened 6 months ago

TripleEmma commented 6 months ago

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 commented 6 months ago

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 commented 6 months ago

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

edzer commented 6 months ago

Fix what?

TripleEmma commented 6 months ago

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 commented 6 months ago

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

TripleEmma commented 6 months ago

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"
TripleEmma commented 6 months ago

@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".