Querying elements from osm and doing practical things with them is really a useful and fun way of code “hacking” for me. In this case, I wanted to show farmers markets with a reactive map like leaflet like I’ve done before in my journal entries. In this case, I wanted to make it easy to browse the times, seasons, and contact info for farmers markets in the Anchorage, Alaska area. I think with some refinement this map could be useful to sellers and buyers to determine which farmers markets they can fit into their schedules.

library(tidyverse)
library(osmdata)
library(leaflet)
library(htmlwidgets)


marketplace <- opq(c(xmin = -150.494625, ymin = 60.771028, xmax = -148.866047, ymax = 61.550720)) %>%
  add_osm_feature("amenity", "marketplace") %>%
  osmdata_sf()

marketplace_points <- marketplace$osm_points

greenLeafIcon <- makeIcon(
  iconUrl = "https://leafletjs.com/examples/custom-icons/leaf-green.png",
  iconWidth = 38, iconHeight = 95,
  iconAnchorX = 22, iconAnchorY = 94,
  shadowUrl = "https://leafletjs.com/examples/custom-icons/leaf-shadow.png",
  shadowWidth = 50, shadowHeight = 64,
  shadowAnchorX = 4, shadowAnchorY = 62
)

content <- paste(sep = "<br/>",
                 paste0("<b><a href='", marketplace_points$website,"'>", marketplace_points$name, "</a></b>"),
                 paste("<i>", marketplace_points$description, "</i>"),
                 #paste(marketplace_points$addr.housenumber, marketplace_points$addr.street),
                 #paste(marketplace_points$addr.city, marketplace_points$addr.state, marketplace_points$addr.postcode),
                 paste("Phone:", marketplace_points$phone),
                 paste("Email:", marketplace_points$email),
                 #paste("Hours:", marketplace_points$opening_hours),
                 paste( paste0("<b><a href='", marketplace_points$website,"'>", "Website", "</a></b>"))
)


(m <- leaflet(marketplace$osm_points) %>%
  addTiles() %>%
  addMarkers(icon = greenLeafIcon,
             popup = lapply(content, htmltools::HTML)))

saveWidget(m, "farmers_markets.html")
#write.csv(marketplace_points, file = "marketplace_points.csv")

Mountain names. A new front in the never-ending cacophony of US political discourse.

Two recent mountain peaks had their names changed by the US Board on Geographical names and there is disagreement in the change sets on whether the name given should be the new indigenous name given by the US Board or the English translation of the US Board. To be clear though, I have not seen the official name from a primary source.

The two peaks were promptly edited to their new names in the local Dena’ina (North and South Yuyanq’ Ch’ex) but within 24 hours changed to what I assume is the English translation of the new Dena’ina name (North and South Heaven’s Breath).

To me it seems like a cut and dry answer based on osm standards what should fall in the “name=” attribute. The editor that changed it to the English translation may feel otherwise though so and to avoid a change set dispute I emailed the Data Working Group to get their decision. I haven’t heard from them in over 10 days however and I haven’t received a message to my follow up email on whether my email got to them. If anyone is on the DWG is interested in bringing me on to handle these issues, I would volunteer. For now though, I’ll hold off on changing the name back and hold faith in the process.

news article:

https://www.adn.com/alaska-news/anchorage/2022/10/14/anchorages-suicide-peaks-receive-new-denaina-name-of-yuyanq-chex/

relevant changesets:

https://www.openstreetmap.org/changeset/127584704

https://www.openstreetmap.org/changeset/127643552

https://www.openstreetmap.org/changeset/127682614

While I never expected I would put so much work in to the Chugach mountains, it does look very nice I believe. But I have to think about where I want to end it or what’s a good place to take a break.

Capturing the biome-like areas from woods to shrubs to heath to scree to bare rock, it’s a bit of a approximation when there are hybrid areas. I’ve browsed different mountain areas around the globe but it’s not clear who else is thinking of capturing all these together.

My first priority is to finish the Eagle River valley and valley in north west Girdwood valley where they follow the Crow Pass trail.

My second priority is to finish up the area around penguin creek and then the missing middle section.

Then, I think I should map the Baldy area in Eagle River from the top to heath and shrubs.

And finally, when I can finish the Peters Creek and Eklutna valleys, I’ll think I’ll feel like I can put a wrap on it.

Good tidings!

My effort to add the natural features to Nancy Lake is finished. https://www.openstreetmap.org/#map=12/61.6291/-150.0966

What was the point of adding these natural features to a wilderness?

Well, besides being one of my favorite spots in Alaska, mostly I’d say it was just another project like mapping Anchorage where I could incrementally improve and feel some accomplishment at the end. At the end of the day, I hope these edits provided value but it’s my goal to be driven by the small neurochemical treat-reward cycle of editing OSM instead of doing nothing which would be less valuable than something.

