a list of use cases

[tim-salabim]

A simple list of existing R workflows where we think that mapedit could be useful

• HistMapR -> function click_sample
• Atlantis ecosystem models -> creating/editing box geometry existing solution is BGM -> SHP -> QGIS->SHP->BGM
• SGAT -> initial path edits
• SOOS field projects database -> Drawing transects or polygons
• MODIS -> getTile.
• digitizing features

[mdsumner]

• Atlantis ecosystem models -> creating/editing box geometry existing solution is BGM -> SHP -> QGIS->SHP->BGM
• SGAT -> initial path edits
• SOOS field projects database -> Drawing transects or polygons

Defining general transect or polygon extraction for localized studies is a really important gap for us. Nuanced edits with interactive information about coastlines and environmental data is critical.

[tim-salabim]

Do you have any links to existing solutions/approaches for any of these?

[mdsumner]

Yes, sorry was caught out with an unfinished post! There are so many uses for this for us, I should be trying it out more and making examples. It's a bit unclear to me what you need here, lots of examples to motivate or just a few key specific tasks to solve? Creating/editing a track in space-time is something I would use constantly and give to others to use.

[fdetsch]

When downloading higher-level MODIS products using runGdal (or getHdf) from MODIS, users are currently required to

• either specify horizontal and vertical tile indices in accordance with the MODIS Sinusoidal grid (i.e., look those up manually)
• or, alternatively, provide a 'Spatial' or 'Raster' object from which IDs of intersecting tiles are identified in the background via getTile.

I think it would be great if we could offer a more flexible, interactive way of how to select tiles, which especially involves

• selecting single tiles by a) displaying the global MODIS grid on a leaflet map, b) clicking relevant polygons and c) returning sf objects back to R
• as well as allowing the user to draw polygons from which to derive corresponding tile IDs under the hood.

[tim-salabim]

Great idea! Can you provide a sample tile grid data set (preferably as sf)?

[edzer]

Cool idea; I've been thinking about the idea of handling tile collections quite a bit, lately.

A set of tiles for sentinel 2 is loaded by:

tiles = st_read("https://sentinel.esa.int/documents/247904/1955685/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml")
Reading layer `Features' from data source `https://sentinel.esa.int/documents/247904/1955685/S2A_OPER_GIP_TILPAR_MPC__20151209T095117_V20150622T000000_21000101T000000_B00.kml' using driver `KML'
converted into: GEOMETRYCOLLECTION
Simple feature collection with 56686 features and 2 fields
geometry type:  GEOMETRYCOLLECTION
dimension:      XYZ
bbox:           xmin: -180 ymin: -83.83595 xmax: 180 ymax: 84.64428
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
plot(tiles[1], xlim = c(0,10), ylim = c(0,10))

These are UTM tiles, laid out in longlat; it would be great if we could handle (subsets of) these in R. My current idea is to do this the way PostGIS raster does this: in a table, with each record containing, or referring, to a single tile, and its metadata (e.g. outline as an sfg).

[mdsumner]

I was just messing around with the sinu grid, getting their largest tiles as per the diagram: https://gist.github.com/mdsumner/19b36c829ec1cbedabff1bf619b4ec18

It would be interesting to abstract that to zoom as well, but I think easiest would be to store empty raster layers the grid contstructs, and use the cell logic to generate them.

I see the long instructions for manually creating a grid for the poles with the .hdr files too, happy to help iron that out if needed.

[fdetsch]

An EPSG:4326 version of the MODIS tile grid is shipped with the package. It originates from before my time, so I am not quite sure where Matteo took it from. But anyway, it can be loaded into R using

library(sf)
library(mapview)

if (!require(MODIS)) 
  devtools::install_github("MatMatt/MODIS", ref = "develop")

shp <- system.file("external", "modis_latlonWGS84_grid_world.shp", package = "MODIS")
grd <- st_read(shp)

m <- mapview(grd)
m

[mdsumner]

