CA_brushes.qmd

---
title: "Tools of mapping"
bibliography: ZoteroDBase.bib
---
::: {.content-hidden}
# CA_Brushes
[[20220227120733]] 01 Welcome to my Zenkettle 
[[20230829120127]] My digital footprint  
[[20240113032207]] My public footprint

**In page links:**

**Related:**
[[gcontext.qmd]] gcontext.qmd

**To Ingest:**
https://r-graph-gallery.com/ggplot2-color.html
https://github.com/EmilHvitfeldt/r-color-palettes
https://gis.stackexchange.com/questions/25099/choosing-colour-ramp-to-use-for-elevation
https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
https://stackoverflow.com/questions/74905034/asymmetric-color-scale-with-colorramppalette-in-r
https://github.com/EmilHvitfeldt/r-color-palettes?tab=readme-ov-file
https://r-charts.com/color-palettes/
https://milospopovic.net/crisp-topography-map-with-r/

## Graveyard

Following my Bob Ross "Paint the landscape" metaphor, there are a few choice brushes I reach for when you ask me to paint a picture.  

Plotting aspatial data can be found [in there respective pages](PD.qmd) and it's a matter of file management and patience to programmatically stitch images together.  Therefore, this page documents to tools and techniques I use to paint my bread and butter.  