My humble pandemic project has been to improve Anchorage, Alaska on osm. In the spirit of Isaac Newton in the time of the plague in London I would try and do a bigger push on a topic and while I didn’t invent the next calculus, I did improve Anchorage quite a bit. As we gently start coming out of the pandemic I figured it’s time to at least put a bow on this as a project and list improvements I’ve made that I’m proud of.

1. Adding buildings using the mapwith.ai
2. adding homeowner’s associations and neighborhoods
3. Adding natural features to the southern front range of Chugach State Park (woods, scrub, heath, scree, bare_rock).
4. Adding all sorts of trails (snowmachine trails north of big lake, hiking trails in chugach state park, connectors that were missed in kinkaid park, trails for smaller parks and HOAs in Anchorage, and little trail easements for elementary school access).
5. Adding retail, commercial, and industrial zoning.
6. Fixing TIGER import oopsies that have been around for a long time. Mostly streets that extend into woods that don’t really exist.
7. school details (playgrounds, ice rinks, sports fields)
8. lots and lots of parking
9. pockets of woods, wetland, and scrub in undeveloped large lots.
10. a pretty comprehensive meta tagging of totem poles with images.
11. natural features near airports
12. religious areas.

Some challenges remain:

1. golf and disc golf courses
2. the Alaska Zoo
3. quasi-neighborhoods that have no name
4. sidewalks, crossings, and footpaths in commerical areas.
5. filling out the zones in kinkaid that are somewhere between woods, scrub, and grass.
6. the woods that innervate neighborhoods.

In time, I feel society will view my additions on an equal merit with Newton’s quarantine achievements. /s

I’ve wanted to merge two ideas for a while. I did so a month or so back and wanted to share via a diary page (which is my sole programming blog outlet at this time).

I want to show a webmap that has overlapping isochrones where it shows the shortest distance between multiple zones. Now, this was sort of possible for me before with vectorized rings but I want it as a raster so I can show the spectrum of time it would take. It makes it more detailed and I think more pretty. Below is an example using the mapbox api for doing so. I assume like most web services it will be deprecated in the future but who knows.

library(mapboxapi)
library(tidyverse)
library(mapdeck)
library(osmdata)
library(leaflet)
library(sf)
library(raster)
library(rgeos)
library(fasterize)
library(rgdal)
library(htmlwidgets)

#mapboxapi::mb_access_token("put_your_own_key_here",
#                           install = TRUE, overwrite = TRUE)


x <- getbb("Anchorage, AK") %>% opq() %>%
  add_osm_feature("shop", "supermarket") %>%
  osmdata_sf()
#y <- getbb("Anchorage, AK") %>% opq() %>%
#  add_osm_feature("shop", "wholesale") %>%
#  osmdata_sf()
#x <- bind_rows(x$osm_polygons, y$osm_polygons)

x <- x$osm_polygons

all_latlon <- c()
for (i in x$osm_id) {
  
  latlon <-
    paste0(x$geometry[[i]][1][[1]][1,1], ", ",
           x$geometry[[i]][1][[1]][1,2])
  
  all_latlon <- c(latlon, all_latlon)
}

data.frame(name = x$name, latlon = all_latlon)



all_isos <- data.frame()

for(j in all_latlon) {
  isos <- mb_isochrone(
    location = j,
    profile = "walking",
    time = 1:20,
  )  
  print(j)
  all_isos <- rbind(all_isos, isos)
  Sys.sleep(1)
}


library(sf)
library(dplyr)

iso_sf <- st_as_sf(all_isos)

iso_union <- iso_sf %>%
  group_by(time) %>%
  summarise()

isos_proj <- st_transform(iso_sf, 32615)

template <- raster(isos_proj, resolution = 100)

iso_surface <- fasterize(isos_proj, template, field = "time", fun = "min")

pal <- colorNumeric("viridis", isos_proj$time, na.color = "transparent")

(map <- leaflet() %>%
  addTiles() %>%
  addRasterImage(iso_surface, colors = pal, opacity = 0.5) %>%
  addLegend(values = isos_proj$time, pal = pal,
            title = "Minutes of Travel") %>%
  addPolygons(data = x, popup = ~name,
              opacity = 1.0, fillOpacity = 0.5, 
              color = "orange"))

writeRaster(iso_surface, filename = "supermarkets.tif")


saveWidget(map, file="grocery_stores.html")

There is an excellent blog post from a very curious and original person that made geographic vonoroni maps in R that I’ve adapted to use with OSM attributes and geometery.

Blog post: https://rstudio-pubs-static.s3.amazonaws.com/97256_125b33cf42414eda8797f53ac632628b.html

Code for doing this in OSM with a boundary shapefile.

library(sp)
library(rgdal)
library(deldir)
library(ggplot2)
library(ggthemes)
library(osmdata)
library(leaflet)
library(dplyr)
library(rgeos)
library(sf)
library(raster)

