Over the years we commonly get asked about how to best support planetary viewsheds (in a world where Earth’s curvature and atmospheric refraction can be often hardwired into the various applications). Some of this is based on an old forum post which has long since disappeared (no nice curvature figures this time around either) and we also get to highlight new implementations just presented at the recent 2021 PDW/PSIDA workshop. Here is the email for calculating a viewshed on Europa.
Yes - you do need to be warry for any viewshed implementation. Testing against known horizon limits is ideal before trusting any application.
Speaking to hardwired curvature, fortunately, most applications do allow you to turn off both curvature and/or atmosphere refraction. But of course, then you don’t get the benefit of curvature in your calculation.
But, first, I don’t recommend coding your own implementation (since it has been implemented).
see Implementation of a Web-Based Viewshed Tool (abstract and recorded talk).
Program from 5th PDW/PSIDA workshop
–Now I’m not sure how available MMGIS is to outside mission team use (or easy to setup for other bodies). But the code is supposed to be open source, although you might need to ask for a more recent version than what I currently see on Github. Also, as it is written in JavaScript and tightly integrated into MMGIS as an interactive tool, it might be harder to extract the code as a stand-alone routine.
Anyway, the authors state they based their implementation on GDAL’s viewshed method, so I would try that first. And while I haven’t proven it to myself, I bet the GDAL implementation is planetary savvy (and of course the code is available to look at). Fortunately, GDAL code has usually been body-size agnostic.
Also, in ArcMap Pro, it looks like newer geodesic viewshed should work for planetary (NOT yet tested).
So back to testing GDAL, if you trust this site, then on a flat expanse on Europa a 1.7m tall person could see ~2.3 km to the horizon. If you get much larger viewshed distances, then you can probably assume it is using Earth (or no curvature at all).
BTW, to install GDAL I recommend Anaconda (see).
So before using GDAL’s method (or really any method – unless truly geodesic), it is best to first optimize the data’s map projection see (1). And if you are forced to use a planar implementation, meaning you needed to turn off curvature see (2).
(1) To stage a DEM for use in GDAL or testing in QGIS, ArcGIS Pro/ArcMap, Global Mapper, test running a viewshed (with any curvature off) using a local region that is in a map projection that minimizes distortions (e.g. orthographic centered on the viewshed center). In this mode, some apps let you limit the max extent to help with no curvature. Setting that should get you close for initial planetary testing).
(2) Now, if you really need curvature, in an application that doesn’t support it, first prepare the same map projection for the DEM. From planned center of the viewshed, calculate a curvature surface using the body’s radius. Subtract that new raster from the DEM and you have a DEM with built-in curvature (which will accurately represent a viewshed using non-curvature method).
Let’s see if I can describe this brute-force curvature method a little more.
(2.1) find your center location on the map (where you plan to initiate your viewshed).
(2.2) map project your DEM (bilinear) to be centered at that point using a suitable map projection to minimize distortions from that point. The could be orthographic, stereographic, Lambert equal area… Note all those will increase in distortion as you move from that center point (but probably not enough to impact your viewshed).
(2.3) Using that “center” point location, create a Euclidean distance raster (Arc toolbox tool for that) at the same cellsize and extent as your DEM (or perhaps gdal_proximity.py — GDAL documentation). The result is a raster with each pixel value growing in radial distance (in the Cartesian or planar system) away from that center. Let’s call that file “EucDist_raster.tif”.
(2.4) then apply a curvature equation in Arc (in the raster calculator or other, e.g. gdal_calc.py) on the newly created Euclidean raster layer – (where R = 1560800.0m is the 2015 IAU definition for Europa’s mean radius):
out_curve.tif = (“EucDist_raster.tif” * “EucDist_raster.tif”) / (R * 2)
(2.5) subtract that from your DEM (Arc’s raster calculator or gdal_calc.py)
dem_curve.tif = “dem.tif” - “out_curve.tif”
(2.6) Do a viewshed on the new DEM (now with curvature emplaced) using the application’s no curvature or “planar” setting (and also turn off atmospheric refraction). Compare to built-in viewshed and horizon calculation.
P.S. Horizon equation from
But just how far away is the horizon anyway? Actually, it’s pretty easy to calculate.
Over the hills and far away…
So assuming your eyes are about 5’7 [1.7m] above the ground (average for a human being), and taking Earth’s radius to be 6378.1 km, that means the horizon you see is 4.66 km away. Perhaps a little closer than you might expect.Mars, on the other hand, is indeed a much smaller planet. The radius of Mars is only 3397 km, which means that from the same height, the martian horizon would appear 3.40 km away.
The Moon is even smaller with a radius of only 1737.4 km, meaning on the moon, the horizon is only 2.43 km away.
