Skip to content

Commit

Permalink
#398 Final touchups, docs, includeSunEarth, bug fixes, utc instead of…
Browse files Browse the repository at this point in the history
… tdb
  • Loading branch information
tariqksoliman committed Sep 6, 2023
1 parent ec2d27e commit e1cb586
Show file tree
Hide file tree
Showing 13 changed files with 668 additions and 294 deletions.
146 changes: 146 additions & 0 deletions docs/pages/Tools/Shade/Shade.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
---
layout: page
title: Shade
permalink: /tools/shade
parent: Tools
---

# Shade

_Shades the ground when line-of-sights to an orbiting target are occluded._

## SPICE

Site administrators are responsible for keeping SPICE kernels up-to-date in `/private/api/spice/kernels` and CHRONOS setup files relevant in `/private/api/spice/chronosSetups`.

There are two SPICE python scripts that require these backend kernel setups:

- `/private/api/spice/chronos.py`
- Converts between time systems.
- Looks for `/private/api/spice/chronosSetups/chronos-{target}.setup` where `{target}` here is filled in as a lowercased ShadeTool variables `"observers"`s `"value"`.
- `/private/api/spice/ll2aer11.py`
- Turns a lnglat and target into a directional azimuth, elevation, range, and lntlat
- Reads in all kernels `/private/api/spice/kernels`.
- `/private/api/spice/getKernelUtils` has some wget scripts as examples for downloadign new kernels (however these resources will quickly become outdated)

## Tool Configuration

### Example

```javascript
{
"dem": "Data/missionDEM.tif",
"data": [
{
"name": "MSL_DEM",
"demtileurl": "pathToDEMTiles/MSL_Gale_DEM_Mosaic_1m_v3/{z}/{x}/{y}.png",
"minZoom": 8,
"maxNativeZoom": 18
}
],
"sources": [
{
"name": "MRO",
"value": "MRO"
},
{
"name": "ODY",
"value": "-53"
},
{
"name": "TGO",
"value": "TGO"
},
{
"name": "MVN",
"value": "MAVEN"
},
{
"name": "The Sun",
"value": "SUN"
}
],
"observers": [
{
"name": "MSL",
"value": "MSL"
}
]
}
```

_**dem**_ - A path to a DEM.tif. This is used to get the current center elevation. This can/should be the same file used for the Measure Tool and the Coordinate's elevation.

_**data**_ - At minimum, the Shade tool requires at least one "data" source. A data source describes a DEM tileset (see /auxiliary/gdal2customtiles or /auxiliary/1bto4b) and allows users to select it by name to generate shade maps over.

_**source**_ - An array of objects with the properties "name" and "value". "name" is the display name for the Source Entity dropdown. "value" is the SPICE spacecraft ID that gets passed to the backend `ll2aerll.py` script. Ensure the right kernels for the configured source entities/targets exist in `/private/api/spice/kernels`.

_**observers**_ - An array of objects with the properties "name" and "value". "name" is the display name for the Source Entity dropdown. "value" is the SPICE spacecraft ID that gets passed to the backend `chronos.py` scripts. Ensure the right kernels for the configured observers exist in `/private/api/spice/kernels` and that there is a proper chronos setup file for each observer's value `private/api/spice/chronosSetups/chronos-{lowercased_observer_value}.setup`.

## Tool Use

**Note:** Terrain beyond the screen's current extent is **not** factored into the displayed visiblity map — only observer-target direction and on-screen terrain is considered. A distant off-screen mountain will **not** cast shadows.

### Interface

- _Time_
- The desired datetime to query. Formatted as `YYYY MMM DD HH:MM:SS` and for example `2023 SEP 06 19:27:05`. Updating this time and pressing 'Enter' will set it as the current time for the ShadeTool and for all of MMGIS. It is both connected to the Observer's local time as well as MMGIS' timeline (expandable via the clock icon in the bottom left of the screen).

#### Source

- _Entity_
- Indicates which spacecraft, orbiter or celestial body to "look towards" and to "shine light back" upon the visible terrain.
- _Include Sun + Earth_
- If true, the relative Sun and Earth positions will also be computed and their directional arrows will be rendered in the bottom azimuth and elevation indicators. In the azimuth and elevation indicators, the Sun is represented by a medium-length yellow arrow and the Earth is represented by a short-length blue-green arrow. These do **not** cast shadows on the visible terrain — only the source entity casts shadows.

#### Observer