Leaflet
[tmap](https://tmap.geocompx.org/)
mapview
ggplot2
[tidyterra](https://dieghernan.github.io/tidyterra/index.html)
GIMP

Given the extreme I've gone to making my digital footprint as seamless as possible, it should come as no surprise I'm a fan of the programatic approach to map making offered by tools like tmaps, mapview, and ggplot2.  While each of these tools might all fall into the "chissle" category, having more than one chissle is critical to the master sculptor.  Leaflet is my favorite because the end user, with sufficient context, will almost naturally explore the data and make a far more meaningful narrative for themselves as they navigate that map than the forced perspective applied by the shitty artist (me) bakes into a static narrative.
[[20240925231616]] Junk Drawer

[Utility belt](20250922022711.md)
[fileio](fileio.qmd)  |  [data sources](datasources.qmd)  |  [templates](templates.qmd)  |  [color](color.qmd)  |  [Indexes](GI.qmd)  |  [Brushes](CA_brushes.qmd)  |  [Junk Drawer](junk_drawer.qmd)  |  [The Grill](FlatGrillDeck.qmd)

:::
```{r, warning = FALSE, echo = FALSE, results='asis'}
library(magrittr)
zenkettle_path <- "~/Dropbox/root/reference/zenkettle" 
path_to_basedata <- "~/Dropbox/root/database/hosted" 
path_to_OWP_hand_outputs <- "~/data/raw/lynker-spatial/FIM/OWP/hand_fim/hand_4_8_7_2"
revelio=TRUE
is_verbose=FALSE
zettlr_to_quarto_text <- function(file_path) {
  # file_path <- file.path("filename.md")
  if(!file.exists(file_path)) { return("I lost a file somewhere (ツ)_/¯") }
  markdown_content <- readLines(file_path) 
  yaml_header <- yaml::yaml.load(paste(markdown_content[(grep("^---$", markdown_content)[1] + 1):grep("^---$", markdown_content)[2] - 1], collapse = "\n"))
  if("private" %in% names(yaml_header)) {
    return("I didn't mean to link to this yet")
  } else if("startgrep" %in% names(yaml_header)) {
    return(cat(markdown_content[(grep(yaml_header$startgrep,x = markdown_content,fixed = TRUE)[[2]]+2):length(markdown_content)], sep = '\n'))
  } else if("starton" %in% names(yaml_header)) {
    return(cat(markdown_content[yaml_header$starton:length(markdown_content)], sep = '\n'))
  }
} # zettlr_to_quarto_text("filename.md")
```

![](./_bob_jim.png)

Following my Bob Ross "Paint the landscape" metaphor, there are a few choice brushes I reach for when you ask me to paint a picture.  

## Coloring

> See also: [data visualization](gcontext.qmd) & [color](color.qmd)

### [DEM](https://jimcoll.github.io/classes/datasources.html#category=Elevation)

Source Files from, for example, [3DEP](datasources_3dep.qmd).

`r if (!revelio) "::: {.content-hidden}"`

::: {.panel-tabset}

## QGIS

I'm fond of the `wkp/country/wiki-scotland` palette here.

![_0b239ec1c312e88716fa80073df97630.png](./_0b239ec1c312e88716fa80073df97630.png)

## R (through ggplot with [tidyterra](https://dieghernan.github.io/tidyterra/))

```{r}
surface <- terra::rast(file.path(path_to_basedata,"demo/HUC6_120901_dem.tif",fsep = .Platform$file.sep))

ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = surface) +
  tidyterra::scale_fill_whitebox_c(name = "Elevation (m)") +
  ggplot2::theme_void()
```

:::

Landscape processing to arrive at these: [Terrain Analysis](TA.qmd)

### Hillshade

::: {.panel-tabset}

## R

```{r}
surface <- terra::rast(file.path(path_to_basedata,"demo/HUC6_120901_dem.tif",fsep = .Platform$file.sep))

## DEM
ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = surface) +
  tidyterra::scale_fill_whitebox_c(name = "Elevation (m)") +
  ggplot2::theme_void()

## SLOPE
# slope <- terra::terrain(surface, "slope", unit = "degrees")
# slope_pal <- RColorBrewer::brewer.pal(9, "YlGnBu")
# ggplot2::ggplot() +
#   tidyterra::geom_spatraster(data = slope) +
#   ggplot2::scale_fill_gradientn(colors = slope_pal, na.value = NA, name = "Slope (°)") +
#   ggplot2::theme_void()

## Combined
slope <- terra::terrain(surface, "slope", unit="radians") 
aspect <- terra::terrain(surface, "aspect", unit="radians")
hillshade <- terra::shade(slope, aspect, 40, 270) 

ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = hillshade) +
  ggplot2::scale_fill_gradient(low = "black", high = "white", na.value = NA, name = "Hillshade") +
  ggnewscale::new_scale_fill() +
  tidyterra::geom_spatraster(data = surface, alpha = 0.5) +
  tidyterra::scale_fill_whitebox_c(name = "Elevation") +
  ggplot2::theme_void()
```

## QGIS

![](_todo.jpg)

## Google Earth Engine

```{javascript, echo = FALSE, eval = FALSE}
// https://www.youtube.com/watch?v=1ZkXZnw-tsc
// https://helpx.adobe.com/photoshop/using/blending-modes.html

//var elevation_data = ee.Image('CGIAR/SRTM90_V4').select('elevation');
var elevation_data = ee.Image('USGS/3DEP/10m').select('elevation');

function radians(img) {
  return img.toFloat().multiply(Math.PI).divide(180);
}
function hillshade(az, ze, slope, aspect) {
  var azimuth = radians(ee.Image(az));
  var zenith = radians(ee.Image(ze));
  return azimuth.subtract(aspect).cos()
    .multiply(slope.sin())
    .multiply(zenith.sin())
    .add(
      zenith.cos().multiply(slope.cos()));
}

var elevation_data_blur0 = elevation_data;
var elevation_data_blur7 = elevation_data.focalMean(7, "circle", "pixels", 1, null);
var elevation_data_blur15 = elevation_data.focalMean(15, "circle", "pixels", 1, null);
var elevation_data_blur30 = elevation_data.focalMean(30, "circle", "pixels", 1, null);

var slope_blur0 = ee.Terrain.slope(elevation_data_blur0);
var slope_blur7 = ee.Terrain.slope(elevation_data_blur7);
var slope_blur15 = ee.Terrain.slope(elevation_data_blur15);
var slope_blur30 = ee.Terrain.slope(elevation_data_blur30);

var aspect_blur0 = ee.Terrain.aspect(elevation_data_blur0);
var aspect_blur7 = ee.Terrain.aspect(elevation_data_blur7);
var aspect_blur15 = ee.Terrain.aspect(elevation_data_blur15);
var aspect_blur30 = ee.Terrain.aspect(elevation_data_blur30);

var hillshade_azimuth = 320; //270
var hillshade_elevation = 45;
var max_opacity_percent = 0.40;
var top_brush = 0.8;
var bottom_brush = 6;

var elevation_data_blur0_hs = ee.Terrain.hillshade(elevation_data_blur0,hillshade_azimuth,hillshade_elevation);
elevation_data_blur0_hs = elevation_data_blur0_hs.mask(
  (elevation_data_blur0_hs.subtract((1-max_opacity_percent)*255))
  .divide(255).clamp(0,1)
  );

var elevation_data_blur7_hs = ee.Terrain.hillshade(elevation_data_blur7,hillshade_azimuth,hillshade_elevation);
elevation_data_blur7_hs = elevation_data_blur7_hs.mask(
  (elevation_data_blur7_hs.subtract((1-max_opacity_percent)*255))
  .divide(255).clamp(0,1)
  );

var elevation_data_blur15_hs = ee.Terrain.hillshade(elevation_data_blur15,hillshade_azimuth,hillshade_elevation);
elevation_data_blur15_hs = elevation_data_blur15_hs.mask(
  (elevation_data_blur15_hs.subtract((1-max_opacity_percent)*255))
  .divide(255).clamp(0,1)
  );
  
var elevation_data_blur30_hs = ee.Terrain.hillshade(elevation_data_blur30,hillshade_azimuth,hillshade_elevation);
elevation_data_blur30_hs = elevation_data_blur30_hs.mask(
  (elevation_data_blur30_hs.subtract((1-max_opacity_percent)*255))
  .divide(255).clamp(0,1)
  );

Map.addLayer(ee.Image(1),{palette:['ffffff']});
Map.addLayer(elevation_data_blur0,{palette:['000000','000000']},'raw dem',false);
Map.addLayer(elevation_data_blur0_hs,{palette:['000000','000000']},'blur0 HS');
Map.addLayer(elevation_data_blur7_hs,{palette:['000000','000000']},'blur7 HS');
Map.addLayer(elevation_data_blur15_hs,{palette:['000000','000000']},'blur15 HS');
Map.addLayer(elevation_data_blur30_hs,{palette:['000000','000000']},'blur30 HS');

slope_blur0 = slope_blur0.mask(slope_blur0.divide(90).clamp(0,1));
Map.addLayer(slope_blur0,{palette:['000000','000000']},'slope0 HS');

var elevation_data_white_highlights_hs = ee.Terrain.hillshade(elevation_data_blur0,hillshade_azimuth,hillshade_elevation);
elevation_data_white_highlights_hs = elevation_data_white_highlights_hs.mask(
  elevation_data_white_highlights_hs.divide(255));
Map.addLayer(elevation_data_white_highlights_hs,{min:255*top_brush, max:255, palette:['ffffff','ffffff']},'slope highlights');

slope_blur7 = slope_blur7.mask(
  (slope_blur7.subtract((1-max_opacity_percent)*90))
  .divide(90).clamp(0,1));
Map.addLayer(slope_blur7,{palette:['000000','000000']},'slope7');
slope_blur15 = slope_blur15.mask(
  (slope_blur15.subtract((1-max_opacity_percent)*90))
  .divide(90).clamp(0,1));
Map.addLayer(slope_blur0,{palette:['000000','000000']},'slope15');
slope_blur30 = slope_blur30.mask(
  (slope_blur30.subtract((1-max_opacity_percent)*90))
  .divide(90).clamp(0,1));
Map.addLayer(slope_blur30,{palette:['000000','000000']},'slope30');

var slope_blur0_hs = ee.Terrain.hillshade(ee.Terrain.slope(elevation_data_blur0),hillshade_azimuth,hillshade_elevation);
var slope_blur0_hs_mask = ee.Image(slope_blur0_hs.divide(255)).interpolate(ee.List([0,0.25,0.75,1]),ee.List([1,0,0,1]),'mask');
slope_blur0_hs = slope_blur0_hs.updateMask(slope_blur0_hs_mask);
Map.addLayer(slope_blur0_hs,{palette:['000000','ffffff']},'slope_blur0_hs');

var slope_blur7_hs = ee.Terrain.hillshade(ee.Terrain.slope(elevation_data_blur7),hillshade_azimuth,hillshade_elevation);
var slope_blur7_hs_mask = ee.Image(slope_blur7_hs.divide(255)).interpolate(ee.List([0,0.25,0.75,1]),ee.List([1,0,0,1]),'mask');
slope_blur7_hs = slope_blur7_hs.updateMask(slope_blur7_hs_mask);
Map.addLayer(slope_blur7_hs,{palette:['000000','ffffff']},'slope_blur7_hs');

var low_shadows_mask = elevation_data_blur0.interpolate(ee.List([20,70]),ee.List([0.4,0]),'clamp');
Map.addLayer(elevation_data_blur0.updateMask(low_shadows_mask),{palette:['000000','ffffff']},'low lands');
```

![_8e7ee7efe11d80763ed7086757b34621.png](./_8e7ee7efe11d80763ed7086757b34621.png)

:::

### Aspect

::: {.panel-tabset}

## R

```{{r}}
pkgset <- '102702'
DEM <- glue::glue("~/data/raw/FIM4/inputs/dems/3dep_dems/10m_5070/20250320/HUC6_{substr(pkgset, 1, 6)}_dem.tif") |> terra::rast()

slope_rad <- terra::terrain(DEM, "slope", unit = "radians")
slope_pct <- tan(slope_rad) * 100
m <- c(0, 5, 10,
       5, 20, 20,
       20, 40, 30,
       40, max(terra::values(slope_pct),na.rm = T), 40)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
slope_rcl <- terra::classify(slope_pct, rclmat, include.lowest=TRUE)
asp <- terra::terrain(DEM, v = "aspect", unit = "degrees")
m <- c(0,22.5,1,
       22.5,67.5,2,
       67.5,112.5,3,
       112.5,157.5,4,
       157.5,202.5,5,
       202.5,247.5,6,
       247.5,292.5,7,
       292.5,337.5,8,
       337.5,max(terra::values(asp),na.rm = T),1)
rclmat <- matrix(m, ncol=3, byrow=TRUE)
asp_rcl <- terra::classify(asp, rclmat, include.lowest=TRUE)
aspect_slope <- asp_rcl + slope_rcl

aspect_slope_colors <- c(
  "11"="#E1E1E1", "12"="#E1E1E1", "13"="#E1E1E1", "14"="#E1E1E1", 
  "15"="#E1E1E1", "16"="#E1E1E1", "17"="#E1E1E1", "18"="#E1E1E1", "19"="#E1E1E1", "20"="#E1E1E1",
  "21"="#ADC0A8", "22"="#ABD0D1", "23"="#B0C8D6", "24"="#B5BDD6", "25"="#B7B1C9", "26"="#C9B1B9", "27"="#D6B8B0", "28"="#D4C4AA",
  "31"="#6DA368", "32"="#68B1B6", "33"="#74A1C2", "34"="#7F8DC2", "35"="#8674A6", "36"="#A67484", "37"="#C28674", "38"="#BBA168",
  "41"="#24821E", "42"="#1D8E94", "43"="#317AA8", "44"="#435AA8", "45"="#543180", "46"="#80314C", "47"="#A84C31", "48"="#9E7A1D"
)

levels_df <- data.frame(ID = as.numeric(names(aspect_slope_colors)), 
                        category = names(aspect_slope_colors))

# 2. Set the raster as categorical without copying the whole object
levels(aspect_slope) <- levels_df

# 2. Use as.factor() inside the aes() call to fix the discrete vs continuous error
ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = aspect_slope) + 
  ggplot2::scale_fill_manual(
    values = aspect_slope_colors, 
    na.value = "transparent",
    guide = "none" 
  ) +
  ggplot2::theme_void() +
  ggplot2::coord_sf(expand = FALSE)

```

## QGIS

![](_todo.jpg)

:::

### Flow Direction

::: {.panel-tabset}

## ggplot



## QGIS

:::

### Blended

::: {.panel-tabset}

## Google Earth Engine

```{javascript, echo = FALSE, eval = FALSE}
// Define a function to convert from degrees to radians.
function radians(img) {
  return img.toFloat().multiply(Math.PI).divide(180);
}

// Define a function to compute a hillshade from terrain data
// for the given sun azimuth and elevation.
function hillshade(az, ze, slope, aspect) {
  // Convert angles to radians.
  var azimuth = radians(ee.Image(az));
  var zenith = radians(ee.Image(ze));
  // Note that methods on images are needed to do the computation.
  // i.e. JavaScript operators (e.g. +, -, /, *) do not work on images.
  // The following implements:
  // Hillshade = cos(Azimuth - Aspect) * sin(Slope) * sin(Zenith) +
  //     cos(Zenith) * cos(Slope)
  return azimuth.subtract(aspect).cos()
    .multiply(slope.sin())
    .multiply(zenith.sin())
    .add(
      zenith.cos().multiply(slope.cos()));
}

var elevation = ee.Image("USGS/SRTMGL1_003");
// var slopeDegrees = ee.Algorithms.Terrain(elevation).select('slope').abs();
var slopeDegrees = ee.Terrain.slope(elevation).select('slope');
var slopePercent = ee.Image(slopeDegrees.tan()).multiply(100).abs();
var aspect = ee.Terrain.aspect(elevation).select('aspect');

var slopereclass = ee.Image(slopePercent)
  .where(slopePercent.gt(0).and(slopePercent.lte(5)), 10)
  .where(slopePercent.gt(5).and(slopePercent.lte(20)), 20)
  .where(slopePercent.gt(20).and(slopePercent.lte(40)), 30)
  .where(slopePercent.gt(40), 40);
var aspectreclass = ee.Image(aspect)
  .where(aspect.gt(0).and(aspect.lte(22.5)), 1)
  .where(aspect.gt(22.5).and(aspect.lte(67.5)), 2)
  .where(aspect.gt(67.5).and(aspect.lte(112.5)), 3)
  .where(aspect.gt(112.5).and(aspect.lte(157.5)), 4)
  .where(aspect.gt(157.5).and(aspect.lte(202.5)), 5)
  .where(aspect.gt(202.5).and(aspect.lte(247.5)), 6)
  .where(aspect.gt(247.5).and(aspect.lte(292.5)), 7)
  .where(aspect.gt(292.5).and(aspect.lte(337.5)), 8)
  .where(aspect.gt(337.5).and(aspect.lte(360)), 1);
var colorTerrain = aspectreclass.add(slopereclass);

// Map.addLayer(slopeDegrees,{},'slope (degrees)',false);
// Map.addLayer(slopePercent,{},'slope (percent)',false);
// Map.addLayer(slopereclass,{},'slope reclass',false);
// Map.addLayer(aspect,{},'aspect',false);
// Map.addLayer(aspectreclass,{},'aspect reclass',false);
Map.addLayer(colorTerrain,{
  min:19,
  max:48,
  palette:['#a1a1a1','#98b581','#72a890','#7c8ead','#8c75a0','#b47ba1',
  '#cb8b8f','#c5a58a','#bdbf89','#8dc458','#3dab71','#5078b6','#77479d',
  '#c04d9c','#e76f7a','#e2a66c','#d6db5e','#84d600','#84d600','#0068c0',
  '#6c00a3','#ca009c','#ff5568','#ffab47','#f4fa00']});
var x = 0.1;
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 1))),{color: "#2f6ecc", width: 0.002}); 
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 2))),{color: "#2f6ecc", width: 0.002}); 
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 3))),{color: "#2f6ecc", width: 0.002});
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 4))),{color: "#2f6ecc", width: 0.002}); 
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 5))),{color: "#2f6ecc", width: 0.002});
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 6))),{color: "#2f6ecc", width: 0.002}); 
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 7))),{color: "#2f6ecc", width: 0.002});
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 8))),{color: "#2f6ecc", width: 0.002}); 
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 9))),{color: "#2f6ecc", width: 0.002});
// Map.addLayer(ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 10))),{color: "#2f6ecc", width: 0.002});

var riv1 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 1))).style({color: "#2f6ecc", width: 1}); 
var riv2 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 2))).style({color: "#2f6ecc", width: 1.2}); 
var riv3 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 3))).style({color: "#2f6ecc", width: 1.4});
var riv4 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 4))).style({color: "#2f6ecc", width: 1.6}); 
var riv5 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 5))).style({color: "#2f6ecc", width: 1.8});
var riv6 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 6))).style({color: "#2f6ecc", width: 2}); 
var riv7 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 7))).style({color: "#2f6ecc", width: 2.2});
var riv8 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 8))).style({color: "#2f6ecc", width: 2.4}); 
var riv9 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 9))).style({color: "#2f6ecc", width: 2.6});
var riv10 = ee.FeatureCollection(Rivers.filter(ee.Filter.eq('RIV_ORD', 10))).style({color: "#2f6ecc", width: 2.8});
Map.addLayer(riv1); 
Map.addLayer(riv2); 
Map.addLayer(riv3); 
Map.addLayer(riv4); 
Map.addLayer(riv5); 
Map.addLayer(riv6); 
Map.addLayer(riv7); 
Map.addLayer(riv8); 
Map.addLayer(riv9); 
Map.addLayer(riv10); 

Map.addLayer(BIA.union().style({color: '#8a5b19', fillColor: '#8a5b1935',width:2}));
Map.addLayer(Aqu.union().style({color: '#1f4bab', fillColor: '#1f4bab45',width:4}));
Map.addLayer(Cities.style({color: '#ecf230'}));
Map.addLayer(SDWells.style({color: '#22e352'}));
Map.addLayer(Uran.style({color: '#cf3221'}));
var State = ee.FeatureCollection("TIGER/2018/States");
var filter = State.filter(ee.Filter.inList("NAME",['South Dakota']).not());
Map.addLayer(filter.geometry().dissolve(), {opacity:0.3}, 'State');
var legend = ui.Panel({
  style: {
    position: 'bottom-right',
    padding: '8px 15px'
  }
});

// Create legend title
var legendTitle = ui.Label({
  value: 'Uranium mines, Wells, Waterways, Terrain, and Populations in the Black Hills of South Dakota',
  style: {
    fontWeight: 'bold',
    fontSize: '18px',
    margin: '0 0 4px 0',
    padding: '0',
    maxWidth:'390px'
  }
});
 
// Add the title to the panel
legend.add(legendTitle);
 
// Creates and styles 1 row of the legend.
var makeRow = function(color, name) {
 
      // Create the label that is actually the colored box.
      var colorBox = ui.Label({
        style: {
          backgroundColor: '#' + color,
          // Use padding to give the box height and width.
          padding: '8px',
          margin: '0 0 4px 0'
        }
      });
 
      // Create the label filled with the description text.
      var description = ui.Label({
        value: name,
        style: {margin: '0 0 4px 6px'}
      });
 
      // return the panel
      return ui.Panel({
        widgets: [colorBox, description],
        layout: ui.Panel.Layout.Flow('horizontal')
      });
};
//  Palette with the colors
var palette =['cf3221', '22e352','ecf230','1f4bab','8a5b19'];
 
// name of the legend
var names = ['Uranium mines (Surficial)','Water wells','Major cities','Inyan Kara Group aquifers','American Indian Land Area Representation'];
 
// Add color and and names
for (var i = 0; i < 5; i++) {
  legend.add(makeRow(palette[i], names[i]));
}  
var legendSources = ui.Label({
  value: 'Sources: The BIA, USGS, SRTM, and the South Dakota Department of Agriculture & Natural Resources',
  style: {
    fontWeight: 'normal',
    fontSize: '12px',
    margin: '0 0 4px 0',
    padding: '0',
    maxWidth:'390px'
  }
});

// add legend to map (alternatively you can also print the legend to the console)



var logo = ee.Image('users/JamesMColl/LargeSlopeAspect').visualize({
    bands:  ['b1', 'b2', 'b3'],
    min: 0,
    max: 255
    });

var thumb = ui.Thumbnail({
    image: logo,
    params: {
        dimensions: '398x382',
        format: 'png'
        },
    style: {height: '398px', width: '382px',padding :'0'}
    });
legend.add(thumb);
Map.add(legend);
// var toolPanel = ui.Panel(thumb, 'flow', {width: '300px'});
// ui.root.widgets().add(toolPanel);
legend.add(legendSources);
Map.setCenter(-103.061848427201,43.811618572829055, 9);
```

![_b7fc57afd097034b03732421b557206b.png](./_b7fc57afd097034b03732421b557206b.png)

:::

### HAND

> See [FIM](FIM_primer.qmd) for more details.

::: {.panel-tabset}

## Leaflet

```{r, warning = FALSE, echo = FALSE, eval = TRUE}
unit_to_map <- "12090301"

# Generic inputs
HAND_path <- file.path(path_to_OWP_hand_outputs,unit_to_map,"branches","0","rem_zeroed_masked_0.tif",fsep = .Platform$file.sep)
HAND <- terra::rast(HAND_path)
HAND_stars <- stars::read_stars(HAND_path)

HAND_values <- terra::values(HAND)
pal <- leaflet::colorNumeric(RColorBrewer::brewer.pal(5, "Blues"), c(min(HAND_values,na.rm = TRUE),max(HAND_values,na.rm = TRUE)),na.color = "transparent")
aoi_bbox <- sf::st_bbox(sf::st_transform(sf::st_as_sfc(sf::st_bbox(HAND)),sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()

map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leaflet::addRasterImage(HAND, colors = pal, opacity = 1,maxBytes=50*1024*1024) %>%
  leafem::addGeoRaster(HAND_stars, colorOptions = leafem::colorOptions(palette = RColorBrewer::brewer.pal(5, "Blues"), domain=c(min(HAND_values,na.rm = TRUE),max(HAND_values,na.rm = TRUE)), na.color = "transparent"), opacity = 1,resolution = 96) %>%
  leaflet::fitBounds(aoi_bbox[3],aoi_bbox[2],aoi_bbox[1],aoi_bbox[4])
map1
```

## ggplot

```{r, warning = FALSE, echo = FALSE, eval = TRUE}
unit_id <- "16060002"
branch_id <- "0"

# Older versions of FIM are in meters, not mm as V4872+ are
rem_version_surface_to_meter_conversion_factor = 1000
# rem_version_surface_to_meter_conversion_factor = 1

rem_file <- file.path(path_to_OWP_hand_outputs,unit_id,"branches",branch_id,glue::glue("rem_zeroed_masked_{branch_id}.tif"),fsep = .Platform$file.sep)
HAND_surface <- terra::rast(rem_file)
mapext <- sf::st_bbox(terra::project(HAND_surface,"EPSG: 3857")) |> sf::st_as_sfc() %>% sf::st_buffer(50)

Branch_HAND_map <- basemaps::basemap_ggplot(mapext,map_service = "osm", map_type = "topographic",verbose = FALSE) +
  ggnewscale::new_scale_fill() +
  tidyterra::geom_spatraster(data = terra::project(HAND_surface,"EPSG: 3857")) +
  tidyterra::scale_fill_whitebox_c(palette = "deep", name = "Elevation (m)") +
  ggplot2::theme_void()
Branch_HAND_map
```

## QGIS

::: 

### View setting

::: {.panel-tabset}

## Leaflet

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

```

## Leaflet zoom control

```{r, warning = FALSE, echo = FALSE, eval = FALSE}
leaflet::leaflet() %>%
  leaflet::addProviderTiles(providers$CartoDB.Positron) %>%
  leaflet::addCircleMarkers(lng=quakes$long,
                   lat=quakes$lat,
                   col="blue",
                   radius=3,
                   stroke=FALSE,
                   fillOpacity = 0.7,
                   group = "Quake Points") %>%                       
  leaflet::addPolygons(data= oceania,
              col="red") %>%                        
  leaflet::groupOptions("Quake Points", zoomLevels = 15:20)      
```

## Mapview

## ggplot

:::

### States (cookiecutter plot)

> Also see [the TIGER data page](datasources_tiger.qmd)

::: {.panel-tabset}

## CONUS

```{r, warning = FALSE, echo = FALSE, eval = FALSE}
states <- sf::st_read(file.path(path_to_basedata,"admin/tiger/TIGER2022/STATE/tl_2022_us_state/tl_2022_us_state.fgb",fsep = .Platform$file.sep),
                      query = glue::glue("SELECT * FROM tl_2022_us_state WHERE GEOID < '57%' AND GEOID NOT LIKE '02%' AND GEOID NOT LIKE '15%'"),
                      quiet=!is_verbose)
```

::: {.panel-tabset}

## Leaflet

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

```

## Mapview

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])
```

:::

## NA

## Leaflet

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])

```

