The web site of lora

#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <assert.h>

#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
#define ABS(a) (((a) < 0) ? (-(a)) : (a))

#define H (4801)
#define W (8401)

int16_t hmap[H][W]; // height map
int vismap[H][W]; // visibility map
int vismap_corr[H][W]; // visibility map corrected for leaflet.js

// upper left corner of the height map in the lon lat space
const double u0lon = 16;
const double v0lat = 49;

// position of the viewer in the lat lon space
const double lon0 = 19.1827218;
const double lat0 = 47.526613;

double lon2u(double lon) {
	return (lon - u0lon) * 1200;
}

double u2lon(double u) {
	return (u + 0.5) / 1200.0 + u0lon;
}

__attribute((constructor)) void test_lon2u() {
	assert((int)lon2u(u0lon) == 0);
	assert((int)(lon2u(u0lon + 2) + 0.5) == 2400);
	assert(ABS(u2lon(lon2u(u0lon)) - u0lon) < 0.001);
}

double lat2v(double lat) {
	return (v0lat - lat) * 1200;
}

double v2lat(double v) {
	return v0lat - (v + 0.5) / 1200.0;
}

__attribute((constructor)) void test_lat2v() {
	assert((int)lat2v(v0lat) == 0);
	assert((int)(lat2v(v0lat - 1) + 0.5) == 1200);
	assert(ABS(v2lat(lat2v(v0lat)) - v0lat) < 0.001);
}

// position of the viewer in the height map
int u0, v0;

// https://ea4eoz.blogspot.com/2015/11/simple-wgs-84-ecef-conversion-functions.html

struct xyz_t {
	double x;
	double y;
	double z;
};

void lla2xyz(struct xyz_t *xyz, double lon, double lat, double alt) {
	lon *= (M_PI / 180.0);
	lat *= (M_PI / 180.0);
	const double a = 6378137.0; // radius
	const double e = 8.1819190842622e-2; // eccentricity
	const double esq = e * e;
	const double sinlat = sin(lat);
	const double coslat = cos(lat);
	const double N = a / sqrt(1.0 - esq * (sinlat * sinlat));
	xyz->x = (N + alt) * coslat * cos(lon);
	xyz->y = (N + alt) * coslat * sin(lon);
	xyz->z = ((1.0 - esq) * N + alt) * sinlat;
}

struct lla_t {
	double lon;
	double lat;
	double alt;
};

void xyz2lla(struct lla_t *lla, double x, double y, double z) {
	const double a = 6378137.0; // radius
	const double e = 8.1819190842622e-2; // eccentricity
	const double asq = a * a;
	const double esq = e * e;
	const double b = sqrt(asq * (1.0 - esq));
	const double bsq = b * b;
	const double epsq = (asq - bsq) / bsq;
	const double p = sqrt(x * x + y * y);
	const double th = atan2(a * z, b * p);
	const double sinth = sin(th);
	const double costh = cos(th);
	lla->lon = 180.0 / M_PI * atan2(y, x);
	lla->lat = 180.0 / M_PI * atan2((z + epsq * b * (sinth * sinth * sinth)), (p - esq * a * (costh * costh * costh)));
	struct xyz_t xyz;
	lla2xyz(&xyz, lla->lon, lla->lat, 0.0);
	const double gm = sqrt(xyz.x * xyz.x + xyz.y * xyz.y + xyz.z *xyz.z);
	const double am = sqrt(x * x + y * y + z * z);
	lla->alt = am - gm;
}

__attribute((constructor)) void test_lla2xyz_xyz2lla() {
	struct xyz_t xyz;
	struct lla_t lla;

	lla2xyz(&xyz, 0, 0, 0);
	xyz2lla(&lla, xyz.x, xyz.y, xyz.z);
	assert(fabs(lla.lon - 0.0) < 0.0001);
	assert(fabs(lla.lat - 0.0) < 0.0001);
	assert(fabs(lla.alt - 0.0) < 0.0001);

	lla2xyz(&xyz, lon0, lat0, 0);
	xyz2lla(&lla, xyz.x, xyz.y, xyz.z);
	assert(fabs(lla.lon - lon0) < 0.0001);
	assert(fabs(lla.lat - lat0) < 0.0001);
	assert(fabs(lla.alt - 0.0) < 0.0001);
}

