Idraluna Archives

The Great Antarctic... Pointcrawl?

Not sure if this is something I'm going to stick with, or just an interesting diversion (a blind alley, if you will, or a dead end, or perhaps a detour?). Solving the associated coding challenges became a single-minded fixation for a week or two. That said, the idea was born from the following considerations:

  1. Antibor has ring-shaped biomes with a vast Mongolia-esque expanse of arid steppe in the middle and temperate forests near the coasts. When mapped with hexes of uniform size, I end up with thousands of empty hexes in remote hinterlands. These don't seem like they'd be fun to traverse.
  2. Since I have GIS software and the skill to use it, why not take full advantage of the tools at my disposal?

Initially, I was inspired by Alex Schroeder's post about treating hexes as units of time to try spacing out points according to travel speed (derived from biome and average terrain ruggedness), then creating Voronoi polygons like Azgaar's Fantasy Map Generator. This would make the map slightly more visually interesting and eliminate the need for some kind of in-game movement-speed-by-terrain table -- you'd move two polygons per day, there are just more of them between you and your destination if you're in, say, a forest.

I did implement this, but ran into an initial snag because the algorithm to delineate the polygons was using Euclidean distance -- so a single polygon could span both sides of an impassable mountain ridge or a lake. To fix this, I ditched the Voronoi command and calculated (conceptually similar) 'cost allocation zones' around each point.

But at that point I thought, 'why am I enforcing an arbitrary rule that every point must be exactly a half-day apart when I have the means to calculate (and record) a precise travel time from any point to any other point?' So then I tried mapping Antibor as a pure pointcrawl.

Methods

1. Placing points of interest

The 10-mile Antibor hexgrid has around 17,000 hexes designated as '(to be) filled' in some way. In an effort to make this less of a cop-out, I decided to place approx. 20,000 points of interest on the pointcrawl map.

This was a chance to re-think how I adapted Luke Gearing's fully random hexfill procedure, specifically in terms of what categories to place and how many of each. This time around I read (and took a lot of inspiration from) Sachagoat's Re-inventing The Wilderness series. He suggests '2% Towns, 22% Scenes, 22% Lairs, 40% Dungeons*, 16% Utilities'. For a compact, sensible (actually good) pointcrawl that's great, but for Antibor I'm going to keep similar ratios for most features but shamelessly pad it with filler landmarks. So for my spin on this schema, I went with:

Points of interest were placed by running a crude site suitability analysis. Along with a modicum of verisimilitude, the aim here was to inject underlying logic into the pattern of settlements and lairs for players to (hopefully) grasp and interact with (e.g. 'if we stay along the river we'll mostly encounter human settlements' or 'we should grab extra supplies before venturing into the highlands, since there will be more monster lairs there'). I was also interested in seeing if I could get natural clusters to form, possibly functioning as sub-pointcrawls akin to how I was using 'Prefectures' before.

Early version of the settlement suitability scores, lighter color = more suitable.

2. Cost allocation zones

Once the points were placed, I calculated a cost distance raster where each pixel records the cost (time or energy) of travel to the nearest POI. This was done by assigning each biome an average miles-per-day value ranging from 30 in steppe to 14 in jungles.

To carve the continuous raster into discrete zones, I used marker-controlled watershed segmentation. This is usually used for watershed basins (like the name suggests) or to pick out individual trees from a 3d tree canopy scan, but it works to segment out anything composed of discrete 'bumps'.

At this point, the result is equivalent to Voronoi polygons, but for a non-Euclidean 'travel time space'. This could be a valid end product, but I also wanted to draw in paths between the POIs.

3. Pseudo-Delaunay paths

The inverse of a Voronoi polygon is a Delaunay Triangulation. These have a bunch of cool properties, but for pointcrawl purposes they're a good method for linking up random points in a way that looks natural and doesn't result in crossing or overlapping paths: if the Voronoi polygons touch, the points get a path.

For this, I returned to the leastcostpath package for R, since it was designed to study ancient roads. The self-imposed rules for plotting a path looked like this:

I also added some math to randomly determine if a connection is a road, a trail, or an unmarked wilderness 'path' -- if both nodes are villages, the odds are highest, ruins and landmarks in the middle, and lairs have the lowest likelihood. Essentially, roll once: if above the threshold, it's a road, other wise roll again. If above the threshold this time, it's a trail, otherwise it's nothing. The result seems ok, but I may manually adjust some to produce more logical networks.

Most POIs are separated by paths of less than 20 miles:

What I most wanted was paths that snake through valleys and mountain passes (non-colored areas are those deemed too steep to traverse normally):

leastcostpath attaches a somewhat inscrutable Cost variable to each path reflecting the value the algorithm minimized when generating it. I ran a linear model and found that the line of best fit relating Euclidian path length (in mi) to Cost had a slope of ~1300. So by dividing Cost by that amount I was able to get something like an 'adjusted miles' value that accounts for a path's difficulty -- I envision using this as its 'length' when running a game.

4. To Markdown (again)

As with the previous GAH post. This time, I simplified the process by using the Json/CSV Importer plugin for Obsidian along with a hex template file -- this allowed me to just export the POIs as a table rather than recycling my code to write 20,000 .md files. To list the path connections and respective travel times, I looped through each point and grabbed the necessary info from the paths file.

Results

When all the code finished crunching I loaded the results into QGIS and started exploring. Had some fun picking out icons for the different POI types, etc.

I was able to get the distance and travel time to display along each route by using the 'Curved' label placement option and the following expression: concat(ceil("Mi_cost"/3), 'h. (',"Length_mi" ,'mi.)').

To revel in non-hexagonal cartographic freedom, I threw together a nice-ish map layout with the full results. (Click here for full size).

Pros & Cons

I like the initial results a lot, but I'm not fully committed to abandoning the hexmap. Even though the actual number of keyed sites is about the same, on some level a pointcrawl just feels less impressive than a continent-scale hexcrawl.

Pros

Cons

Swapping back and forth between pointcrawl and hexcrawl is viable, either using hexes as a 'an interesting type of ruler that's overlaid directly atop the map' or by keying up the pointcrawl sites and then spatially joining them back to the hexmap later. This latter option might be a viable way to utilize a smaller hex grid, like 6 or even 3 miles, where the paths (possibly snapped to hex centers) help orient players toward the keyed locations.

Appendix: Code

library(terra)
library(tidyverse)
library(tidyterra)

mapdir <- '../Maps'
tabledir <- '../Tables'
chapterdir <- '../Chapters'
hexdir <- '../Antibor'
set.seed(42069)

Biome_map <- data.frame(
  'Code' = c(0, 1,2,4,5,6,7,8,10,11,12,13,14,15, 16),
  'Biome_name' = c('Water', 'Jungle', 'Forest (warm)', 'Forest (temperate)', 'Conifers (temperate)', 'Taiga', 'Savannah', 'Steppe', 'Tundra', 'Tundra', 'Chaparral', 'Desert', 'Urban', 'Citadel', 'Xenoformed'),
  'Hab_score' = c(NA, 36, 94, 63, 63, 25, 94, 15, 5, 5, 125, 30, NA, NA, 5),  # scale of 1-10
  'Friction_mult_trail' = c(500, 3, 2, 2, 2, 2.5, 1, 1, 1.5, 1.5, 1, 1.1, 1, 1, 10),
  'Friction' = c(500, 90, 35, 35, 35, 35, 10, 20, 50, 50, 10, 20, 0, 0, 100),
  'MPD' = c(NA, 14, 18, 20, 20, 20, 26, 30, 30, 30, 24, 30, NA, NA, 12)
)

# Load some useful layers
biome_base.r <- rast(file.path(mapdir, 'biomes2.tif')) %>% trim()
elevation_full.r <- rast(file.path(mapdir, 'Antarctica_rebound_sealevelrise.tif'))
elevation.r <- elevation_full.r %>% crop(biome_base.r)
hexes.shp <- vect(file.path(mapdir, 'Full_hexes.gpkg'))
citadels.r <- hexes.shp %>% filter(Biome=='Citadel') %>%
  rasterize(biome_base.r)
xeno.shp <- hexes.shp %>% filter(Biome=='Xenoformed') 
xeno.r <- xeno.shp %>%
  rasterize(biome_base.r) %>%
  focal(w=5, fun=function(m){sample(m, 1)})
big_rivers.shp <- vect(file.path(mapdir, 'rivers_big_smooth.gpkg'))
big_rivers.r <- big_rivers.shp %>%
  rasterize(biome_base.r)
small_rivers.r <- vect(file.path(mapdir, 'rivers_smooth.gpkg')) %>%
  rasterize(biome_base.r)
water_bodies.r <- ifel(is.na(biome_base.r), 1, NA) %>% patches() %>%
  zonal(cellSize(., unit='km'),.,sum, as.raster=T)
ocean.r <- ifel(water_bodies.r>10000000, 1, NA)
lakes.r <- ifel(water_bodies.r<10000000 & water_bodies.r>25, 1, NA)
slope.r <- elevation.r %>% terrain()
slope_barriers.r <- ifel(slope.r<7.5, 0, 1)
terrain_var.r <- elevation.r %>%
  focal(w=11, fun = 'sd') %>%  # a diameter of 11 km is around 5-6 mi, or 1/3rd day walking
  crop(biome_base.r, mask=T)
terrain_var.r[slope_barriers.r==1] <- 0
mpd_penalty.r <- (1000-terrain_var.r)/1000
big_roads.shp <- vect(file.path(mapdir,'Prefec_roads.gpkg')) 
big_roads.r <- big_roads.shp %>%
  rasterize(biome_base.r)
biome.r <- biome_base.r
biome.r[is.na(biome.r)] <- 0
biome.r[big_rivers.r==1] <- 0
biome.r[xeno.r==1] <- 16
biome.r[citadels.r==1] <- 15
biomes.shp <- as.polygons(biome.r) %>% left_join(Biome_map, by=c('class'='Code'))

MPD.r <- biome.r %>%
  classify(rcl = data.frame('From'=Biome_map$Code, 'To'=Biome_map$MPD))
MPD.r <- MPD.r * mpd_penalty.r
MPD.r[slope_barriers.r==1] <- NA
MPD.r[big_roads.r==1] <- 32
MPD.r[big_rivers.r==1&!is.na(biome_base.r)] <- 1
plot(MPD.r)
mPD.r <- MPD.r * 1608  # meters per day

###############################################################
### Compile Points of Interest ################################
###############################################################
# Named hex locations -- snap to cell centers 
named_hex_pts.shp <- hexes.shp %>% filter(!is.na(Name)) %>% centroids() %>% select(Name, Region, Biome, Subregion, Prefecture) %>%
  mutate(Type = 'Imported')
nhp_cells.shp <- named_hex_pts.shp %>% rasterize(biome.r) %>% as.points()
named_hex_pts.shp <- named_hex_pts.shp %>% snap(nhp_cells.shp, 2000)

# Villages
hab_score.r <- biome.r %>%
  classify(rcl = data.frame('From'=Biome_map$Code, 'To'=Biome_map$Hab_score)) 

dist_freshwater.r <- costDist(ifel(big_rivers.r==1|lakes.r==1, -99, 1/mPD.r), target=-99)
dist_ocean.r <- costDist(ifel(is.na(mPD.r), -99, 1/mPD.r), target=-99)
dist_road.r <- costDist(ifel(big_roads.r==1, -99, 1/mPD.r), target=-99)

town_raster <- named_hex_pts.shp %>% filter(Biome=='Urban') %>% rasterize(biome_base.r)
dist_town.r <- costDist(ifel(town_raster==1, -99, 1/mPD.r), target=-99)

par(mfrow=c(1,3));plot(dist_freshwater.r);plot(dist_road.r);plot(dist_town.r)
dev.off()

dist_freshwater_quantile.r <- dist_freshwater.r %>%
  classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))
