R/openmap.R

In mapmiscTerra: Utilities for Producing Maps

osmTiles = function(name, xyz, suffix) {
  result = c(
    osm = "http://tile.openstreetmap.org",
    'osm-admin' = 'http://korona.geog.uni-heidelberg.de/tiles/adminb',
    'osm-roads-grey' = 'http://korona.geog.uni-heidelberg.de/tiles/roadsg/',
    'osm-roads' = 'http://korona.geog.uni-heidelberg.de/tiles/roads',
    'osm-semitransparent' = 'http://korona.geog.uni-heidelberg.de/tiles/hybrid/',
    "osm-no-labels"="http://c.tiles.wmflabs.org/osm-no-labels/",
    "osm-de"="http://c.tile.openstreetmap.de/tiles/osmde/",
    "osm-ru" = "http://a.tiles.wmflabs.org/osm-multilingual/ru,_/",
    "osm-transport"="http://tile2.opencyclemap.org/transport/",
    "stamen-toner" = "https://stamen-tiles-d.a.ssl.fastly.net/toner/",
    "stamen-watercolor" = "https://tiles.stadiamaps.com/styles/stamen_watercolor/",
    "bw-mapnik"="http://b.tiles.wmflabs.org/bw-mapnik2/",
#			mapquest="http://otile1.mqcdn.com/tiles/1.0.0/osm/",
#			"mapquest-sat"="http://otile1.mqcdn.com/tiles/1.0.0/sat",
#      "mapquest-labels"='http://otile3.mqcdn.com/tiles/1.0.0/hyb/',
    'osm-cyclemap' = 'http://a.tile.opencyclemap.org/cycle/',
    'osm-seamap' = 'http://tiles.openseamap.org/seamark/',
    'osm-fr' = 'http://a.tile.openstreetmap.fr/osmfr/',
    'landscape'="http://tile.opencyclemap.org/landscape/",
    'hyda' = 'http://c.tile.openstreetmap.se/hydda/full/',
    'hyda-base' = 'http://c.tile.openstreetmap.se/hydda/base/',
    'hyda-roads' = 'http://c.tile.openstreetmap.se/hydda/roads_and_labels/',
    "opentopomap" = "http://opentopomap.org/",
    "maptoolkit"="http://tile2.maptoolkit.net/terrain/",
    waze="https://worldtiles2.waze.com/tiles/",
    'waze-us'='https://livemap-tiles2.waze.com/tiles/',
    humanitarian="http://a.tile.openstreetmap.fr/hot/",
    cartodb='http://c.basemaps.cartocdn.com/light_all/',
    'cartodb-dark'='http://c.basemaps.cartocdn.com/dark_all/',
#  historical='http://www.openhistoricalmap.org/ohm_tiles/',
    nrcan = 
    'http://geoappext.nrcan.gc.ca/arcgis/rest/services/BaseMaps/CBMT_CBCT_GEOM_3857/MapServer/tile/',
    'nrcan-text' = 
    'http://geoappext.nrcan.gc.ca/arcgis/rest/services/BaseMaps/CBMT_TXT_3857/MapServer/tile/',
    'nrcan-text-fr' = 
    'http://geoappext.nrcan.gc.ca/arcgis/rest/services/BaseMaps/CBCT_TXT_3857/MapServer/tile/',
    spinal = 'http://c.tile.thunderforest.com/spinal-map/',
    neighbourhood = 'https://a.tile.thunderforest.com/neighbourhood/',
    pioneer = 'https://b.tile.thunderforest.com/pioneer/',
    'mobile-atlas'='https://b.tile.thunderforest.com/mobile-atlas/',
    wikimedia = 'https://maps.wikimedia.org/osm-intl/',
    'sputnik' = 'http://tiles.maps.sputnik.ru/'
  )

  # toronto
  #https://gis.toronto.ca/arcgis/rest/services/basemap/cot_topo/MapServer/tile/9/186/142
# https://map.toronto.ca/maps/map.jsp?app=TorontoMaps_v2

#	skobbler="http://tiles3.skobbler.net/osm_tiles2/",	
#		"osm2world"="http://tiles.osm2world.org/osm/pngtiles/n/",
#		bvg="http://mobil.bvg.de/tiles/",
#	landshaded="http://tiles.openpistemap.org/landshaded/",
#		"osm-retina"="http://tile.geofabrik.de/osm_retina/",
#      'osm-rail' = 'http://a.tiles.openrailwaymap.org/standard/',
# rail is 512 insstead of 256 tiles
#			hill="http://www.toolserver.org/~cmarqu/hill/",
#	eu="http://alpha.map1.eu/tiles/",
# 'esri' = 'http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer/tile/',
# 'esri-grey' = 'http://services.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/',
# 'esri-transport'='http://server.arcgisonline.com/ArcGIS/rest/services/Reference/World_Transportation/MapServer/tile/',
# 'esri-topo' = 'http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/'
  
  
# http://server.arcgisonline.com/arcgis/rest/services/Ocean/World_Ocean_Base/MapServer/tile/2/1/1.jpg
  
  # language labels don't appear to be working
  languages = c("en","fr","de", "it","es","ru")
  toadd =	paste("http://a.www.toolserver.org/tiles/osm-labels-", languages,"/", sep="")
  names(toadd) = paste("osm-labels-", languages, sep="")
#	result = c(result, toadd)
  
  
  result = c(result, toadd)
  
  
  
  if(!missing(name)) {
    if(all(name %in% names(result), na.rm=TRUE)) {
      result = result[name]
    } else {
      result = name
    }
  }
  if(!missing(xyz))
    attributes(result)$tileNames = xyz	
  if(!missing(suffix))
    attributes(result)$suffix = suffix	
  
  result
  
}