struct xyz_t xyz0;

int compute_visibility(int u1, int v1) {
	const double lon1 = u2lon(u1);
	const double lat1 = v2lat(v1);
	struct xyz_t xyz1;
	lla2xyz(&xyz1, lon1, lat1, hmap[v1][u1] + 2);
//	fprintf(stderr, "x1=%f y1=%f z1=%f\n", xyz1.x, xyz1.y, xyz1.z);

	const int du = ABS(u1 - u0);
	const int dv = ABS(v1 - v0);
	const int dmax = MAX(du, dv);

	if (!dmax) {
		return 1;
	}

	for (int ti = 1; ti <= dmax; ++ti) {
		const double t = (double)ti / dmax;
		const double x = xyz0.x + (xyz1.x - xyz0.x) * t;
		const double y = xyz0.y + (xyz1.y - xyz0.y) * t;
		const double z = xyz0.z + (xyz1.z - xyz0.z) * t;
		struct lla_t lla;
		xyz2lla(&lla, x, y, z);
		const int u = lon2u(lla.lon);
		const int v = lat2v(lla.lat);
		assert(u >= 0 && u < W);
		assert(v >= 0 && v < H);
		const int h = hmap[v][u];
		if (lla.alt < h) {
			return 0;
		}
	}
	return 1;
}

int main() {
	fprintf(stderr, "[OK]  All tests are passed.\n");

	// position of the viewer in the height map
	u0 = lon2u(lon0);
	v0 = lat2v(lat0);

	// load height map
	int prev = 0;
	for (int v = 0; v < H; ++v) {
		for (int u = 0; u < W; ++u) {
			int8_t msb = getchar();
			uint8_t lsb = getchar();
			int16_t h = msb * 256 + lsb;
			hmap[v][u] = (h == -32768) ? prev : h;
			prev = hmap[v][u];
		}
	}

	// compute view point position
	fprintf(stderr, "u0=%d v0=%d h0=%d\n", u0, v0, hmap[v0][u0] + 2);
	lla2xyz(&xyz0, lon0, lat0, hmap[v0][u0] + 2);
	fprintf(stderr, "x0=%f y0=%f z0=%f\n", xyz0.x, xyz0.y, xyz0.z);

	fprintf(stderr, "Computing visibility map...\n");
	for (int v = 0; v < H; ++v) {
		fprintf(stderr, "\r%d%%", v * 100 / H);
		for (int u = 0; u < W; ++u) {
			vismap[v][u] = compute_visibility(u, v);
		}
	}
	fprintf(stderr, "\r100%%\n");
/*
	int vismin = 5, vismax = 0;
	for (int v = 0; v < 1201; ++v) {
		for (int u = 0; u < 2401; ++u) {
			vismin = MIN(vismin, vismap[v][u]);
			vismax = MAX(vismax, vismap[v][u]);
		}
	}
	fprintf(stderr, "vismin=%d vismax=%d\n", vismin, vismax);
*/
	// distort it for leaflet.js
	for (int v = 0; v < H; ++v) {
		double a = (H/2.0 - v) / (H/2.0);
		double dv = (1.0 - a * a) * 0.0375 * 1200;
		int v2 = v - dv;
		for (int u = 0; u < W; ++u) {
			vismap_corr[v][u] = vismap[v2][u];
		}
	}

	fprintf(stderr, "Outputing visibility map...\n");
	printf("P3\n");
	printf("%d %d\n", W, H);
	printf("%d\n", 255);
	for (int v = 0; v < H; ++v) {
		fprintf(stderr, "\r%d%%", v * 100 / H);
		for (int u = 0; u < W; ++u) {
			int r = MAX(vismap_corr[v][u], 0) * 255;
			int g = 0;//MIN(vismap[v][u], 0) * 255 / vismin;
			int b = 0;
			printf("%d %d %d  ", r, g, b);
		}
		puts("");
	}
	fprintf(stderr, "\r100%%\n");
}