dist_road_quantile.r <- dist_road.r %>%
  classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))
dist_town_quantile.r <- dist_town.r %>%
  classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))
ruggedness_quantile.r <- terrain_var.r %>%
  classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))

dist_road_quantile.r[is.na(dist_road_quantile.r)&!is.na(hab_score.r)] <- 10
dist_town_quantile.r[is.na(dist_town_quantile.r)&!is.na(hab_score.r)] <- 10
dist_freshwater_quantile.r[is.na(dist_freshwater_quantile.r)&!is.na(hab_score.r)] <- 10

village_suitability.r <- hab_score.r - 3*dist_freshwater_quantile.r - 2*dist_town_quantile.r - 3*dist_road_quantile.r - ruggedness_quantile.r + 10*dist_ocean.r**-0.75
plot(village_suitability.r, col=viridis::magma(100))

village_suitability_quantile.r <- village_suitability.r %>% classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))
plot(village_suitability_quantile.r, col=viridis::magma(10))

village_sites.shp <- village_suitability_quantile.r %>% as.numeric() %>%
  as.data.frame(xy=T) %>%
  mutate(prob = class**1.25) %>%  # best areas around 11x as likely as worst
  slice_sample(n = 3531, weight_by = prob) %>%
  mutate(Type = 'Settlement') %>%
  vect(geom=c('x', 'y'), crs=crs(biome_base.r))