- _Entity_
- Which observing spacecraft/orbiter to use. This is only used for formatting and converting the upcoming 'Time' parameter. The true observer position is always the visible map's center longitude and latitude value (represented by a green circle) and always facing north with zero tilt.
- _Time_
- Offers the ability to set the current working time using a mission/spacecraft's custom date type.
- _Height_
- Height in meters above the surface to use when calculating line-of-sight shading. For instance, a point on the surface (0m) may not be visible to a 'Source Entity', say the Mars Reconnaissance Orbiter (MRO), but 2m above that point may be. This value does not _only_ apply to the center longtitude and latitude but to all points on the visible terrain. Gradually increaing this value shows the shade map n-meters above the surface.

#### Shaded Region Options

- _Color_
- The color to shade the shadowed regions on the map.
- _Opacity_
- The opaqueness to shade the shadowed regions on the map. A value of 0 is fully transparent and a value of 1 is fully opaque.
- _Resolution_
- MMGIS downloads terrain data needed for the shading alogrithm. Increasing the resolution improves the quality of the shade map and the cost of download and render speed. Each higher option is 4x the resolution of the previous one (i.e. 'ultra' is 4x more terrain data than 'high' and 16x more data than 'medium'). To save on performance, if the resolution is 'high' or 'ultra', the Shade Tool will no longer regenerate the shaded map whenever any parameter changes and instead 'Generate/Regenerate' must manually be pressed.
- _Elevation Map_

- Specifies the terrain dataset to use.

- _Generate/Regenerate_
- Submits a request to generate a shade map with the provided parameters. Note that if the resolution is 'high' or 'ultra', the Shade Tool will not regenerate the shaded map whenever any parameter changes and instead 'Generate/Regenerate' must manually be pressed.

#### Results

- _Azimuth_: The compass-angle in (0 -> 360) degrees clockwise from north of the direction of the 'Source Entity' as seen from the map's center longitude and latitude. 0 = North, 90 = East, 180 = South, 270 = West.
- _Elevation_: The angular height (-90 -> 90) between the horizon and the 'Source Entity'. -90 = Straight Down, 0 = Level with the Horizon, 90 = Straight Overhead.
- _Range_: The straight-line distance in kilometers between the map's center longitude, latitude and terrain elevation and the 'Source Entity'.
- _Longitude_: The map's center longitude value used in the computation.
- _Latitude_: The map's center latitude value used in the computation.
- _Altitude_: The distance in kilometers above the map's center position's tangential plane and the 'Source Entity'. In other words, in a 3D cartesian coordinate-system where the Z-axis goes through both the center of the visible map and the center of the planet, this 'Altitude' is the Z distance between that center and the 'Source Entity'.

#### Indicators

- _Azimuth_: A top-down birds-eye view of the surface with north up. The long yellow-orange arrow visualizes the azimuthal direction towards the 'Source Entity'. If 'Include Sun + Earth' is on, shorter Sun and Earth arrows will also appear in the indicator with the respective yellow and green-blue colors.
- _Elevation_: A horizontal and half-submerged side view of the surface. The long yellow-orange arrow visualizes the elevational direction towards the 'Source Entity'. If 'Include Sun + Earth' is on, shorter Sun and Earth arrows will also appear in the indicator with the respective yellow and green-blue colors. Note that elevation values only goes from -90 -> 90 but that the rendered elevation arrow can be drawn between 0 -> 360. This is because, while only half a circle is needed, the elevation arrow will choose whether to draw in the left or right half circle depending on which half-circle the azimuth value is in. Azimuth values from 0 -> 180 will result in an elevation arrow drawn in the right half-circle and azimuth values from 180 -> 360 will results in an elevation arrow drawn in the left half-circle. This is to aid in visualizing the 'Source Entity's 3D direction.

### Algorithm

1. The following are taken and fed into SPICE:
- The longitude, latitude, and elevation location at the center of the map (the observer)
- The current date
- The target/source-entity (which may be an orbiter, the Sun, etc.)
2. The following are returned
- The azimuth, elevation, and range from that location to the target
- Source location is assumed to be facing north with no tilt
- The longitude, latitude on the map directly under the target and its elevation
3. All elevation values from the current screen extent and queried
4. The target's longitude, latitude, elevation are projected onto a plane tangential to the observer
5. The screen elevation values are placed in an xy grid and the, from the previous values, the target's respective x,y,elev is computed and run through a modified version of [_Generating Viewsheds without Using Sightlines_](https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf) by _Jianjun Wang, Gary J. Robinson, and Kevin White_
Binary file added private/api/spice/chronos
Binary file not shown.
12 changes: 10 additions & 2 deletions private/api/spice/chronos.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
import os
import subprocess
import shlex
import platform


