Skip to content

Instantly share code, notes, and snippets.

@Dan-Q Dan-Q/geohash-pcwater.js
Last active Aug 22, 2018

Embed
What would you like to do?
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/
/*
* 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}%`);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.