def toImage()

in analysis/webservice/algorithms/Tomogram3D.py [0:0]

• 104 lines of code
• 3 McCabe index (conditional complexity)

    def toImage(self):
        _, _, _, view_azim, view_elev = 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=view_elev, azim=view_azim)

        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')
            s = 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()

        logger.info('Writing plot to buffer')
        plt.savefig(buffer, format='png', facecolor='white')

        buffer.seek(0)
        plt.close(fig)
        return buffer.read()