if(FALSE){
 x = vect(as.matrix(expand.grid(seq(-5e6,-1e6,len=100), seq(-3e6,3e6,len=100))), crs='EPSG:3573')
 xLL = project(x, crsLL)
 zoom=4
 theExt = ext(-6e6,6e6,-6e6,6e6)
 xx = openmap(
  x=rast(theExt, res = (xmax(theExt)-xmin(theExt))/1200, crs=crs('EPSG:3031')),
  path="https://stamen-tiles-d.a.ssl.fastly.net/toner/",  
#  path="https://tiles.stadiamaps.com/styles/stamen_watercolor/", suffix='.jpg',
  zoom=2, verbose=TRUE, fact=2)

}


openmap = function(
  x, 
  zoom, 
  path="http://tile.openstreetmap.org/",
  maxTiles = 9,
  crs=terra::crs(x),
  buffer=0, fact=1,
  verbose=getOption('mapmiscVerbose'),
  cachePath=getOption('mapmiscCachePath')
) {



  verbose = max(c(0, verbose))
  

  NtestCols = 100

  
  if(!is.null(attributes(x)$ellipse) ) {
    # to do: check for ellipses
    # x is a crs object
    # get the ellipse
    crs = x
    toCrop = attributes(x)$ellipse
    x = attributes(x)$regionLL
  } else {
    toCrop = NULL
  }
  


  crsOut=crs

  crsIn = terra::crs(x)

  if(all(is.na(crsIn))) {
    if(is.vector(x)){
      crsIn=crsLL
    } else{
      crsIn = crs 
    }
  }
  

# get extent of output

# get output raster
  testRast = rast(terra::ext(x), res = (terra::xmax(x) - terra::xmin(x))/NtestCols, crs = crs(x))
  testPoints = vect(terra::xyFromCell(testRast, 1:terra::ncell(testRast)), crs=terra::crs(testRast))

  testPointsMerc = project(testPoints, crsMerc)
  testPointsOut = terra::project(testPoints, crsOut)
  outExtent= terra::ext(testPointsOut)

# buffer
  if(terra::is.lonlat(crsOut) & any(buffer > 90)) {
          # assume buffer is in km
          # transform to merc , buffer, transform back
          outExtentMerc = terra::extend(terra::ext(testPointsMerc), buffer)
          outPointsMerc = vect(matrix(as.vector(outExtentMerc), ncol=2), crsMerc)
          outPointsLL = project(outPointsMerc, crsOut)
          outExtent = terra::ext(outPointsLL)

    } else {
          outExtent = terra::extend(outExtent, buffer)
    }

  testRast = rast(outExtent, res = (terra::xmax(outExtent) - terra::xmin(outExtent))/NtestCols, crs = crsOut)
  testPoints = vect(terra::xyFromCell(testRast, 1:terra::ncell(testRast)), crs=terra::crs(testRast))
  testPointsMerc = project(testPoints, crsMerc)

  if(missing(zoom)) {


    # get zoom

    zoom = 0
    Ntiles = 0
    while(Ntiles <= maxTiles & zoom <= 18) {
      zoom = zoom + 1
      Ntilesm1 = Ntiles
      Ntiles = length(unique(terra::cellFromXY(.getRasterMerc(zoom), terra::crds(testPointsMerc))))
    }
    zoom = zoom - 1
    if(verbose) cat("zoom is ", zoom, ", ", Ntilesm1, "tiles\n")
  }



  # create out raster
  # find average area of pixels in downloaded tiles

    mercHere = .getRasterMerc(zoom)
   if(identical(crsOut, crsMerc)) {
        # output crs is mercator, return tiles as-is
        outraster = terra::crop(mercHere, testRast, snap='out')
        outraster = terra::disagg(outraster, 256)

    } else {
    theTable = as.data.frame(table(terra::cellFromXY(mercHere, terra::crds(testPointsMerc))))
    theTable$cell = as.numeric(as.character(theTable[,1]))
 
    mercHere = terra::crop(mercHere, 
      terra::ext(rep(terra::xyFromCell(mercHere, theTable[which.max(theTable$Freq), 'cell']), each=2) + 
        0.6*rep(terra::res(mercHere), each=2)*c(-1,1,-1,1)))
    # each tile is 256 x 256
    mercHere = terra::disagg(mercHere, 256)

    # width of the cell with the most test points in it
    cellWidthMerc = quantile(terra::values(terra::cellSize(mercHere, unit='m')), 0.5)


    areaRast = terra::cellSize(testRast, unit='m')
    cellWidthRast = quantile(unlist(terra::spatSample(areaRast,  size=200)), prob=0.5)


    areaRatio = cellWidthRast/cellWidthMerc

    newNumberOfCells = fact*NtestCols*sqrt(areaRatio)


    outraster = rast(outExtent, res = (terra::xmax(outExtent) - terra::xmin(outExtent))/newNumberOfCells, crs = crsOut)
  } # end not merc

# cache


  if(is.null(cachePath)) {
    cachePath = tempdir()
  }
  if(!nchar(cachePath)) {
    cachePath = '.'
  }
  
  alltiles = osmTiles()
  pathOrig = path
  pathisname = gsub("-", ".", pathOrig) %in% 
  gsub("-", ".", names(alltiles))
  path[pathisname] = alltiles[path[pathisname]]
  
  if(length(grep("[[:punct:]]$", path, invert=TRUE)))
    path[ grep("[[:punct:]]$", path, invert=TRUE)] =
  paste(path[ grep("[[:punct:]]$", path, invert=TRUE)], 
    "/", sep="")
  
  if(length(grep("^http[s]*://", path, invert=TRUE)))
    path[ grep("^http[s]*://", path, invert=TRUE)] = 
  paste("http://", 
    path[ grep("^http[s]*://", path, invert=TRUE)], sep="")
  names(path) = pathOrig
  




  
  Dpath = names(path)[1]
  Durl = path[1]

  if(verbose){
    cat(Dpath, '\n')
    cat(Durl, '\n')
  }

  if(length(grep(
    'nrcan\\.gc\\.ca|gov\\.bc\\.ca', Durl))
){
    suffix = ''
    tileNames = 'zyx'
  } else if(
    length(grep(
      '[.]arcgisonline[.]com', Durl
    ))) {
    suffix='.jpg'
    tileNames = 'zyx'
  } else if(
    length(grep(
      'stamen.watercolor', Durl
    ))) {
    suffix='.jpg'
    tileNames = 'zyx'
  } else if(
    length(grep(
      'heidelberg.de/tiles/(hybrid|adminb|roadsg|roads)/?$', 
      Durl)) |
    length(grep(
      '&$',Durl))
  ){
    tileNames = 'xyz='
    suffix = ''
  } else {
    suffix = '.png'
    tileNames = 'zxy'
  }

  if(length(attributes(pathOrig)$tileNames))
    tileNames = attributes(pathOrig)$tileNames
  if(length(attributes(pathOrig)$suffix))
    suffix = attributes(pathOrig)$suffix




  result = try(
    getTiles(outraster, 
      zoom=zoom,
      path=Durl,
      verbose=verbose,
      suffix=suffix,
      tileNames = tileNames,
      cachePath = cachePath),
    silent=!verbose)

  if(any(class(result)=="try-error")){
    message(paste(Durl, "not accessible"))
    # create an empty raster
    result = outraster
    attributes(result)$openmap = list(
      tiles=NA,
      message=result,
      path=path,
      pathOrig=pathOrig,
      zoom=zoom
    )
  } else {
    attributes(result)$openmap = list(
      path=path,
      pathOrig=pathOrig,
      zoom=zoom
    )
  }



  result
}