SPointsDF_to_voronoi_SPolysDF <- function(sp) {

  # tile.list extracts the polygon data from the deldir computation
  vor_desc <- tile.list(deldir(sp@coords[,1], sp@coords[,2]))

  lapply(1:(length(vor_desc)), function(i) {

    # tile.list gets us the points for the polygons but we
    # still have to close them, hence the need for the rbind
    tmp <- cbind(vor_desc[[i]]$x, vor_desc[[i]]$y)
    tmp <- rbind(tmp, tmp[1,])

    # now we can make the Polygon(s)
    Polygons(list(Polygon(tmp)), ID=i)

  }) -> vor_polygons

  # hopefully the caller passed in good metadata!
  sp_dat <- sp@data

  # this way the IDs _should_ match up w/the data & voronoi polys
  rownames(sp_dat) <- sapply(slot(SpatialPolygons(vor_polygons),
                                  'polygons'),
                             slot, 'ID')

  SpatialPolygonsDataFrame(SpatialPolygons(vor_polygons),
                           data=sp_dat)

}

### load Anchorage bourdary file
anc_boundary <- st_read(
  "Municipality_of_Anchorage_Boundary.shp")
anc_boundary <- as(anc_boundary, Class = "Spatial") # conver sf to sp

# query overpass API
x <- opq(bbox = c(-150.090972, 61.470444, -149.321814, 60.034781)) %>% 
  add_osm_feature(key = 'amenity', value = "place_of_worship", value_exact = FALSE) %>%
  osmdata_sf()

### Make vonroni polygons
pow_points <- st_centroid(x$osm_polygons)
pow_points <- as(pow_points, Class = "Spatial")
vor <- SPointsDF_to_voronoi_SPolysDF(pow_points)
vor <- intersect(anc_boundary, vor)
vor_df <- fortify(vor)

### Make Popup
y <- x$osm_polygons %>% 
  mutate(name = ifelse(!is.na(website), paste0('<a href="', website, '">', name, '</a>'), name),
         religion = paste0("religion: ", religion),
         denomination = ifelse(denomination != "none", paste0("denomination: ", denomination), ""),
         phone = ifelse(!is.na(website), paste0("phone #: ", phone), ""),
         email = ifelse(!is.na(website), paste0("email: ", email), "")) %>%
  mutate(for_popup = paste(name, phone, email, sep = "<br>"))

plot(intersect(anc_boundary, vor))


(map <- leaflet() %>% addTiles() %>% 
  setView(lng = -149.845, lat = 61.15, zoom = 12) %>%
  addPolygons(data = vor, fill = TRUE, fillOpacity = 0.1, weight = 1) %>% 
  addPolygons(data = y, color = "black", popup = ~for_popup, stroke = TRUE, fill = TRUE, weight = 3, 
              fillOpacity = 0.5, fillColor = "white",  options = popupOptions(autoPan = TRUE, closeButton = FALSE)))


library(htmlwidgets)

saveWidget(map, file = "anc_places_of_worship.html")

I’ve done a good amount of work to add basic geometry to Anchorage in the past. Now, I would like to store more data inside the polygons for buildings that have been created.

It is a difficult task to do quality control on the wide number of churches, schools, businesses, etc. in Anchorage because it is just too big. 100’s of buildings in each category where there are also points or multipolgons as well. To accomplish a definitive test of which polygons of a certain type have a data tag, I’ll need to use the overpass API. Fortunately, the work has mostly been done to use it as a tool by relying on R and the library osmdata.

For example, I’ve been able to query all of the amenities with “place_of_worship” to see what tags are available, and map the polygons on top of an osm basemap with the following code.

library(osmdata)
library(leaflet)
library(dplyr)
x <- opq(bbox = c(-150.092146, 61.483194, -148.981596, 60.814537)) %>% 
  add_osm_feature(key = 'amenity', value = "place_of_worship", value_exact = FALSE) %>%
  osmdata_sf()
y <- x$osm_polygons %>% filter(is.na(denomination))
leaflet(y) %>% addTiles() %>% addPolygons(popup = ~y$name)
y %>% as.data.frame() %>% select(name, email, website, phone, denomination) 

…This will show me which places of worship still need to have a denomination added. And I reuse this for different tags I want to add.

I set out to use the rapid.ai interface to quickly bring in building data into openstreetmap and I’ve achieved the goal I sent out for. The greater majority of Anchorage buildings have been added.

Along the way I’ve added trails, schools, HOAs and apartment complexes and more features.

I’ll still add things here and there but I’m going to relax the speed and vigilance I’ve done to add features since I’ve accomplished the goal I set out for.

Add walking paths, enclosures, etc. to Alaska Zoo

• HOA areas.

• South Addition Alley paved/unpaved metadata

small trails
ethnic restaurants
American restaurants
chain restaurants
parking lots
snowmachine trails (nancy lake region)
buildings
state and city offices