Commit ad359acf authored by Eva Schiffer's avatar Eva Schiffer
Browse files

a first draft at moving to cartopy

parent 30b811f6
This diff is collapsed.
#!/usr/bin/env python
# encoding: utf-8
"""
$Id: graphics.py,v 1.3 2008/05/08 16:46:02 rayg Exp $
Library of routines for plotting HSR data on a map
This file was completely rewritten in 2019 to use cartopy objects (instead of basemap)
to plot data on the Earth.
This file is based on the Keoni.map.graphics but has a significantly changed
interface to support more modular interactions with the basemap object.
~ Eva Schiffer, August 20th, 2009
~ Eva Schiffer, Aug. 2019
"""
from mpl_toolkits.basemap import Basemap, shiftgrid
from numpy import arange, array, reshape, concatenate, nan
import logging
# the value that will denote "bad" longitudes and latitudes
badLonLat = 1.0E30
from cartopy import crs
import cartopy.feature as cfeature
def create_basemap (lon, lat=None, axis=None, projection='lcc', resolution='i') :
import numpy
from numpy import argmax, ma
LOG = logging.getLogger(__name__)
def _cap_extent (my_extent, projection_object, ) :
"""
Create an instance of basemap using either the specified axis info or the
specified lon and lat info to pick the viewing area.
the format of the axis is [lon min, lon max, lat min, lat max]
where the min and max are for the entire area that you wish to show.
Note: There are known viewing area problems with conic projections that
may cause "rectangular" data to be clipped.
Make sure my extent doesn't go beyond the available extent for the projection.
"""
if lat is None:
# then assume (lon,lat) pairs in a generator
perimeter = list(lon)
lon = array([x[0] for x in perimeter])
lat = array([x[1] for x in perimeter])
# make sure the axis is the correct shape
# and fill in any missing values
if axis is None:
axis = [None,None,None,None]
if axis[0] is None:
axis[0] = lon.min()
if axis[1] is None:
axis[1] = lon.max()
if axis[2] is None:
axis[2] = lat.min()
if axis[3] is None:
axis[3] = lat.max()
# pull out the longitude/latitude info
lon_left = axis[0]
lat_bottom = axis[2]
lon_right = axis[1]
lat_top = axis[3]
lon_mid = (lon_left + lon_right ) / 2.
lat_mid = (lat_top + lat_bottom) / 2.
# make our basemap
m = None
if projection is 'ortho' :
# orthographic projections require this call
m = Basemap(resolution=resolution, area_thresh=10000., projection=projection,
lat_0=lat_mid, lon_0=lon_mid)
else :
# most of the other projections use this call
m = Basemap(llcrnrlon=lon_left,llcrnrlat=lat_bottom,urcrnrlon=lon_right,urcrnrlat=lat_top,
resolution=resolution, area_thresh=10000., projection=projection,
lat_1=lat_mid,lon_0=lon_mid)
return m, axis
def draw_basic_features(baseMapInstance, axis) :
global_bounds = projection_object.boundary.bounds
g_extent = (global_bounds[0], global_bounds[2], global_bounds[1], global_bounds[3])
return [max(my_extent[0], g_extent[0]), min(my_extent[1], g_extent[1]), max(my_extent[2], g_extent[2]), min(my_extent[3], g_extent[3]),]
def make_or_find_good_mask(data_array, good_mask=None,) :
"""
Draw the basic outlines of the earth's features.
If we don't have a good mask, then try to find or make one.
"""
# draw the basic physical and geopolitical features
baseMapInstance.drawcoastlines()
baseMapInstance.drawcountries()
baseMapInstance.drawstates()
baseMapInstance.drawmapboundary()
# pull out the longitude/latitude info
lon_left = axis[0]
lat_bottom = axis[2]
lon_right = axis[1]
lat_top = axis[3]
# draw the parallels and meridians
parallels = arange(-80.,90.,abs(lat_top - lat_bottom) / 4.0)
baseMapInstance.drawparallels(parallels,labels=[1,0,0,1])
meridians = arange(0., 360.,abs(lon_left - lon_right) / 4.0)
baseMapInstance.drawmeridians(meridians,labels=[1,0,0,1])
return
def show_lon_lat_data(lon, lat, baseMapInstance, data=None, levelsToUse=None, **kwargs) :
# if there is no good mask, find one
if good_mask is None :
# check to see if we can get a mask from a masked array
if data_array is not None and type(data_array) is ma :
good_mask = data_array.mask
# if we still don't have a mask, make a simple one
if good_mask is None :
good_mask = numpy.ones(data_array.shape, dtype=numpy.bool,)
return good_mask
def lon_subtract_360(lon1, lon2) :
"""
Show data corresponding to the longitude and latitude set provided on the earth using the provided basemap.
levelsToUse is a list of numbers representing data ranges that will be used
return lon1 - lon2, correctly handling the wraparound case
inputs can be in [-180, 180) form or [0, 360) form
output is in the [0, 360) range
"""
x, y = baseMapInstance(lon, lat) # translate into the coordinate system of the basemap
return show_x_y_data(x, y, baseMapInstance, data, levelsToUse, **kwargs)
return (lon1 - lon2) % 360.
def show_x_y_data(x, y, baseMapInstance, data=None, levelsToUse=None, **kwargs) :
def lon_add_180(lon1, lon2) :
"""
Show data corresponding to a given x, y using the provided basemap.
levelsToUse is a list of numbers representing data ranges that will be used
return lon1 + lon2, correctly handling the wraparound case
inputs can be in [-180, 180) form or [0, 360) form
output is in the [-180, 180) range
"""
artistsAdded = [ ]
# only try to plot the data if there is some
if data is not None:
newX = make_2D_version_if_needed(x, badLonLat)
newY = make_2D_version_if_needed(y, badLonLat)
newData = make_2D_version_if_needed(data, nan)
if levelsToUse is not None :
p = baseMapInstance.contourf(newX, newY, newData, levelsToUse, **kwargs)
else :
p = baseMapInstance.contourf(newX, newY, newData, **kwargs)
return (lon1 + lon2 + 180.) % 360. - 180.
def calculate_lon_mid(lon_left, lon_right) :
"""
calculate midpoint between two longitudes
"""
return lon_add_180(lon_left, (lon_subtract_360(lon_right, lon_left) / 2.))
def calculate_lon_range(longitude_data, good_mask, longitude_data_b=None, good_mask_b=None, ) :
"""
calculate range of longitudes.
Can handle data that spans the anti-meridian
returns the westmost and eastmost longitude that define the range, in that order
Note: This assumes that the ranges of the main data and b data are mostly overlapping.
This strategy will fail when the two data sets are disparate in longitude, but if that much
data is missing even though the base lon/lat data matches Glance is going to have some
problems drawing your image correctly anyway.
Note: This assumes that data you pass in has a valid good mask.
"""
# select only the good data values and flatten out our data for ease of use
good_data = longitude_data[good_mask].flatten()
if longitude_data_b is not None :
numpy.append(good_data, longitude_data_b[good_mask_b].flatten())
# if we have nothing, just stop here
if good_data.size == 0 :
return None, None
# sort the lons, smallest to largest
good_data.sort()
# find the largest gap in lons. OK to subtract lons because values are sorted from -180 to 180
lonDiffs = good_data[1:] - good_data[:-1]
i = argmax(lonDiffs)
lonMax = good_data[i]
lonMin = good_data[i+1]
LOG.debug("naive gap detected: lonMin, lonMax, gap = " + str(lonMin) + ", " + str(lonMax) + ", " + str(lonMin - lonMax))
# check for largest gap spanning anti-meridian
if lon_subtract_360(good_data[0], good_data[-1]) > lon_subtract_360(lonMin, lonMax):
lonMin, lonMax = good_data[0], good_data[-1]
LOG.debug("anti-meridian gap detected: lonMin, lonMax, gap = " + str(lonMin) + ", " + str(lonMax) + ", " + str(lonMin - lonMax))
return lonMin, lonMax
def calculate_lat_range(latitude_data, good_mask, latitude_data_b=None, good_mask_b=None,) :
"""
Calculate the range of latitudes.
If a B data set is given, include the range of the latitude_data and the latitude_data_B.
returns the min latitude and max latitude that define the range in that order
Note: This assumes that data you pass in has a valid good mask.
"""
# select only the good data values and flatten out our data for ease of use
good_data = latitude_data[good_mask].flatten()
if latitude_data_b is not None:
numpy.append(good_data, latitude_data_b[good_mask_b].flatten())
# if we have nothing, just stop here
if good_data.size == 0:
return None, None
return numpy.min(good_data), numpy.max(good_data)
def get_extents(longitude_data, latitude_data,
lon_good_mask=None, lat_good_mask=None,
longitude_data_b=None, latitude_data_b=None,
lon_good_mask_b=None, lat_good_mask_b=None,) :
"""
Calculate the extents that will include all this data.
returns extents in the format [lon westmost, lon eastmost, lat northmost, lat southmost]
Note: The b data is entirely optional, but masks must be included
for the data provided, eithere as explicit numpy arrays or by making
the original lat/lon masked arrays instead of simple numpy arrays.
"""
print ("lat type: " + str(type(latitude_data)))
print ("lon type: " + str(type(longitude_data)))
if latitude_data is None:
# then assume (lon,lat) pairs in a generator
perimeter = list(longitude_data)
longitude_data = numpy.array([x[0] for x in perimeter])
latitude_data = numpy.array([x[1] for x in perimeter])
if longitude_data_b is not None and latitude_data_b is None:
# again, assume this is a generator
perimeter_b = list(longitude_data_b)
longitude_data_b = numpy.array([x[0] for x in perimeter_b])
latitude_data_b = numpy.array([x[1] for x in perimeter_b])
# make sure we aren't processing lon or lat values that are way outside the acceptable range
lat_valid_mask = (latitude_data >= -90.0) & (latitude_data <= 90.0)
lon_valid_mask = (longitude_data >= -180.0) & (longitude_data <= 360.0)
lat_valid_mask_b = (latitude_data_b >= -90.0) & (latitude_data_b <= 90.0) if latitude_data_b is not None else None
lon_valid_mask_b = (longitude_data_b >= -180.0) & (longitude_data_b <= 360.0) if longitude_data_b is not None else None
# double check that all our good masks exist
lat_good_mask = make_or_find_good_mask(latitude_data, good_mask=lat_good_mask,)
lon_good_mask = make_or_find_good_mask(longitude_data, good_mask=lon_good_mask,)
lat_good_mask_b = make_or_find_good_mask(latitude_data_b, good_mask=lat_good_mask_b,) if latitude_data_b is not None else None
lon_good_mask_b = make_or_find_good_mask(longitude_data_b, good_mask=lon_good_mask_b,) if longitude_data_b is not None else None
# make some combined masks
a_mask = lon_good_mask & lon_valid_mask & lat_good_mask & lat_valid_mask
b_mask = lon_good_mask_b & lon_valid_mask_b & lat_good_mask_b & lat_valid_mask_b if latitude_data_b is not None else None
# get the lon and lat bounds for our extents
lon_westmost, lon_eastmost = calculate_lon_range(longitude_data, good_mask=a_mask,
longitude_data_b=longitude_data_b, good_mask_b=b_mask,)
lat_northmost, lat_southmost = calculate_lat_range(latitude_data, good_mask=a_mask,
latitude_data_b=latitude_data_b, good_mask_b=b_mask,)
return [lon_westmost, lon_eastmost, lat_northmost, lat_southmost,]
def create_cartopy_proj_objects (extents, out_proj4_str=None,) :
"""
Create cartopy "crs" projection objects for the input data and output projections.
This is only a first draft of how Glance will probably use projections in the future.
For now, we will assume that all data comes in with lat/lon and therefore
corresponds to a PlateCarree crs object.
The format of the extents is [lon westmost, lon eastmost, lat northmost, lat southmost]
where these extents define the whole area you wish to show.
If the out projection string is None we will be using PlateCarree as our output projection.
TODO, other projections are not currently supported, but should be someday.
Note: You will need to call get_extents on your own now. And the extents will be in lat/lon,
so you will need to use the input_proj_obj when you're setting your plot extents later.
"""
LOG.debug("Using extents: " + str(extents))
if out_proj4_str is not None :
LOG.debug("Input proj4 string provided, but this code cannot use it yet: " + out_proj4_str)
# FUTURE, this is not actually useful here now, but hang on to it for bounds setting later in our display process
# # pull out the longitude/latitude info
# lon_left = extents[0]
# lon_right = extents[1]
# lat_bottom = extents[2]
# lat_top = extents[3]
# lon_mid = calculate_lon_mid(lon_left, lon_right) # handle lon wraparound case
# lat_mid = (lat_top + lat_bottom) / 2.
# for now just use PlateCarree, in
# FUTURE we will want to pick projections based on our input so we could use x and y when desired
input_proj_obj = crs.PlateCarree()
# FUTURE, we will probably want different projections; select based on the lat/lon range?
# for reference, here's what cartopy supports: https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html#cartopy-projections
output_proj_obj = crs.PlateCarree()
return input_proj_obj, output_proj_obj
def draw_basic_features(in_proj_obj, axes_object, extents,) :
"""Draw the basic outlines of the earth's features.
"""
# set the extents (Note: we may have problems setting extents while plotting full disks TODO test that)
axes_object.set_extent(extents, crs=in_proj_obj)
# draw some features on our map
line_color = 'black' # Glance always uses white background and black lines
# draw some lat/lon gridlines
axes_object.gridlines(linestyle=':', zorder=2, edgecolor=line_color, )
# TODO, do we need to add anything to label the gridlines?
# draw physical and geopolitical stuff
axes_object.coastlines(resolution='50m', zorder=3, color=line_color, )
axes_object.add_feature(cfeature.BORDERS, linewidth=0.5, zorder=4, edgecolor=line_color, ) # this makes country borders
axes_object.add_feature(cfeature.LAKES, linestyle=':', linewidth=0.5, facecolor='none', edgecolor=line_color, zorder=5, )
# add states and provinces
# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
axes_object.add_feature(states_provinces, edgecolor=line_color, linewidth=0.25, zorder=6, )
# return the original x and y so the caller can match any external data in shape
return baseMapInstance, x, y
return
def contourf_on_map(lon, lat, data, axes_object, in_projection, levelsToUse=None, **kwargs ) :
"""
Contour plot the data on a map. Assume the input coordinates are in the in_projection form
and the output should be plotted in out_projection.
optional: levelsToUse is a list of numbers representing data ranges that will be used
"""
if levelsToUse is not None:
img_temp = axes_object.contourf(lon, lat, data, levelsToUse, transform=in_projection, **kwargs)
else:
img_temp = axes_object.contourf(lon, lat, data, transform=in_projection, **kwargs)
return img_temp
def show_quiver_plot (lon, lat, baseMapInstance, (uData, vData)=(None,None), colordata=None, **kwargs) :
def quiver_plot_on_map(lon, lat, axes_object, in_projection, (uData, vData)=(None, None), colordata=None, **kwargs) :
"""
Show a quiver plot of the given vector data at the given longitude and latitude
"""
x, y = baseMapInstance(lon, lat) # translate into the coordinate system of the basemap
# show the quiver plot if there is data
if (uData is not None) and (vData is not None) :
if (uData is not None) and (vData is not None):
if colordata is None:
p = baseMapInstance.quiver(x, y, uData, vData, **kwargs)
else :
p = baseMapInstance.quiver(x, y, uData, vData, colordata, **kwargs)
# return the original x and y so the caller can match any external data in shape
return baseMapInstance, x, y
img_temp = axes_object.quiver(lon, lat, uData, vData, transform=in_projection, **kwargs)
else:
img_temp = axes_object.quiver(lon, lat, uData, vData, colordata, transform=in_projection, **kwargs)
def make_2D_version_if_needed(originalArray, fillValue) :
return img_temp
def plot_on_map(lon, lat, axes_object, in_projection, marker_const='o', color_const='#000000', marker_size=1, **kwargs) :
"""
if the array isn't already 2D, add some fillValued points and reshape it
to make it so
Show a simple point plot at the given longitude and latitude.
"""
arrayToChange = originalArray
# if the array has too few dimentions, we'll need to reshape it
if len(arrayToChange.shape) < 2 :
# if there are an odd number of points, add a null point
if (arrayToChange.size % 2) > 0 :
arrayToChange = concatenate((arrayToChange, [fillValue]))
# now fold the array in half
newDim = arrayToChange.size / 2
arrayToChange = reshape(arrayToChange, (newDim, 2))
return arrayToChange
def mapshow(lon, lat=None, data=None, axis=None, projection='lcc',resolution='i', levelsToUse=None, **kwargs):
"""Like matshow(), but requiring longitude and latitude
optional axis is in the form of [left-longitude, right-longitude, bottom-latitude, top-latitude]
axis values left as None will be auto-maximized
returns basemap object, x-transform, y-transform
if lat is None, then assume that lon is a sequence of (lon,lat) points
"""
# make the basemap object
baseMapInstance, axis = create_basemap(lon, lat, axis, projection, resolution)
# draw the minimal map features
draw_basic_features(baseMapInstance, axis)
# plot the actual data on the map
baseMapInstance, x, y = show_lon_lat_data(lon, lat, baseMapInstance, data, levelsToUse, **kwargs)
return baseMapInstance, x, y
# show the plot if there is data
if (lon.size > 0) and (lat.size > 0) and (lon.size == lat.size) :
img_temp = axes_object.plot(lon, lat, marker_const, color=color_const, markersize=marker_size, transform=in_projection, **kwargs)
return img_temp
......@@ -7,7 +7,7 @@ Created by evas Oct 2011.
Copyright (c) 2011 University of Wisconsin SSEC. All rights reserved.
"""
import sys, os.path, logging
import sys, logging, traceback
from PyQt5.QtWidgets import QApplication
......@@ -52,7 +52,7 @@ class GlanceGUIController (object) :
self.figs = gui_figs.GlanceGUIFigures(self.model)
# how about an icon?
from pkg_resources import resource_filename;
from pkg_resources import resource_filename
from PyQt5.QtGui import QIcon
self.qtApp.setWindowIcon(QIcon(resource_filename(__name__, "glance-icon.png")))
......@@ -73,7 +73,7 @@ class GlanceGUIController (object) :
start up the GUI
"""
self.view.show()
return(self.qtApp.exec_()) # this is needed to keep the gui from exiting prematurely
return self.qtApp.exec_() # this is needed to keep the gui from exiting prematurely
################# methods to handle user input reporting #################
......@@ -247,7 +247,7 @@ class GlanceGUIController (object) :
"""
the user has asked for a plot
"""
try :
self.figs.spawnPlot()
except (IncompatableDataObjects, ValueError) as idove :
......
......@@ -25,9 +25,8 @@ import numpy as np
import glance.data as dataobjects
import glance.figures as figures
import glance.gui_model as model
import glance.graphics as maps
from glance.gui_constants import *
from glance.plotcreatefns import select_projection
LOG = logging.getLogger(__name__)
......@@ -65,17 +64,17 @@ CAN_BE_MAPPED = {
# which data sets the plot needs
NEEDED_DATA_PER_PLOT = \
{
ORIGINAL_A : set([A_CONST]),
ORIGINAL_B : set([ B_CONST]),
ABS_DIFF : set([A_CONST, B_CONST]),
RAW_DIFF : set([A_CONST, B_CONST]),
HISTOGRAM_A : set([A_CONST]),
HISTOGRAM_B : set([ B_CONST]),
HISTOGRAM : set([A_CONST, B_CONST]),
MISMATCH : set([A_CONST, B_CONST]),
SCATTER : set([A_CONST, B_CONST]),
D_SCATTER : set([A_CONST, B_CONST]),
HEX_PLOT : set([A_CONST, B_CONST]),
ORIGINAL_A : {A_CONST, },
ORIGINAL_B : { B_CONST,},
ABS_DIFF : {A_CONST, B_CONST,},
RAW_DIFF : {A_CONST, B_CONST,},
HISTOGRAM_A : {A_CONST, },
HISTOGRAM_B : { B_CONST,},
HISTOGRAM : {A_CONST, B_CONST,},
MISMATCH : {A_CONST, B_CONST,},
SCATTER : {A_CONST, B_CONST,},
D_SCATTER : {A_CONST, B_CONST,},
HEX_PLOT : {A_CONST, B_CONST,},
}
class GlanceGUIFigures (object) :
......@@ -147,7 +146,7 @@ class GlanceGUIFigures (object) :
# get the red, green, and blue data
canGetData = self.dataModel.makeSureVariablesAreAvailable(filePrefix, [RED_VAR_NAME, GREEN_VAR_NAME, BLUE_VAR_NAME])
if not canGetData : # if the basic rgb data doesn't exist, stop now
"", None, ""
return "", None, ""
_, rDataObj, _ = self._getVariableInformation(filePrefix, variableName=RED_VAR_NAME, doCorrections=False)
_, gDataObj, _ = self._getVariableInformation(filePrefix, variableName=GREEN_VAR_NAME, doCorrections=False)
_, bDataObj, _ = self._getVariableInformation(filePrefix, variableName=BLUE_VAR_NAME, doCorrections=False)
......@@ -171,7 +170,7 @@ class GlanceGUIFigures (object) :
"""
if appropriate for the image type, build the difference object, otherwise return None
this method may rase an IncompatableDataObjects exception if the two data objects it's given can't be compared
this method may raise an IncompatableDataObjects exception if the two data objects it's given can't be compared
"""
diffObject = None
......@@ -188,9 +187,9 @@ class GlanceGUIFigures (object) :
return diffObject
def _load_and_analyse_lonlat (self, listOfFilePrefixes=[A_CONST, B_CONST], lonNames=None, latNames=None, stopIfComparisonFails=False) :
def _load_and_analyse_lonlat (self, listOfFilePrefixes=list({A_CONST, B_CONST,}), lonNames=None, latNames=None, ) :
"""
load information on the longidue and latitude,
load information on the longitude and latitude,
if there are multiple file prefixes given:
find the shared range
analyse how different the navigation is between the files
......@@ -199,66 +198,83 @@ class GlanceGUIFigures (object) :
lonNames and latNames should be dictionaries giving the names of the longitude and latitude variables indexed by the file prefixes
This method may raise an IncompatableDataObjects exception if multiple file prefixes are passed in the listOfFilePrefixes
and the longitude and latidues for those files can not be compared.
and the longitude and latitudes for those files can not be compared.
"""
lonlatData = { }
lonRange = None
latRange = None
# load and compare stuff for each file prefix
for filePrefix in listOfFilePrefixes :
extents = None
# this now explicitly only works for 1 or 2 data sets
if len(listOfFilePrefixes) > 2 or len(listOfFilePrefixes) <= 0 :
LOG.debug("Somehow you've gotten into a state where you're trying to analyze lonlat for a meaningless number of files.")
return { }, None, None
# load and process stuff for each file prefix
for filePrefix in listOfFilePrefixes:
# get information on the lon/lat from the current file
currentLonObj, currentLatObj, currentLonRange, currentLatRange = self._load_lonlat(filePrefix, lonNames[filePrefix], latNames[filePrefix])
currentLonObj, currentLatObj = self._load_lonlat(filePrefix, lonNames[filePrefix], latNames[filePrefix])
currentLonObj.self_analysis()
currentLatObj.self_analysis()
# expand our lon/lat ranges if we need to
if lonRange is None :
lonRange = currentLonRange
else :
lonRange[0] = min(currentLonRange[0], lonRange[0])
lonRange[1] = max(currentLonRange[1], lonRange[1])
if latRange is None:
latRange = currentLatRange
else :
latRange[0] = min(currentLatRange[0], latRange[0])
latRange[1] = max(currentLatRange[1], latRange[1])
# we can't use longitude and latitude that don't match in size
if currentLonObj.data.shape != currentLatObj.data.shape :
raise ValueError ("Longitude and Latitude for file " + filePrefix + " are different shapes." +
"\nCannot match differently shaped navigation data.")
# compare this file to whatever other data we have