Why and how to fix Mapview not show all points and st_buffer merges some areas R?
I'm building a map with mapview package. However, when I use st_sample and generate points, some don't get mapped. In addition, st_buffer merges some areas (where they should be separate).
library(sf)
library(mapview)
library(magrittr)
set.seed(145)
chemins.32188 = st_zm(chemins2) %>% st_transform(crs = 32188)
rd.pts.sentiers <- st_sample(x = chemins.32188, size = 10)
buff.pt = st_buffer(x = rd.pts.sentiers,
dist = units::set_units(x = 100,
value = "m"))
mapviewOptions(fgb = FALSE)
mapview(chemins.32188) +
mapview(rd.pts.sentiers) +
mapview(buff.pt)
When I use plot, it works... but the buffer is still merged.
plot(chemins.32188$geom)
plot(buff.pt, add= TRUE, col = scales::alpha("blue", .5))
plot(rd.pts.sentiers, add= TRUE, pch = 19, col = "red")
This is the data
chemins2 = structure(list(Name = "Path Measure", description = NA_character_,
timestamp = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), begin = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), end = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), altitudeMode = NA_character_, tessellate = -1L, extrude = 0L,
visibility = -1L, drawOrder = NA_integer_, icon = NA_character_,
sentier = NA_character_, layer = "path Parc — Path Measure",
path = "~/Desktop/path Parc.kml|layername=Path Measure|geometrytype=LineString25D|uniqueGeometryType=yes",
geom = structure(list(structure(list(structure(c(299434.130314458,
300531.43655987, 299434.130314458, 5046816.92873926, 5047257.53966348,
5046816.92873926), dim = 3:2), structure(c(299434.130314458,
300401.157828111, 299434.130314458, 5046816.92873926, 5046249.12770288,
5046816.92873926), dim = 3:2)), class = c("XY", "MULTILINESTRING",
"sfg"))), n_empty = 0L, crs = structure(list(input = "EPSG:32188",
wkt = "PROJCRS[\"NAD83 / MTM zone 8\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"MTM zone 8\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-73.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9999,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",304800,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting (E(X))\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (N(Y))\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Canada - Quebec between 75°W and 72°W.; Canada - Ontario - east of 75°W.\"],\n BBOX[44.98,-75,62.53,-72]],\n ID[\"EPSG\",32188]]"), class = "crs"), class = c("sfc_MULTILINESTRING",
"sfc"), precision = 0, bbox = structure(c(xmin = 299434.130314458,
ymin = 5046249.12770288, xmax = 300531.43655987, ymax = 5047257.53966348
), class = "bbox"))), row.names = 1L, sf_column = "geom", agr = structure(c(Name = NA_integer_,
description = NA_integer_, timestamp = NA_integer_, begin = NA_integer_,
end = NA_integer_, altitudeMode = NA_integer_, tessellate = NA_integer_,
extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_,
icon = NA_integer_, sentier = NA_integer_, layer = NA_integer_,
path = NA_integer_), class = "factor", levels = c("constant",
"aggregate", "identity")), class = c("sf", "data.frame"))
Solution
In your reprex you have MULTILINESTRING of 2 LINESTRINGs, st_sample() returns 2 MULTIPOINTs and with that particular seed, the first one happens to be empty (i.e. all points are from the same sub-LINESTRING):
set.seed(145)
chemins.32188 = st_zm(chemins2) %>% st_transform(crs = 32188)
st_as_text(chemins.32188$geom)
#> [1] "MULTILINESTRING ((299434.1 5046817, 300531.4 5047258, 299434.1 5046817), (299434.1 5046817, 300401.2 5046249, 299434.1 5046817))"
rd.pts.sentiers <- st_sample(x = chemins.32188, size = 10)
rd.pts.sentiers
#> Geometry set for 2 features (with 1 geometry empty)
#> Geometry type: MULTIPOINT
#> Dimension: XY
#> Bounding box: xmin: 299516.7 ymin: 5046255 xmax: 300390.5 ymax: 5046768
#> Projected CRS: NAD83 / MTM zone 8
#> MULTIPOINT EMPTY
#> MULTIPOINT ((300390.5 5046255), (299516.7 50467...
Apparently mapview (or leaflet) is not too happy about empty (point) geometries and you should also receive the following warning from mapview():
Warning message:
In validateCoords(lng, lat, funcName) :
Data contains 1 rows with either missing or invalid lat/lon values and will be ignored
And when you create a buffer from MULTIPOINT rd.pts.sentiers, each feature will be a MULTIPOLYGON, so overlaps within a single feature will cause a polygon merge.
To fix, you might want to check samples for empty geometries and cast samples to POINTs before creating buffers:
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(mapview)
library(magrittr)
set.seed(145)
chemins.32188 <- st_zm(chemins2) %>% st_transform(crs = 32188)
rd.pts.sentiers <- st_sample(x = chemins.32188, size = 10)
rd.pts.sentiers
#> Geometry set for 2 features (with 1 geometry empty)
#> Geometry type: MULTIPOINT
#> Dimension: XY
#> Bounding box: xmin: 299516.7 ymin: 5046255 xmax: 300390.5 ymax: 5046768
#> Projected CRS: NAD83 / MTM zone 8
#> MULTIPOINT EMPTY
#> MULTIPOINT ((300390.5 5046255), (299516.7 50467...
# keep only non-empty samples, cast from MULTIPOINT to POINT
rd.pts.sentiers <-
rd.pts.sentiers[!st_is_empty(rd.pts.sentiers)] %>%
st_cast("POINT")
rd.pts.sentiers
#> Geometry set for 10 features
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 299516.7 ymin: 5046255 xmax: 300390.5 ymax: 5046768
#> Projected CRS: NAD83 / MTM zone 8
#> First 5 geometries:
#> POINT (300390.5 5046255)
#> POINT (299516.7 5046768)
#> POINT (300002.8 5046483)
#> POINT (299737 5046639)
#> POINT (300193.1 5046371)
# get 10 POLYGON features instead of 2(1) MULTIPOLYGONs
buff.pt = st_buffer(x = rd.pts.sentiers,
dist = units::set_units(x = 100,
value = "m"))
buff.pt
#> Geometry set for 10 features
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 299416.7 ymin: 5046155 xmax: 300490.5 ymax: 5046868
#> Projected CRS: NAD83 / MTM zone 8
#> First 5 geometries:
#> POLYGON ((300490.5 5046255, 300490.3 5046250, 3...
#> POLYGON ((299616.7 5046768, 299616.5 5046763, 2...
#> POLYGON ((300102.8 5046483, 300102.7 5046478, 3...
#> POLYGON ((299837 5046639, 299836.9 5046634, 299...
#> POLYGON ((300293.1 5046371, 300293 5046366, 300...
mapviewOptions(fgb = FALSE)
mapview(chemins.32188) +
mapview(rd.pts.sentiers) +
mapview(buff.pt)
Created on 2024-04-23 with reprex v2.1.0
- How to get the selected tab's name of a Quarto Tabset Panel?
- Separting alphanumeric string using tidyr separate wider regex
- Convert JSON type scraped string to dataframe in R
- Testing private methods in R6 classes in R
- How can I align the tag with the displayed text?
- Date labels in the x axis are misaligned by one month
- How can I force a shiny app to switch colors after a specific time?
- pretty_breaks not working with geom_spatraster
- pretty_breaks does not working with geom_sf
- overwrite a facet panel with custom plot
- How to draw a shapefile with ggpolar
- Exported file from R using `haven` cannot be opened by SAS
- How to install libpng-dev on Mac
- Ignoring NAs while aggregating character vector
- Shinyselected row in a striped DT::datatable() disappears when style = "bootstrap5" and RowGroup extension is used
- ggplot2 secondary axis tidyeval
- parse dynamic case_when based on thresholds values
- dbWriteTable or dbAppendTable fail for date column when appending to DB2-table
- How can I use a loop in RMarkdown to programmatically plot figures and text on individual pages in a document?
- Conditional recursive formula in dplyr
- Current development in how to adjust individual x axis title using ggplot facet_wrap
- How to make a probability figure using a logistic regression?
- dplyr - Need to add the CSV filename after importing all files from a directory
- How to reorder a stacked bar chart in ggplot2?
- how can I render a quarto website to a single pdf
- Assign specific values to unique values in all dataframe columns in R
- Unable to install tidyverse package for r kernel in jupyter
- How to replace values in a matrix based on a conditional variable
- How can I make my enormous dataframe from wide to long in R?
- Group variable stacking when using geom_ridgeline_gradient()