## Mapview

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])
```

:::

## Exporting


::: {.panel-tabset}

## ggplot

[ [[20250916093105]] [ggplot2 export helper](20250916093105.md) ]{.content-hidden}
```{r, warning = FALSE, echo = FALSE,results='asis'}
zettlr_to_quarto_text("20250916093105.md")
```

## Leaflet

[ [[20250916093246]] [leaflet export helper](20250916093246.md) ]{.content-hidden}
```{r, warning = FALSE, echo = FALSE,results='asis'}
zettlr_to_quarto_text("20250916093246.md")
```

## Mapview

```{r, warning = FALSE, echo = FALSE, eval = FALSE}

 %>%
  leaflet::fitBounds(hf_bounds[3],hf_bounds[2],hf_bounds[1],hf_bounds[4])
```

:::

## tidyterra

```{r}
#| warning: false
#| echo: false
#| eval: false

# Pilfered from https://dieghernan.github.io/202210_tidyterra-hillshade/
basedir <- "~/data/raw/nohrsc/owp_files/nwm/nwm_parameters/NWM_coastal_parameters.tar/coastal_parameters/coastal/atlgulf/"
my_map_window

bbox <- sf::st_bbox(sf::st_transform(my_map_window,sf::st_crs('EPSG:4326'))) %>% unlist() %>% unname()
dem_meters <- terra::rast(dem_path)
dem_meters_values = terra::values(dem_meters)
range =  max(dem_meters_values,na.rm = TRUE) - min(dem_meters_values,na.rm = TRUE)

