Perform pixel-analysis on OpenStreetMap tiles to estimate water coverage of a graticule. More info: https://danq.me/2018/08/21/graticule-water-coverage-calculator/

geohash-pcwater.js

/*
 * More details can be found at:
 * https://danq.me/2018/08/21/graticule-water-coverage-calculator/
 * 
 * Given a graticule (e.g. 51 -1), returns the percentage water cover
 * of that graticule based on pixel colour sampling of OpenStreetMap
 * tile data. Change the zoomLevel to sample with more (higher) or less
 * (lower) granularity: this also affects the run duration. Higher
 * granularity improves accuracy both by working with a greater amount
 * of data AND by minimising the impact that artefacts (e.g. text,
 * borders, and ferry lines, which are usually detected as land) have
 * on the output.
 *
 * Expects a Unix-like system. Requires grep, wc, wget, and "convert"
 * (from the ImageMagick suite). And a Javascript interpreter (e.g.
 * node), of course. On a Debian/Ubuntu-like distro, all non-node
 * dependencies can probably be met with:
 *   sudo apt install -y wget imagemagick
 *
 * (c) Dan Q [danq.me] 2018; no warranty expressed or implied; distribute
 *     freely under the MIT License (https://opensource.org/licenses/MIT)
 *
 * Sample outputs:
 *   $ node geohash-pcwater.js 51 -1 # Swindon, Oxford (inland)
 *   ...
 *   Water ratio: 0.68%
 *
 *   $ node geohash-pcwater.js 49 -2 # Channel Islands (islandy!)
 *   ...
 *   Water ratio: 93.13%
 */
const { execSync } = require('child_process');

const lngToTile = (lng, zoom)=>(Math.floor((Number(lng)+180)/360*Math.pow(2,zoom)));
const latToTile = (lat, zoom)=>(Math.floor((1-Math.log(Math.tan(lat*Math.PI/180) + 1/Math.cos(lat*Math.PI/180))/Math.PI)/2 *Math.pow(2,zoom)));
const urlRange = (lng1, lng2, lat1, lat2, zoom)=>{
  let urls = [];
  const x1 = lngToTile(lng1, zoom), x2 = lngToTile(lng2, zoom), y1 = latToTile(lat1, zoom), y2 = latToTile(lat2, zoom);
  for(let x = Math.min(x1, x2); x <= Math.max(x1, x2); x++){
    for(let y = Math.min(y1, y2); y <= Math.max(y1, y2); y++){
      const server = String.fromCharCode(Math.floor(Math.random() * 3) + 97);
      const url = `https://${server}.tile.openstreetmap.org/${zoom}/${x}/${y}.png`;
      urls.push(url);
    }
  }
  return urls;
}

if(process.argv.length < 4){
  console.log('Syntax: node geohash-pcwater.js 51 -1 (where 51 -1 is your graticule)');
  process.exit();
}

const graticule = [process.argv[2], process.argv[3]];
const zoomLevel = 10; // OpenStreetMap zoom level; impacts granularity of sampling: each time you add one you'll approximately quadruple the number of tiles to download for a given graticule - 10 seems nice

graticuleTop = (graticule[0][0] == '-' ? parseInt(graticule[0]) : parseInt(graticule[0]) + 1);
graticuleBottom = (graticule[0][0] == '-' ? parseInt(graticule[0]) - 1 : parseInt(graticule[0]));
graticuleLeft = (graticule[1][0] == '-' ? parseInt(graticule[1]) - 1 : parseInt(graticule[1]));
graticuleRight = (graticule[1][0] == '-' ? parseInt(graticule[1]) : parseInt(graticule[1]) + 1);

const images = urlRange(graticuleLeft, graticuleRight, graticuleTop, graticuleBottom, zoomLevel);
console.log(`${images.length} images must be processed...`)

let pxTotal = 0, pxWater = 0;
for(let url of images){ // for each tile...
  console.log(`Fetching ${url}:`);
  execSync(`wget ${url} -qO tmp.png`); // use wget to download the tile
  console.log(' > extracting data');
  execSync('convert tmp.png tmp.txt'); // use imagemagick to extract the data as text
  console.log(' > analysing')
  pxTotal += (parseInt(execSync('cat tmp.txt | wc -l').asciiSlice().trim()) - 1); // wc/grep the text
  pxWater += (parseInt(execSync('grep -Ei "#(abd3df|aad3df)" tmp.txt | wc -l').asciiSlice().trim())); // abd3df and aad3df are the hex codes for the two colours I've seen of water
}
pcWater = Math.round((pxWater / pxTotal) * 10000) / 100;
console.log(`Water ratio: ${pcWater}%`);