Thanks @fdetsch - that is helpful - you can make more sense of the plot when you do this:

plot(st_transform(grd, st_crs("+proj=sinu"))

and then to get it into mapview with sufficient density of vertices along edges, this works but there's only so far we can push the sense of mapping in Web Mercator here:

stz <- st_segmentize(st_set_crs(grd, NA), dfMaxLength=.05)
mapview(st_set_crs(stz, 4326))

(You can't segmentize with the CRS set because sf will follow great circles, so distance calc now has to occur in native space for non-orthodromes).

Note that this regular tiling in the sinusoidal projection is anathemic to the integrized sinusoidal grid used by "L3bin" products. I don't expect that will come up in this context, but it's a terminology-clash to watch for.

[fdetsch]

Thanks for pointing this out @mdsumner, very useful stuff!

[tim-salabim]

I had a quick play w the MODIS example and this seems to work (Note that we need to set a class character grouping in order to get this working)

grd$cat = as.character(grd$cat)

m = leaflet() %>% 
    addTiles() %>% 
    addPolygons(data = grd, weight = 1, group = ~cat)

selectMap(m,
           styleFalse = list(weight = 1),
           styleTrue = list(weight = 4))

[tim-salabim]

With 4dfd8d94edfe5378e4b53c36003b63c04ee5a439 we can now do

library(mapview)
library(sf)

subs = selectFeatures(breweries, styleTrue = list(fillColor = "red"))
## or
subs = selectFeatures(trails, styleTrue = list(color = "red"))
## or
subs = selectFeatures(franconia, styleTrue = list(fillColor = "red"))
## which is similar to 
subs = selectFeatures(grd, styleTrue = list(fillColor = "red")) # taking grd from the example above

All of these will return a sf object, hence

mapview(subs)

should work for all of the above.

This is, however, just a rough first implementation as there are still many assumptions and it's not at all flexible as of now. Still, this may be one road we could take. Feedback is more than welcome! Especially from @timelyportfolio

[timelyportfolio]

Can't wait to check it out. Thanks everyone!

[fdetsch]

Looks great, thanks everyone. I'll create a new MODIS branch at https://github.com/MatMatt/MODIS and add this example to getTile(). This will definitely find its way into our talk at this year's useR! (:

[timelyportfolio]

@tim-salabim, are you happy with the styleTrue argument? I have not found any better pattern, and I love the flexibility, but this does require some knowledge of CSS.

[tim-salabim]

@timelyportfolio yes, people should be familiar with the CSS parameters from the likes of leaflet::pathOptions and especially leaflet::highlightOptions. List semantic is also fine. We can provide a suitable function to build the list at any stage without problem

[timelyportfolio]

@fdetsch, any feedback on @tim-salabim's very nice first approach? Thanks again for this very helpful use case. As I tried to think of a very simple illustration within the sf ecosystem, this was one simplified example that might demonstrate the functionality.

# example ?st_make_grid
library(sf)
library(mapview)
library(mapedit)

grd <- st_make_grid(what="polygons")
grd_select <- selectFeatures(
  st_sf(grd)
)

plot(grd)
plot(grd_select, col="pink", add=TRUE)

[timelyportfolio]

Here similar to above example with st_voronoi.

# example ?st_voronoi
x = st_multipoint(matrix(runif(30),,2))
box = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
if (sf_extSoftVersion()["GEOS"] >= "3.5.0") {
  v = st_sfc(st_voronoi(x, st_sfc(box)))
  plot(v, col = 0, border = 1, axes = TRUE)
  plot(box, add = TRUE, col = 0, border = 1) # a larger box is returned, as documented
  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
  plot(st_intersection(st_cast(v), box)) # clip to smaller box
  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
}
# select voronoi polygons
voronoi_selected <- selectFeatures(
  st_sf(st_cast(v)),
  styleFalse = list(fillColor="transparent", weight=1),
  styleTrue = list(fillColor="pink", weight=3)
)
plot(voronoi_selected)

[mdsumner]

EDIT: discussion moved out to https://github.com/r-spatial/mapedit/issues/25

[fdetsch]

@timelyportfolio, I've opened a new MODIS branch called mapedit and implemented Tim's approach in getTile(). It works every now and then, but then again, it fails at least as often without obvious reason. But see for yourselves:

[timelyportfolio]

thanks @fdetsch for the screenshot. I will explore and debug. Sorry for the trouble.

[fdetsch]

Sounds good, thank you! In addition, in order to expand interactive tile selection a little bit, it would also be great if we could have editMap() return an object of class 'sf' (eg polygons, lines) based on which we would then be able to identify intersecting MODIS tiles.

[timelyportfolio]

@fdetsch, editMap() returns sf as of https://github.com/r-spatial/mapedit/pull/18. However, I did not realize when I wrote the geojson -> sf conversion that there was already a way in sfr to convert (see issue comment). This might offer a more robust solution, but I have not had the opportunity to test.

[timelyportfolio]

@fdetsch, @tim-salabim promoted selectFeatures to its own issue #24.

[tim-salabim]

@mdsumner I moved your 'fresh' thouhgts to #25

[mdsumner]

With a colleague we just put together a workflow that changed a very manual tiling and image-exploration task into a straightforward set of steps, useable for teaching and real work. This started as a question about how to tile an image with a buffer, and 45 minutes later we had a working prototype ready for some final tweaks. :)

Here is all the code. The list of files is not reproducible, it's just a set of GeoTIFFs on a local drive (in UTM of various zones).

The steps are

1. load an image and produce a tile index of it at given grain (500m)
2. set up a sampled subset of the tiles to be intensively searched
3. set up the mapview object with the right setting (the image scaled appropriately, and the tile region overlad)
4. user interaction, the user interacts with one of the sampled tiles to search the image for features (seals in this case), they can see the tile border as well as a 50m buffer of the larger image outside of that tile, in order to avoid edge-effects and possibly missing counting a feature
5. follow up with the drawn output, plot independently in another context

Thanks! This is so awesome.

library(mapedit)
library(mapview)
library(raster)
library(spex)
img.list <- dir(pattern=".TIF")
p1 <- raster(img.list[1])

## 1. -----------------------------------------------------------------------------
## grid of tiles at 500m grain
grain <- 500
buf <- 50
tiles <- qm_rasterToPolygons_sp(raster(buffer_extent(extent(p1), grain), 
                                       crs = projection(p1), res = grain))

## 2. -----------------------------------------------------------------------------
## we treat this as a sampling exercise, a subset of the tiles
## wll be searched intensively to illustrate issues of study and sampling design
nsamples <- 10
tiles$sample <- FALSE
tiles$sample[sample(nrow(tiles), nsamples)] <- TRUE
#plot(tiles, col = factor(tiles$sample + 1))

# 3. -----------------------------------------------------------------------------
## presentation of the map
ss <- which(tiles$sample)
i <- 1  ## of sum(tiles$sample)
tile_image <-crop(p1, extent(tiles[ss[i], ]) + buf)
col_map <- grey(seq(0, 1, length = 256))
val_map <- seq(cellStats(p1, min), cellStats(p1, max))
map1 <- mapview(tile_image, maxpixels = ncell(tile_image), 
                col.regions = col_map, 
                values = val_map) + mapview(as(tiles[ss[i], ], "SpatialLines"))

# 4. user interaction -----------------------------------------------------------------------------
## the editing session, user sees the "focus tile" as well as the image within
## that tile, but also an extra border of image to width "buf"
## the user's task is to "find seals" and click on where they are
em <- map1 %>% editMap()

# 5. follow up -----------------------------------------------------------------------------
## remember we get the drawn result in 4326 so need to transform
tile_i_points <- sf::st_transform(em$drawn, projection(p1))
plot(tile_image)
plot(tile_i_points, add = TRUE, col = "hotpink")