slope <- terra::terrain(dem_meters, "slope", unit = "radians")
aspect <- terra::terrain(dem_meters, "aspect", unit = "radians")
hill <- terra::shade(slope, aspect, 30, 270)
names(hill) <- "shades"  # normalize names
pal_greys <- grDevices::hcl.colors(1000, "Grays")

tidyterra::autoplot(dem_meters) +
  ggplot2::theme_minimal()
ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = hill) +
  ggplot2::scale_fill_gradientn(colors = pal_greys, na.value = NA)

index <- hill %>%
  dplyr::mutate(index_col = scales::rescale(shades, to = c(1, length(pal_greys)))) %>%
  dplyr::mutate(index_col = round(index_col)) %>%
  dplyr::pull(index_col)
vector_cols <- pal_greys[index]
hill_plot <- ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = hill, fill = vector_cols, maxcell = Inf,alpha = 1)

hill_plot
grad_hypso <- tidyterra::hypso.colors2(10, "dem_poster")

tidyterra::autoplot(dem_meters) +
  ggplot2::scale_fill_gradientn(colours = grad_hypso, na.value = NA)
# Try some options, but we need to be aware of the values of our raster

r_limits <- terra::minmax(dem_meters) %>% as.vector()
r_limits <- c(-5, 140) # Rounded to lower and upper 500
# r_limits <- pmax(r_limits, 0) # And making min value to 0.

