[R-sig-Geo] Error in MtrxSet ... polygons not (all) closed

Wed Aug 21 19:57:04 CEST 2024

Thanks for the update.
Indeed, I have noticed that this polygon is on the edge of Sentinel-2 tiles. As a result, the AOI gets covered with two different tiles on every overpass, but not completely (bad luck for a 4.5 ha field!). Even the images that display 100% area coverage actually miss a tiny, tiny fraction of it (approximately 4e-10).
catalog <- SearchCatalog(aoi = aoi, from = "2023-06-01", to = "2023-08-31", collection = "sentinel-2-l2a”, with_geometry = TRUE, client = OAuthClient)
100 - subset(catalog, areaCoverage > 99)$areaCoverage
##  [1] 4.085194e-10 4.085194e-10 4.085194e-10 4.085194e-10 4.085194e-10 ...
I'm not sure if this could be the source of the problem, but it's good that you found a workaround.
Best regards,
Zivan

On 21 Aug 2024, at 18:59, Micha Silver <tsvibar using gmail.com> wrote:

Follow up:
After some gnawing of my teeth, I found the problem: Nothing to do with R, sf, nor the CDSE package. My date range included one date where the Sentinel-2 tile from Copernicus data space was damaged/clipped/mangled. And my study area was "hanging over the edge".  I found the problem image using the Copernicus Dataspace browser, and going thru date by date...
I'm not sure why exactly that caused the error, but once I changed the date range to exclude that problem image, the procedure worked as expected.

Thanks for your help.

On 19/08/2024 21:29, Roger Bivand wrote:
> Adding the output of sessionInfo() and sf_extSoftVersion() might also help. The difficulty Zivan points out is that the structure as read is valid, but before conversion to text by dput() may have actually had non-identical first and last coordinates, because of the number of digits used. Also try putting the file created by saveRDS(aoi, ...) on a website with a link here. It would be useful to see the code in the workflow at the point of failure (you point to Error in MtrxSet, maybe also run traceback() after the failure to give more precision in where this is occurring, maybe in stars?
> 
> Hop this helps,
> 
> Roger
> 
> --
> Roger Bivand
> Emeritus Professor
> Norwegian School of Economics
> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway
> Roger.Bivand using nhh.no <mailto:Roger.Bivand using nhh.no>
> 
> ________________________________________
> From: R-sig-Geo <r-sig-geo-bounces using r-project.org> <mailto:r-sig-geo-bounces using r-project.org> on behalf of Zivan Karaman <zivan.karaman using gmail.com> <mailto:zivan.karaman using gmail.com>
> Sent: 19 August 2024 18:40
> To: Micha Silver
> Cc: R-sig-Geo using r-project.org <mailto:R-sig-Geo using r-project.org>
> Subject: Re: [R-sig-Geo] Error in MtrxSet ... polygons not (all) closed
> 
> Could you please provide an example of code that actually raises this error?
> Best regards,
> Zivan
> 
> 
> On 19 Aug 2024, at 16:23, Micha Silver <tsvibar using gmail.com> <mailto:tsvibar using gmail.com> wrote:
> 
> Searching for the error in the subject line returns some discussions from a few years ago. But I'm not able to overcome this error.
> 
> Here's the polygon in question:
> 
>> dput(aoi)
> structure(list(x =
> 
> structure(list(structure(list(structure(c(11.0127043, 11.0127061, 11.0134526, 11.0134554, 11.0134626, 11.0134692, 11.0134776, 11.0134813, 11.0157343, 11.0157373, 11.0157437, 11.0157495, 11.0157546, 11.0157585, 11.015761, 11.0157617, 11.0157616, 11.0157613, 11.0153992, 11.0153983, 11.0153956, 11.0153922, 11.0153869, 11.0153804, 11.0153729, 11.0153694, 11.0127478, 11.0127447, 11.0127379, 11.012732, 11.0127269, 11.0127232, 11.0127212, 11.01272, 11.0127197, 11.0126956, 11.0126958, 11.0126965, 11.0126978, 11.0127004, 11.0127043, 46.8484069, 46.8484048, 46.8476584, 46.847656, 46.8476511, 46.8476481, 46.8476461, 46.8476456, 46.8475213, 46.8475214, 46.847522, 46.8475235, 46.8475266, 46.8475311, 46.8475365, 46.8475425, 46.8475489, 46.8475519, 46.8495623, 46.8495656, 46.8495729, 46.8495788, 46.8495833, 46.8495856, 46.849587, 46.8495874, 46.8496081, 46.8496079, 46.8496069, 46.849605, 46.8496014, 46.8495965, 46.8495906, 46.8495839, 46.8495808, 46.8484309, 46.8484282, 46.8484218, 46.8484166, 46.848412, 46.8484069), dim = c(41L, 2L), dimnames = list( NULL, c("X", "Y")))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = 11.0126956, ymin = 46.8475213, xmax = 11.0157617, ymax = 46.8496081), class = "bbox"), crs = structure(list( input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = 1L, class = c("sf", "data.frame"), sf_column = "x", agr = structure(integer(0), class = "factor", levels = c("constant", "aggregate", "identity"), names = character(0)))
> 
> This polygon is used as area of interest to crop Copernicus imagery that I am downloading using the {CDSE} package. (I have run this workflow successfully many times with several other aoi polygons)
> 
> I have tried:
> 
> 1- sf::st_make_valid()
> 
> 2- transforming to a UTM CRS
> 
> 3- extracting coordinates and recreating the polygon (i.e.:
> aoi_p <- st_polygon(list(st_coordinates(aoi)[,1:2]))
>     aoi3 <- aoi_p |>
>       st_sfc() |>
>       st_as_sf()
>     st_crs(aoi3) <- "EPSG:4326"
> 
> 4- buffering by a small amount
> 
> 
> The above error recurs in all cases (only with this problem polygon). Any suggestions?
> 
> 
> Thanks
> 
> 
> Micha Silver
> Ben Gurion Univ.
> Sde Boker, Remote Sensing Lab
> cell: +972-523-665918
> https://orcid.org/0000-0002-1128-1325
> 
> --
> Micha Silver
> Ben Gurion Univ.
> Sde Boker, Remote Sensing Lab
> cell: +972-523-665918
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org <mailto:R-sig-Geo using r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org <mailto:R-sig-Geo using r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918

	[[alternative HTML version deleted]]