plot(village_sites.shp, add=T, col='green')

# Landmarks
landmark_sites.shp <- village_suitability_quantile.r %>% as.numeric() %>%
  as.data.frame(xy=T) %>%
  mutate(prob = class+1) %>%  # best areas around 2.2x as likely as worst
  slice_sample(n = 8007, weight_by = prob) %>%
  mutate(Type = 'Landmark') %>%
  vect(geom=c('x', 'y'), crs=crs(biome_base.r))

plot(landmark_sites.shp, add=T, col='blue')

# Lairs
lair_suitability.r <- 0.33*hab_score.r + dist_town_quantile.r + dist_road_quantile.r - dist_freshwater_quantile.r + 3*ruggedness_quantile.r + ifel(!is.na(xeno.r), 30, 0) + 5*dist_ocean.r**-0.75

lair_suitability_quantile.r <- lair_suitability.r %>% classify(t(global(., quantile, probs=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1), na.rm=TRUE)))
plot(lair_suitability_quantile.r)

lair_sites.shp <- lair_suitability_quantile.r %>% as.numeric() %>%
  as.data.frame(xy=T) %>%
  mutate(prob = class+1) %>%  
  slice_sample(n = 4042, weight_by = prob) %>%
  mutate(Type = 'Lair') %>%
  vect(geom=c('x', 'y'), crs=crs(biome_base.r))

plot(lair_sites.shp, add=T, col='red')

# Ruins
ruin_sites.shp <- lair_suitability_quantile.r %>% as.numeric() %>%
  as.data.frame(xy=T) %>%
  mutate(prob = class+2) %>%  
  slice_sample(n = 4444, weight_by = prob) %>%
  mutate(Type = 'Ruin') %>%
  vect(geom=c('x', 'y'), crs=crs(biome_base.r))

plot(ruin_sites.shp, add=T, col='gray')

all_sites.shp <- rbind(named_hex_pts.shp, village_sites.shp, landmark_sites.shp, lair_sites.shp, ruin_sites.shp) %>% 
  distinct(geometry, .keep_all = TRUE)
writeVector(all_sites.shp, file.path(mapdir, 'Pointcrawl_sites.gpkg'), overwrite=T)

###############################################################
### Generate Cost Allocation Map ##############################
###############################################################
POI.r <- all_sites.shp %>% rasterize(mPD.r)

cd_from_POI.r <- costDist(ifel(!is.na(POI.r), -99, 1/mPD.r), target=-99)
plot(cd_from_POI.r)

point_ids.shp <- sf::st_as_sf(all_sites.shp %>% mutate(treeID = row_number()))

seg.r <- ForestTools::mcws(
  treetops = point_ids.shp,
  CHM = cd_from_POI.r*-1+10
)

plot(seg.r)
writeRaster(seg.r, file.path(mapdir, 'seg_pointcrawl.tif'), overwrite=T)

###############################################################
### Generate linkages #########################################
###############################################################
gc()
seg_poly.shp <- seg.r %>% as.polygons()

seg_adjacencies.df <- adjacent(seg_poly.shp, pairs=T, symmetrical=T)
library(leastcostpath)
friction_stack.r <- c(ifel(is.na(mPD.r), 0.1, mPD.r), add_dem_error(crop(elevation_full.r, mPD.r, mask=T), rmse = 50, type='u'))

paths.shp <- vect()