# Now lets have some fun with scales from tidyterra

elevt_test <- ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = dem_meters)

# Create a helper function

# plot_pal_test <- function(pal) {
#   elevt_test +
#     scale_fill_hypso_tint_c(
#       limits = r_limits,
#       palette = pal
#     ) +
#     ggtitle(pal) +
#     theme_minimal()
# }
# 
# plot_pal_test("etopo1_hypso")
# plot_pal_test("dem_poster")
# plot_pal_test("spain")
# plot_pal_test("pakistan")
# plot_pal_test("utah_1")
# plot_pal_test("wiki-2.0_hypso")

base_plot <- hill_plot +
  tidyterra::geom_spatraster(data = dem_meters) +
  tidyterra::scale_fill_hypso_tint_c(limits = r_limits,palette = "dem_poster",alpha = 0.4,labels = scales::label_comma(),
    # For the legend I use custom breaks
    breaks = c(-5,10,20,40,60,80,100,120,140)
  )

base_plot

base_text_size <- 24
ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = hill, fill = vector_cols, maxcell = Inf,alpha = 1) +
  tidyterra::geom_spatraster(data = dem_meters) +
  tidyterra::scale_fill_hypso_tint_c(limits = r_limits,palette = "dem_poster",alpha = 0.4,labels = scales::label_comma(),
                                     breaks = c(-5,10,20,40,60,80,100,120,140)
  ) +
  ggplot2::guides(fill = ggplot2::guide_legend(
    title = "   m",
    direction = "horizontal",
    nrow = 1,
    keywidth = 1.75,
    keyheight = 0.5,
    label.position = "bottom",
    title.position = "right",
    override.aes = list(alpha = 1)
  )) +
  ggplot2::labs(
    title = "Elevation from USGS 3DEP",
    subtitle = "Domain context",
    caption = paste0(
      "Data: ClimateR and nwmTools"
    )
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    plot.background = ggplot2::element_rect("grey97", colour = NA),
    plot.margin = ggplot2::margin(20, 20, 20, 20),
    plot.caption = ggplot2::element_text(size = base_text_size * 0.5),
    plot.title = ggplot2::element_text(face = "bold", size = base_text_size * 1.4),
    plot.subtitle = ggplot2::element_text(
      margin = ggplot2::margin(b = 10),
      size = base_text_size
    ),
    axis.text = ggplot2::element_text(size = base_text_size * 0.7),
    legend.position = "bottom",
    legend.title = ggplot2::element_text(size = base_text_size * 0.8),
    legend.text = ggplot2::element_text(size = base_text_size * 0.8),
    legend.key = ggplot2::element_rect("grey50"),
    legend.spacing.x = ggplot2::unit(0, "pt")
  )
