in analysis/webservice/algorithms/Tomogram3D.py [0:0]
def toGif(self):
orbit_elev, orbit_step, frame_duration, _, _ = self.render_params
results = self.results()
xyz = results[['lon', 'lat', 'elevation']].values
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev=orbit_elev, azim=0)
if not self.hide_basemap:
min_lat, min_lon, max_lat, max_lon = self.bounds
# Originally used mpl_toolkits.basemap.Basemap.aspect for this, but basemap installs an LGPL3 data package
# so we can't have that in a release. Thankfully, manual computation ends up being quite simple
aspect = (max_lat - min_lat) / (max_lon - min_lon)
basemap_size = 512
params = dict(
bbox=f'{min_lon},{min_lat},{max_lon},{max_lat}',
bboxSR=4326, imageSR=4326,
size=f'{basemap_size},{int(aspect * basemap_size)}',
dpi=2000,
format='png32',
transparent=True,
f='image'
)
url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/export'
logger.info(f'Pulling basemap ({url}?{urlencode(params)})')
try:
elevations = results[['elevation']].values
z_range = np.nanmax(elevations) - np.nanmin(elevations)
z_coord = np.nanmin(elevations) - (0.1 * z_range)
r = requests.get(url, params=params)
r.raise_for_status()
buf = BytesIO(r.content)
img = Image.open(buf)
img_data = np.flipud(np.array(img))
lats = np.linspace(min_lat, max_lat, num=img.height)
lons = np.linspace(min_lon, max_lon, num=img.width)
X, Y = np.meshgrid(lons, lats)
Z = np.full(X.shape, z_coord)
logger.info('Plotting basemap')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=img_data / 255)
except:
logger.error('Failed to pull basemap, will not draw it')
if self.render_type == 'scatter':
logger.info('Plotting tomogram data')
if self.style == 'db':
v = dict(vmin=-30, vmax=-10)
cbl = 'Tomogram (dB)'
elif self.style == 'percentile':
v = dict(vmin=0, vmax=1)
cbl = 'Tomogram Cumulative value (%/100)'
else:
elevs = results[['elevation']].values
v = dict(vmin=np.nanmin(elevs), vmax=np.nanmax(elevs))
cbl = 'Voxel elevation'
s = ax.scatter(
xyz[:, 0], xyz[:, 1], xyz[:, 2],
marker=',',
c=results[['tomo_value']].values if self.style != 'percentile_elevation'
else results[['elevation']].values,
zdir='z',
depthshade=True,
cmap=self.cmap,
**v
)
cbar = fig.colorbar(s, ax=ax)
cbar.set_label(cbl)
else:
logger.info('Collecting data by lat/lon')
lats = np.unique(results['lat'].values)
lons = np.unique(results['lon'].values)
data_dict = {}
for r in results.itertuples(index=False):
key = (r.lon, r.lat)
if key not in data_dict:
data_dict[key] = ([r.elevation], [r.tomo_value])
else:
data_dict[key][0].append(r.elevation)
data_dict[key][1].append(r.tomo_value)
logger.info('Determining peaks')
vals = np.empty((len(lats), len(lons)))
for i, lat in enumerate(lats):
for j, lon in enumerate(lons):
if (lon, lat) in data_dict:
elevs, tomo_vals = data_dict[(lon, lat)]
i_max = np.argmax(np.array(tomo_vals))
vals[i, j] = elevs[i_max]
else:
vals[i, j] = np.nan
X2, Y2 = np.meshgrid(lons, lats)
logger.info('Plotting peak surface')
ax.plot_surface(
X2,
Y2,
vals,
rstride=1, cstride=1,
color='xkcd:leaf'
)
ax.set_ylabel('Latitude')
ax.set_xlabel('Longitude')
ax.set_zlabel('Elevation w.r.t. dataset reference (m)')
plt.tight_layout()
buffer = BytesIO()
with TemporaryDirectory() as td:
for azim in range(0, 360, orbit_step):
logger.info(f'Saving frame for azimuth = {azim}')
ax.view_init(azim=azim)
plt.savefig(join(td, f'fr_{azim}.png'))
with contextlib.ExitStack() as stack:
logger.info(f'Combining frames into final GIF')
imgs = (stack.enter_context(Image.open(join(td, f'fr_{a}.png'))) for a in range(0, 360, orbit_step))
img = next(imgs)
img.save(buffer, format='GIF', append_images=imgs, save_all=True, duration=frame_duration, loop=0)
buffer.seek(0)
plt.close(fig)
return buffer.read()