for (i in 1:nrow(seg_adjacencies.df)){
  if(i%%5000==0){gc()
    writeVector(paths.shp, file.path(mapdir, 'pointcrawl_paths.gpkg'), overwrite=T)}
  print(paste(i, 'of', nrow(seg_adjacencies.df)))
  
  cnx_pair <- seg_adjacencies.df[i,]
  zones_to_link.shp <- seg_poly.shp[cnx_pair,]
  
  travel_zones.r <- seg.r %>% crop(ext(zones_to_link.shp)) %>%
    ifel(.%in%zones_to_link.shp$lyr.1, ., NA)
  #plot(travel_zones.r)
  
  points_to_link.shp <- point_ids.shp %>% filter(treeID %in% zones_to_link.shp$lyr.1) %>% vect() %>%
    mutate(Road_chance = ifelse(Type%in%c('Settlement', 'Imported'), .4, 
                                ifelse(Type%in%c('Ruin', 'Landmark'), .25, .05)))
  
  elev_subset.r <- friction_stack.r[[2]] %>% crop(travel_zones.r, mask=T)
  cs.r <- create_slope_cs(elev_subset.r)#; plot(cs.r)
  cs.r <- update_values(x=cs.r, sf=sf::st_as_sf(xeno.shp), FUN = function(j) { j * 0.33})  # xenoformed 3x as hard to traverse
  cs.r <- update_values(x=cs.r, sf=sf::st_as_sf(big_rivers.shp), FUN = function(j) { j * .01})  # rivers 100x as hard
  cs.r <- update_values(x=cs.r, sf=sf::st_as_sf(big_roads.shp), FUN = function(j) { j * 2})
  
  if(inherits(tryCatch(check_locations(x=cs.r, locations=points_to_link.shp), error=function(e){e}), 'error')){
    print(points_to_link.shp$TreeID)
    path.shp <- as.lines(points_to_link.shp) 
    path.shp$origin <- zones_to_link.shp$lyr.1[1]
    path.shp$destination <- zones_to_link.shp$lyr.1[2]
    path.shp$Road_chanceA <- 0
    path.shp$Road_chanceB <- 0
    path.shp$Type <- 'Other'
    paths.shp <- rbind(paths.shp, path.shp)
    next
  }
  
  path.shp <- create_lcp(cs.r, origin = points_to_link.shp[1,], destination = points_to_link.shp[2,], cost_distance = T)
  #plot(travel_zones.r);plot(path, add=T)
  
  # figure out if the connection should be a road, trail, or other
  path.shp <- path.shp %>% mutate(
    origin = zones_to_link.shp$lyr.1[1],
    destination = zones_to_link.shp$lyr.1[2],
    Road_chanceA = points_to_link.shp$Road_chance[1],
    Road_chanceB = points_to_link.shp$Road_chance[2]
  )
  
  paths.shp <- rbind(paths.shp, path.shp)
}

paths.shp <- paths.shp %>%
  mutate(Road_chance = Road_chanceA+Road_chanceB) %>%
  rowwise() %>%
  mutate(Type = ifelse(runif(1)<Road_chance**2, 'Road', 
                       ifelse(runif(1)<Road_chance**2, 'Trail', 'None')))

paths.shp <- paths.shp %>% ungroup() %>%
  mutate(Length_mi = as.integer(perim(.)/1609)) %>% mutate(Mi_cost = ceiling(cost/1300))

writeVector(paths.shp, file.path(mapdir, 'pointcrawl_paths.gpkg'), overwrite=T)

###############################################################
### Generate markdown #########################################
###############################################################

POIs.shp <- vect(file.path(mapdir, 'Pointcrawl_sites.gpkg')) %>%
  mutate(Name = ifelse(is.na(Name), paste('POI', ID), Name))
POI_paths.shp <- vect(file.path(mapdir, 'pointcrawl_paths.gpkg'))
POI_paths.df <- as.data.frame(POI_paths.shp) %>%
  mutate(Road_type = tolower(ifelse(Type == 'None', 'through wilderness', Type)))

hexes.shp <- vect(file.path(mapdir, 'Full_hexes.gpkg'))
biome.r <- rast(file.path(mapdir, 'biomes_complete.tif'))

to_md.df <- as.data.frame(POIs.shp, geom='XY') %>%
  mutate(coord_id = paste0(x, '.', y))

## calculate neighbors from roads
for(POI.id in to_md.df$ID){
  print(POI.id)
  paths <- POI_paths.df %>%
    filter(origin==POI.id|destination==POI.id) %>%
    mutate(other_POI = ifelse(origin==POI.id, destination, origin)) %>%
    left_join(to_md.df, by=c('other_POI'='ID')) %>%
    mutate(phrase = paste('\n-', Mi_cost, 'adj. mi.', Road_type, 'to', paste0('[[', Name, ']]')))
  
  to_md.df[to_md.df$ID==POI.id, 'Neighbors'] <- paste(paths$phrase, collapse='')
}

write.csv(to_md.df, file.path(mapdir, 'pointcrawl_POI_export.csv'))

#DIY #GIS #antibor #lore24