```

```{r}
#| warning: false
#| echo: false
#| eval: false

google_aoi <- c(44.28,-68.3)
distance_to_buffer = .6*10000

my_map_window <- sf::st_bbox(sf::st_buffer(sf::st_sfc(sf::st_point(c(google_aoi[2],google_aoi[1])),crs = sf::st_crs('EPSG: 4326')),units::as_units(distance_to_buffer,'meter'))) |> sf::st_as_sfc() |> sf::st_as_sf()
mapview::mapview(my_map_window)

NOMADS_from_nwmTools <- nwmTools::get_nomads_urls(
  output = 'total_water',
  domain = 'atlgulf',
  num = 18,
  version = 'v3.0',
  outdir = file.path("~/data/temp/fun/",fsep = .Platform$file.sep)) %>%
  nwmTools::get_timeseries(
    index_id = c('SCHISM_hgrid_node_x', 'SCHISM_hgrid_node_y'),
    varname  = 'elevation'
  )
NOMADS_from_nwmTools_sf <- sf::st_as_sf(NOMADS_from_nwmTools, coords = c("SCHISM_hgrid_node_x", "SCHISM_hgrid_node_y")) |>
  sf::st_set_crs(sf::st_crs('EPSG:6349')) |>
  sf::st_transform(sf::st_crs('EPSG: 4326'))
target_nodes <- NOMADS_from_nwmTools_sf[my_map_window,]
mapview::mapview(target_nodes)

# or
dap.catalog =  arrow::read_parquet("https://github.com/mikejohnson51/climateR-catalogs/raw/refs/heads/main/docs/catalog.parquet")
dem_catalog = dplyr::filter(dap.catalog,id == "USGS 3DEP",asset=="10m CONUS DEM")
DEM <- opendap.catalog::dap(URL = dem_catalog$URL, AOI = my_map_window)

target_nodes_z <- terra::extract(DEM,sf::st_transform(target_nodes,sf::st_crs(DEM)))
sf::st_drop_geometry(target_nodes)

for(i in 1:ncol(sf::st_drop_geometry(target_nodes))) {
  target_nodes[,i] <- target_nodes[,i] - target_nodes_z$USGS_Seamless_DEM_13
}


low_ocean_point <- do.call(pmin, c(sf::st_drop_geometry(target_nodes), list(na.rm=TRUE))) |> min(na.rm = TRUE)
high_ocean_point <-do.call(pmax, c(sf::st_drop_geometry(target_nodes), list(na.rm=TRUE))) |> max(na.rm = TRUE)

  
## Create hillshade effect
slope <- terra::terrain(DEM, "slope", unit = "radians")
aspect <- terra::terrain(DEM, "aspect", unit = "radians")
hill <- terra::shade(slope, aspect, 30, 270)
names(hill) <- "shades"

# Hillshading, but we need a palette
pal_greys <- hcl.colors(1000, "Grays")
grad <- tidyterra::hypso.colors(10, "dem_poster")

# Try some options, but we need to be aware of the values of our raster
r_limits <- minmax(DEM) %>% as.vector()
r_limits <- c(floor(r_limits[1] / 500), ceiling(r_limits[2] / 500)) * 500

# And making min value to 0.
r_limits <- pmax(r_limits, 0)

# Compare
minmax(r) %>% as.vector()
#> [1]    0 2481
r_limits
#> [1]    0 2500


# Now lets have some fun with scales from tidyterra

elevt_test <- ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = DEM)

# Create a helper function

plot_pal_test <- function(pal) {
  elevt_test +
    tidyterra::scale_fill_hypso_tint_c(
      limits = r_limits,
      palette = pal
    ) +
    ggtitle(pal) +
    theme_minimal()
}

plot_pal_test("etopo1_hypso")
plot_pal_test("dem_poster")
plot_pal_test("spain")
plot_pal_test("pakistan")
plot_pal_test("utah_1")
plot_pal_test("wiki-2.0_hypso")
```

## Leaflet

```{r}
#| warning: false
#| echo: false
#| eval: false