try:
from urllib.parse import unquote
Expand All @@ -22,9 +24,15 @@
def chronos(target, fromFormat, fromtype, to, totype, time):
package_dir = os.path.dirname(os.path.abspath(__file__)).replace('\\','/')

cmd = os.path.join(package_dir + '/', 'chronos.exe')
plt = platform.system()

if plt == "Windows":
cmd = os.path.join(package_dir + '/', 'chronos.exe')
else:
cmd = os.path.join(package_dir + '/', 'chronos')

target = target.replace('\\','').replace('/','')
setup = os.path.join(package_dir + '/', f'chronos/chronos-{target}.setup')
setup = os.path.join(package_dir + '/', f'chronosSetups/chronos-{target}.setup')
fullCmd = shlex.split(f'{cmd} -setup {setup} -from {fromFormat} -fromtype {fromtype} -to {to} -totype {totype} -time {time} -NOLABEL')

result = subprocess.run(fullCmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True)
Expand Down
File renamed without changes.
155 changes: 30 additions & 125 deletions private/api/spice/ll2aerll.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,10 @@
import sys
import json
import os
import subprocess
import shlex
import uuid
import numpy as np

import math
import spiceypy
import pymap3d as pm
import re


# chronos.exe -setup ./chronos/chronos-msl.setup -from utc -fromtype scet -to lst -totype lst -time "2023-07-27 23:16:05.644"
Expand All @@ -26,7 +21,7 @@
except ImportError:
from urllib import unquote

def ll2aerll(lng, lat, height, target, time):
def ll2aerll(lng, lat, height, target, time, includeSunEarth):
# Load kernels
package_dir = os.path.dirname(os.path.abspath(__file__)).replace('\\','/')

Expand All @@ -37,30 +32,6 @@ def ll2aerll(lng, lat, height, target, time):
# Prints only text file present in My Folder
kernels_to_load.append(x)

# Fill out dynamic template
with open (os.path.join(package_dir, 'pinpoint/dynamic-template.setup'), 'r' ) as f:
content = f.read()
content_new = re.sub('{{LNG}}', str(lng), content, flags = re.M)
content_new = re.sub('{{LAT}}', str(lat), content_new, flags = re.M)
content_new = re.sub('{{HEIGHT_KM}}', str(height / 1000), content_new, flags = re.M)

fileTf = open(os.path.join(package_dir + '/kernels/', 'dynamic.setup'), "w")
fileTf.write(content_new)
fileTf.close()

# Create dynamic
dynamicSpkFile = f'dynamic-{str(uuid.uuid4())}.bsp'
dynamicFkFile = f'dynamic-{str(uuid.uuid4())}.tf'
cmd = os.path.join(package_dir + '/', 'pinpoint.exe')
setup = os.path.join(package_dir + '/kernels/', 'dynamic.setup')
pck = os.path.join(package_dir + '/kernels/', 'mars_iau2000_v1.tpc')
dynamicSpk = os.path.join(package_dir + '/kernels/', dynamicSpkFile)
dynamicFk = os.path.join(package_dir + '/kernels/', dynamicFkFile)
subprocess.call(shlex.split(f'{cmd} -def {setup} -pck {pck} -spk {dynamicSpk} -fk {dynamicFk}'))

kernels_to_load.append(dynamicSpkFile)
kernels_to_load.append(dynamicFkFile)

for k in kernels_to_load:
spiceypy.furnsh( os.path.join(package_dir + '/kernels/', k) )

Expand Down Expand Up @@ -90,28 +61,28 @@ def ll2aerll(lng, lat, height, target, time):

razel = output[0]


# B
"""
obs = "-654321"
ref = "IAU_MARS"
try:
state = spiceypy.spkezr( target, et, ref, abcorr, obs )
razel = spiceypy.recazl( [state[0][0], state[0][1], state[0][2]], azccw, elplsz )
except Exception as e:
# Unload kernels
for k in kernels_to_load:
spiceypy.unload( os.path.join(package_dir + '/kernels/', k) )
# Remove temporary dyanmic files
os.remove(dynamicSpk)
os.remove(dynamicFk)
raise e
"""
az_output = razel[1] * spiceypy.dpr()
el_output = razel[2] * spiceypy.dpr()
range_output = razel[0] * 1000
range_output = razel[0]

sun_az_output = None
sun_el_output = None
earth_az_output = None
earth_el_output = None
if includeSunEarth == 'true':
target = 'SUN'
sun_obspos = spiceypy.georec( lng * spiceypy.rpd(), lat * spiceypy.rpd(), height / 1000, radii[0], flattening)
sun_output = spiceypy.azlcpo( method, target, et, abcorr, azccw, elplsz, sun_obspos, obsctr, obsref)
sun_razel = sun_output[0]
sun_az_output = sun_razel[1] * spiceypy.dpr()
sun_el_output = sun_razel[2] * spiceypy.dpr()