pal <- leaflet::colorNumeric(c("forestgreen","yellow","tan","brown"), dem_meters_values,na.color = "transparent")
map1 <- leaflet::leaflet(options = leaflet::leafletOptions(preferCanvas = TRUE)) %>%
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  # leafem::addFeatures(data = sf::st_transform(levelpaths,sf::st_crs('EPSG:4326')),stroke = TRUE,opacity = 1,fillOpacity = 1,weight = 1,color="#0327cd") %>% 
  leafem::addFeatures(data = sf::st_transform(my_map_window,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 0,weight = 2,color="#000000") %>% 
  # leafem::addFeatures(data = sf::st_transform(nearby_hucs,sf::st_crs('EPSG:4326')),stroke=TRUE,fill=FALSE,opacity = 1,fillOpacity = 1,weight = 1,color="#000000") %>% 
  # leafem::addGeoRaster(dem_meters,colorOptions = leafem:::colorOptions(palette = pal,
  #                                                                      # breaks = as.numeric(c(0:xx$max_ffreq)),
  #                                                                      domain = c(dem_meters_values_min,dem_meters_values_max))) %>%
  # leafem::addGeoRaster(dem_meters,color= pal) %>%
  leaflet::addRasterImage(dem_meters, colors = pal, opacity = 1,maxBytes=42*1024*1024) %>%
  leaflet::fitBounds(bbox[3],bbox[2],bbox[1],bbox[4])
# addLegend(pal = pal, values = terra::values(dem_meters), title = "Elevation data for Santander (mts)")
# map1
mapview::mapshot2(map1,file=file.path(trash_path,glue::glue("test.png"),fsep = .Platform$file.sep))
magick::image_write(
  magick::image_read(
    cropcircles::crop_circle(file.path(trash_path,glue::glue("test.png"),fsep = .Platform$file.sep), border_size = 4,just = "center")
  ),file.path(trash_path,glue::glue("test_cropped.png"),fsep = .Platform$file.sep))
unlink(file.path(trash_path,"*",fsep = .Platform$file.sep))
```

## Cartography with digital tools

If you've spent any time in this domain, you might find it becomes pretty easy to guess what tool was used to make a graph/map.  Tools like matplotlib or ggplot are extensible but the styles tend to be pretty cookie-cutter.  The following are my own cookie cutters that (should hopefully?) make slapping a map together using a standardized design language _easier_.

## Making a map

#### The most basic map

```{r}
#| echo: false
#| eval: false
leaflet::leaflet(options = leafletOptions(preferCanvas = TRUE)) %>% 
      leaflet::addProviderTiles("OpenStreetMap",group = "OpenStreetMap") %>%
      leaflet::addProviderTiles("Stamen.Toner",group = "Stamen.Toner") %>%
      leaflet::addProviderTiles("Stamen.Terrain",group = "Stamen.Terrain") %>%
      leaflet::addProviderTiles("Esri.WorldStreetMap",group = "Esri.WorldStreetMap") %>%
      leaflet::addProviderTiles("Wikimedia",group = "Wikimedia") %>%
      leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
      leaflet::addProviderTiles("Esri.WorldImagery",group = "Esri.WorldImagery") %>%
      leafgl::addGlPoints(unprocessed_nodes,fillColor = 'red',group = "pts") %>%
      leafgl::addGlPoints(processed_nodes,fillColor = 'green', group = "pts") %>%
      leafem::addFeatures(ahulls,opacity = 1,fillOpacity = 0,weight = 2,color = 'black', group = "huc") %>% 
      # leafgl::addGlPolygons(ahulls,fillOpacity = 0,stroke = TRUE,color = 'black',popup="huc10", group = "huc") %>%
      leaflet::addLegend("bottomright",colors = c("green","red"),
                         labels = c(paste0("SCHISM nodes: n=",processed_nodes_count), paste0("outside of huc bounds: n=",unprocessed_nodes_count)),
                         title = "SCHISM Domain (thinned)",opacity = 1) %>%
      leaflet::addLayersControl(
        baseGroups = c(
          "OpenStreetMap", "Stamen.Toner",
          "Stamen.Terrain", "Esri.WorldStreetMap",
          "Wikimedia", "CartoDB.Positron", "Esri.WorldImagery"
        ),
        position = "topleft",
        overlayGroups = c("pts","huc"))
```

#### A little more 

```{r}
#| echo: false
#| eval: false

my_map_window <- AOI::geocode("Fort Collins",bbox = TRUE)

## ==========
mapview::mapview(my_map_window)
## ==========
# Dark, positron, voyager, mapnik, 
ggplot2::ggplot() +
  basemapR::base_map(bbox = sf::st_bbox(my_map_window), basemap = 'mapnik', increase_zoom = 2, nolabels=FALSE) +
  ggplot2::geom_sf(data = my_map_window, fill = NA)
## ==========
bm <- basemaps::basemap_terra(ext = sf::st_bbox(sf::st_transform(my_map_window,sf::st_crs("EPSG:3857"))),map_service = "esri", map_type = "world_imagery", verbose = FALSE) %>%
  as("SpatRaster")
ggplot2::ggplot() +
  tidyterra::geom_spatraster_rgb(data = bm) +
  ggplot2::geom_sf(data = my_map_window, fill = NA)
## ==========
leaflet::leaflet() %>% 
  leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
  leafem::addFeatures(my_map_window,opacity = 1,fillOpacity = 0,weight = 2,color = 'black')
```

### Adding vector layers

```{r}
#| echo: false
#| eval: false
leaflet::leaflet() %>%
    leaflet::addProviderTiles(provider = providers$Stamen.TonerLabels) %>%
    leafgl::addGlPolygons(sf::st_transform(sf::st_cast(features_to_process,"POLYGON"),sf::st_crs("EPSG:4326")),fillOpacity = 0,stroke = TRUE,color = 'black', group = "huc") %>%
    leafgl::addGlPoints(sf::st_transform(valid_nodes,sf::st_crs("EPSG:4326")),fillColor = 'green',group = "pts") %>%
    leafgl::addGlPolygons(sf::st_transform(sf::st_cast(feature,"POLYGON"),sf::st_crs("EPSG:4326")),fillOpacity = 0,stroke = TRUE,color = 'black', group = "feat") %>%
    addLayersControl(
      overlayGroups = c("pts","huc","feat"),
      options = layersControlOptions(collapsed = FALSE)
    ) 
  

features_to_process <- features_still_needed[schsim_nodes,]

  leafgl::addGlPolygons(sf::st_transform(shcism_ahull_poly,sf::st_crs("EPSG:4326")),fillOpacity = 0,stroke = TRUE,color = 'black', group = "huc") %>%
    leafgl::addGlPolygons(sf::st_transform(feature_ahull,sf::st_crs("EPSG:4326")),fillOpacity = 0,stroke = TRUE,color = 'red', group = "huc")
```

### Adding raster layers

That's exactly what COG is; it has overviews built in and acts as a tile cache. COG has grown wildly in the last few years. COG use _range requests_ to extract part of a raster file out of cloud storage, this includes those small overviews that act as a tile cache. The above map uses a 1.5GB National Land Cover raster and the associated color table. Only overviews or small sections are used to draw the map layer depending on zoom and extent. The entire file is never downloaded. The initial load of the above map downloads < 1M of raster data.

```{r}
#| echo: false
#| eval: false
library("sf")
library("rgdal")
library("gdalUtilities")
library("shiny")
library("raster")
library("terra")
library("stars")
library("leaflet")
library("r-spatial/leafem")
library("r-spatial/leafgl")
library("mapview")

# Works (draws in map)
mapview::mapview(raster::raster('C:/Users/James.Coll/Desktop/data/102500/COG102500catchmask_p.tif'))   
# Does not work
mapview::mapview(stars::read_stars('C:/Users/James.Coll/Desktop/data/102500/COG102500catchmask_p.tif'))    
# How to use local COG here???
leaflet() %>% 
  addTiles(group = "osm") %>% 
  #addProviderTiles("Esri.WorldImagery", group = "esri") %>%
  addMapPane("cog", zIndex = 500) %>%
  setView(lng = 11.35, lat = 46.5, zoom = 12) %>%
  leafem:::addCOG(
    url = raster_lcog
    , group = "corine"
    , opacity = 0.8
    , options = list(pane = "cog")
    , resolution = 128
    , autozoom = FALSE
  ) 

## A working example of above:
hosturl <- 'http://oin-hotosm.s3.amazonaws.com/56f9b5a963ebf4bc00074e70/0/56f9c2d42b67227a79b4faec.tif'
leaflet() %>%
  addTiles() %>%
  leafem:::addCOG(
    url = hosturl
    , opacity = 1
    , options = list(pane = "cog")
    , resolution = 128
    , autozoom = FALSE
  )

url <- 'http://oin-hotosm.s3.amazonaws.com/56f9b5a963ebf4bc00074e70/0/56f9c2d42b67227a79b4faec.tif'
leaflet() %>% 
  addTiles(group = "osm") %>% 
  #addProviderTiles("Esri.WorldImagery", group = "esri") %>%
  addMapPane("cog", zIndex = 500) %>%
  setView(lng = 168.5, lat = -17.5, zoom = 12) %>%
     leafem:::addCOG(
         url = url
         , group = "corine"
         , opacity = 0.8
         , options = list(pane = "cog")
         , resolution = 128
         , autozoom = FALSE
       ) 
```

### Adding map flair

```{r}
#| echo: false
#| eval: false
leaflet() %>%
  addTiles() %>%
  leafem:::addFgb(
    url = url
    # , group = "counties"
    # , label = "NAME"
    , popup = TRUE
    , fill = TRUE
    , fillColor = "blue"
    , fillOpacity = 0.6
    , color = "black"
    , weight = 1
  ) %>%
  addLayersControl(overlayGroups = c("counties")) %>%
  addMouseCoordinates() %>%
  setView(lng = -105.644, lat = 51.618, zoom = 3)
```

### Writing her out

```{r}
#| echo: false
#| eval: false
unlink(file.path(output_folder_name,"SCHISM_based_FIM_domain_files"),recursive=TRUE)
    unlink(file.path(output_folder_name,"SCHISM_based_FIM_domain.html"))
    # m = leaflet(height=1000, width=1000) %>% 
    m = leaflet::leaflet(options = leafletOptions(preferCanvas = TRUE)) %>% 
      leaflet::addProviderTiles("OpenStreetMap",group = "OpenStreetMap") %>%
      leaflet::addProviderTiles("Stamen.Toner",group = "Stamen.Toner") %>%
      leaflet::addProviderTiles("Stamen.Terrain",group = "Stamen.Terrain") %>%
      leaflet::addProviderTiles("Esri.WorldStreetMap",group = "Esri.WorldStreetMap") %>%
      leaflet::addProviderTiles("Wikimedia",group = "Wikimedia") %>%
      leaflet::addProviderTiles("CartoDB.Positron",group = "CartoDB.Positron") %>%
      leaflet::addProviderTiles("Esri.WorldImagery",group = "Esri.WorldImagery") %>%
      leafgl::addGlPoints(unprocessed_nodes,fillColor = 'red',group = "pts") %>%
      leafgl::addGlPoints(processed_nodes,fillColor = 'green', group = "pts") %>%
      leafem::addFeatures(ahulls,opacity = 1,fillOpacity = 0,weight = 2,color = 'black', group = "huc") %>% 
      # leafgl::addGlPolygons(ahulls,fillOpacity = 0,stroke = TRUE,color = 'black',popup="huc10", group = "huc") %>%
      leaflet::addLegend("bottomright",colors = c("green","red"),
                         labels = c(paste0("SCHISM nodes: n=",processed_nodes_count), paste0("outside of huc bounds: n=",unprocessed_nodes_count)),
                         title = "SCHISM Domain (thinned)",opacity = 1) %>%
      leaflet::addLayersControl(
        baseGroups = c(
          "OpenStreetMap", "Stamen.Toner",
          "Stamen.Terrain", "Esri.WorldStreetMap",
          "Wikimedia", "CartoDB.Positron", "Esri.WorldImagery"
        ),
        position = "topleft",
        overlayGroups = c("pts","huc"))
    
mapview::mapshot(m,file.path(output_folder_name,'SCHISM_based_FIM_domain.html'))
```

### Sharing

Sometimes Google thinks I'm shipping people viruses when trying to email HTML, if you zip it those tend to do a little better with the email filters.

## Making a dashboard

![](_todo_moose.jpg)

## Resources:

https://tomjenkins.netlify.app/tutorials/r-leaflet-introduction/

## TMAP

https://rpubs.com/ials2un/dem_analysis

## USGS Topo Maps

From the [USGS map store](https://store.usgs.gov/map-locator).

`r if (!revelio) ":::"`

#### Interactive

* [Ancient Earth Globe](https://dinosaurpictures.org/ancient-earth#750): What the earth looked like _n_ years ago

### Other resources

https://nationalequityatlas.org/lab/get-started-data-visualization
https://handsondataviz.org/
https://blog.mapbox.com/7-best-practices-for-mapping-a-pandemic-9f203576a132
https://rpubs.com/KingaHill/947744
https://environmentalinformatics-marburg.github.io/mapview/popups/html/popups.html
https://flatgeobuf.org/examples/leaflet/filtered.html
https://stackoverflow.com/questions/50111566/applying-leaflet-map-bounds-to-filter-data-within-shiny
https://plugins.qgis.org/planet/user/29/tag/uncategorized/#:~:text=From%20the%20Print%20Layout%20toolbar,to%20a%20circular%20shape!).

## Explanation

### About Brutalist Cartography

![fcc64f58866e4c4e53ea54437b3f21b2.png](./_fcc64f58866e4c4e53ea54437b3f21b2.png)

Brutalist Cartography is a term of my own making which I use to describe the way I prepare the QGIS window for a screenshot.  Even if you've mastered the toolset and have the eye for aesthetics, good cartography can take quite a while to assemble and present.  That artistic flair is usually appreciated by all, but can slow the iteration process down quite substantially and is secondary to the main task of that map (for most of _my_ use cases: see [[20240508224456]] Defining roles and goals]), which is to place a tractable example of a theoretical state of the data in front of someone else with the goal of exampling problems, placing extra geographic context on a piece of data, or simply adding a better looking (in my opinion) screenshot.

To accomplish the Brutalist cartography [Aesthetic](https://www.youtube.com/watch?v=fjiQX74gqpk&ab_channel=Windows96), one simply turns off most of the QGIS user interface while more prominently accentuation an organized legend, an overview window, and some of the more common map decoration items like the corner attribution image and the scale bar.  

As a helpful note to future me, the panels that are typically turned on include:

Panels:

- [ ] Browser
- [ ] Layers

Toolbars:

- [ ] Annotations
- [ ] Attributes
- [ ] Data Source Manager
- [ ] Digitizing
- [ ] Help
- [ ] Label
- [ ] Manage layers
- [ ] Map Navigation
- [ ] Mesh digitization
- [ ] Plugin
- [ ] Project
- [ ] Selection
- [ ] Vector
- [ ] Web

#### The Final extra Step

Go to settings and turn font size to something larger like 24, close and reopen.  The downside is that you have to close and reopen to change this back again, and unfortunately this makes many of the menus unusable since the monitor resolution remains fixed.

## Old holdover

Although we could spend the day dissecting all the minute cartographic and stylistic errors the authors may have made in there initial map, the larger issue is the disingenuous presentation of the data.  The map in the initial assessment, shown below and recreated for demonstration purposes, was the map provided for the impact assessment.

![SLT_EJ_upscaled.png](./_SLT_EJ_upscaled.png)

![image (1).png](./_lfkimpact.png)

The assessment makes it appear that only 1 census block of predominantly minority populations are affected, when looking at the breakdown of impacted population as a total of the study domain we see those are far closer to a 50/50 distribution.