target = 'EARTH'
earth_obspos = spiceypy.georec( lng * spiceypy.rpd(), lat * spiceypy.rpd(), height / 1000, radii[0], flattening)
earth_output = spiceypy.azlcpo( method, target, et, abcorr, azccw, elplsz, earth_obspos, obsctr, obsref)
earth_razel = earth_output[0]
earth_az_output = earth_razel[1] * spiceypy.dpr()
earth_el_output = earth_razel[2] * spiceypy.dpr()

try:
visibilities = [] #lltarg2vistimes(1,2,3,4,5)
Expand All @@ -122,9 +93,6 @@ def ll2aerll(lng, lat, height, target, time):
for k in kernels_to_load:
spiceypy.unload( os.path.join(package_dir + '/kernels/', k) )

# Remove temporary dyanmic files
os.remove(dynamicSpk)
os.remove(dynamicFk)

# Compute ll position on surface directly under target/orbiter
target_ll = pm.aer2geodetic(az_output, el_output, range_output * 1000, lat, lng, height, None, True)
Expand All @@ -141,86 +109,23 @@ def ll2aerll(lng, lat, height, target, time):
"longitude": target_ll[1],
"latitude": target_ll[0],
"altitude": target_ll[2],
"horizontal_altitude": horizontal_altitude
"horizontal_altitude": horizontal_altitude,
"ancillary": {
"sun_az": sun_az_output,
"sun_el": sun_el_output,
"earth_az": earth_az_output,
"earth_el": earth_el_output
}
})

"""
def lltarg2vistimes(lng, lat, height, target, startTime):
package_dir = os.path.dirname(os.path.abspath(__file__))
# Set vars
MAXWIN = 750
TIMFMT = "YYYY-MON-DD HR:MN:SC.###### (TDB) ::TDB ::RND"
TIMLEN = 41
GS_RESULT = spiceypy.stypes.SPICEDOUBLE_CELL(2 * MAXWIN)
cnfine = spiceypy.stypes.SPICEDOUBLE_CELL(2)
relate = ">"
crdsys = "LATITUDINAL"
coord = "LATITUDE"
targ = "MRO"
obsrvr = "-654321"
frame = "IAU_MARS"
abcorr = "NONE"
#Store the time bounds of our search interval in
#the cnfine confinement window.
begtim = spiceypy.str2et( "2023 MAY 01" )
endtim = spiceypy.str2et( "2023 MAY 05" )
spiceypy.wninsd( begtim, endtim, cnfine )
#This search uses a step size of four hours since the
#time for all declination zero-to-max-to-zero passes
#within the search window exceeds eight hours.
#The example uses an 83 degree elevation because of its
#rare occurrence and short duration.
step = int((4.0/24.0) * spiceypy.spd())
adjust = 0.0
refval = 65.0 * spiceypy.rpd()
#List the beginning and ending points in each interval
#if result contains data.
spiceypy.gfposc( targ, frame, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step, MAXWIN, cnfine, GS_RESULT )
spiceypy.wnfetd(GS_RESULT, 0)
#Creates an array of times in which it passes over target
GS_times = []
The_times = [GS_times]
The_results = [GS_RESULT]
for i in range(len(The_times)):
for ii in range(spiceypy.wncard(The_results[i])):
The_times[i].append(np.linspace(*spiceypy.wnfetd(The_results[i], ii), step, endpoint=False))
visibilities = []
for time in GS_times:
t = spiceypy.timout(time, TIMFMT)
s = time[0]
st = time[-1]
duration = (st -s)
start = t[0]
stop = t[-1]
closest_approach = spiceypy.timout((st+s)/2, TIMFMT)
visibilities.append({"start": start, "stop": stop, "closest_approach": closest_approach})
return visibilities
"""

lng = float(sys.argv[1])
lat = float(sys.argv[2])
height = float(sys.argv[3])
target = unquote(sys.argv[4])
time = unquote(sys.argv[5])
includeSunEarth = sys.argv[6]

try:
print(ll2aerll(lng, lat, height, target, time))
print(ll2aerll(lng, lat, height, target, time, includeSunEarth))
except:
print(json.dumps({"error": True, "message": 'Error: ' + str(sys.exc_info()[0])}))
Binary file removed private/api/spice/pinpoint.exe
Binary file not shown.
Loading

0 comments on commit e1cb586

Please sign in to comment.