diff --git a/clients/python/sliderule/io.py b/clients/python/sliderule/io.py index 91f19133f..4268120a4 100644 --- a/clients/python/sliderule/io.py +++ b/clients/python/sliderule/io.py @@ -31,15 +31,20 @@ import json import warnings import datetime -import geopandas import numpy as np -from sliderule import logger +from sliderule import logger, SLIDERULE_EPSG # imports with warnings if not present try: import scipy.io except ModuleNotFoundError as e: sys.stderr.write("Warning: missing packages, some functions will throw an exception if called. (%s)\n" % (str(e))) +try: + import geopandas + from geopandas.io.arrow import _geopandas_to_arrow + from geopandas._compat import import_optional_dependency +except ModuleNotFoundError as e: + sys.stderr.write("Warning: missing packages, some functions will throw an exception if called. (%s)\n" % (str(e))) try: import h5py except ModuleNotFoundError as e: @@ -337,6 +342,25 @@ def winding(x,y): wind = np.sum([(x[i+1] - x[i])*(y[i+1] + y[i]) for i in range(npts - 1)]) return wind +# geodesic area of a polygon in square kilometers +def area(lon, lat, a=6378.1370): + assert len(lon) == len(lat), "Longitude and Latitude arrays must be the same length" + assert len(lon) > 2, "Need at least 3 points to calculate area" + # factor for converting to radians + dtr = np.pi/180.0 + # initialize area + area = 0.0 + npts = len(lon) + # for each point in the polygon + for i in range(npts): + # wrap coordinates + ii = (i + 1) % npts + # calculate geodesic area in square kilometers + area += (a**2)*(lon[ii]*dtr - lon[i]*dtr) * \ + (2.0 + np.sin(lat[i]*dtr) + np.sin(lat[ii]*dtr))/2.0 + # return absolute value of area + return np.abs(area).squeeze() + # fix longitudes to be -180:180 def wrap_longitudes(lon): phi = np.arctan2(np.sin(lon*np.pi/180.0),np.cos(lon*np.pi/180.0)) @@ -361,7 +385,7 @@ def from_region(polygon): # convert geodataframe vector object to a list of sliderule regions def from_geodataframe(gdf): # verify that geodataframe is in latitude/longitude - geodataframe = gdf.to_crs(epsg=4326) + geodataframe = gdf.to_crs(SLIDERULE_EPSG) # create a list of regions regions = [] # for each region @@ -375,7 +399,11 @@ def to_json(filename, **kwargs): # set default keyword arguments kwargs.setdefault('parameters',None) kwargs.setdefault('regions',[]) - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) + # import optional dependencies + pyproj = import_optional_dependency( + "pyproj", extra="pyproj is required for CRS conversion." + ) # add each parameter as an attribute SRparams = ['H_min_win', 'atl08_class', 'atl03_quality', 'ats', 'cnf', 'cnt', 'len', 'maxi', 'res', 'sigma_r_max', 'srt', 'yapc'] @@ -388,7 +416,7 @@ def to_json(filename, **kwargs): except: pass # save CRS to JSON - crs = geopandas.tools.crs.CRS.from_string(kwargs['crs']) + crs = pyproj.CRS.from_string(kwargs['crs']) output['crs'] = crs.to_string() # save each region following GeoJSON specification output['type'] = 'FeatureCollection' @@ -437,12 +465,115 @@ def from_json(filename, **kwargs): # return the sliderule parameters and regions return (parms, regions) +def to_parquet(gdf, filename, **kwargs): + # set default keyword arguments + kwargs.setdefault('index',None) + kwargs.setdefault('compression','snappy') + kwargs.setdefault('schema_version',None) + kwargs.setdefault('parameters',dict()) + kwargs.setdefault('regions',[]) + kwargs.setdefault('crs',SLIDERULE_EPSG) + # import optional dependencies + pyproj = import_optional_dependency( + "pyproj", extra="pyproj is required for CRS conversion." + ) + parquet = import_optional_dependency( + "pyarrow.parquet", extra="pyarrow is required for Parquet support." + ) + # output metadata + output = {} + # for each adjustable sliderule parameter + [kwargs['parameters'].pop(p) for p in ['poly']] + for p,val in kwargs['parameters'].items(): + # try to convert the parameter if available + try: + output[p] = attributes_encoder(val) + except: + pass + # save CRS to JSON + crs = pyproj.CRS.from_string(kwargs['crs']) + output['crs'] = crs.to_string() + # save each region following GeoJSON specification + output['type'] = 'FeatureCollection' + output['features'] = [] + for i,poly in enumerate(kwargs['regions']): + lon, lat = from_region(poly) + lon = attributes_encoder(lon) + lat = attributes_encoder(lat) + geometry=dict(type="polygon", coordinates=[]) + geometry['coordinates'].append([[ln,lt] for ln,lt in zip(lon,lat)]) + output['features'].append(dict(type="Feature", geometry=geometry)) + # dump the attributes to encoded JSON-format + sliderule_metadata = json.dumps(output).encode('utf-8') + # convert geodataframe to arrow table + table = _geopandas_to_arrow(gdf, + index=kwargs['index'], + schema_version=kwargs['schema_version'] + ) + # store sliderule specific file-level metadata + metadata = table.schema.metadata + metadata.update({b"sliderule": sliderule_metadata}) + # replace schema metadata with updated + table = table.replace_schema_metadata(metadata) + # write arrow table to parquet file + parquet.write_table(table, filename, + compression=kwargs['compression'] + ) + +def from_parquet(filename, **kwargs): + # set default keyword arguments + kwargs.setdefault('crs',SLIDERULE_EPSG) + kwargs.setdefault('return_parameters',False) + kwargs.setdefault('return_regions',False) + # import optional dependencies + parquet = import_optional_dependency( + "pyarrow.parquet", extra="pyarrow is required for Parquet support." + ) + # read arrow table from parquet file + # and convert to coordinate reference system + gdf = geopandas.read_parquet(filename).to_crs(kwargs['crs']) + # if not returning the query parameters or polygon + if not (kwargs['return_parameters'] or kwargs['return_regions']): + # return geodataframe + return gdf + # create tuple with returns + output = (gdf,) + # get parquet file metadata + metadata = parquet.read_metadata(filename).metadata + # validate sliderule metadata + if b'sliderule' not in metadata.keys(): + logger.error("No sliderule metadata found in Parquet file") + return output + # decode sliderule metadata from JSON + parms = json.loads(metadata[b'sliderule'].decode('utf-8')) + # create a list of regions + regions = [] + # for each feature in the JSON file + for feature in parms.pop('features'): + # for each coordinate set in the feature + for coords in feature['geometry']['coordinates']: + # append to sliderule regions + regions.append([{'lon':ln,'lat':lt} for ln,lt in coords]) + # if returning the parameters + if kwargs['return_parameters']: + # add parameters to output tuple + [parms.pop(v) for v in ('version','commit','type')] + output += (parms,) + # if returning the regions + if kwargs['return_regions']: + # add regions to output tuple + output += (regions,) + # return the combined tuple + return output + # output geodataframe to netCDF (version 3) def to_nc(gdf, filename, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default keyword arguments kwargs.setdefault('parameters',None) kwargs.setdefault('regions',[]) - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('units','seconds since 2018-01-01T00:00:00') @@ -492,7 +623,7 @@ def to_nc(gdf, filename, **kwargs): # save geodataframe coordinate system fileID.crs = kwargs['crs'] # add geospatial attributes - if (kwargs['crs'] == 'EPSG:4326'): + if kwargs['crs'] in ('EPSG:4326', SLIDERULE_EPSG): fileID.geospatial_lat_units = \ attributes['geospatial_lat_units'] fileID.geospatial_lon_units = \ @@ -535,8 +666,10 @@ def to_nc(gdf, filename, **kwargs): # input geodataframe from netCDF (version 3) def from_nc(filename, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default crs - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('index_key','time') @@ -631,7 +764,7 @@ def to_hdf(gdf, filename, **kwargs): kwargs.setdefault('driver','pytables') kwargs.setdefault('parameters',None) kwargs.setdefault('regions',[]) - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('units','seconds since 2018-01-01T00:00:00') @@ -658,10 +791,12 @@ def to_hdf(gdf, filename, **kwargs): # write pandas dataframe to pytables HDF5 def write_pytables(df, filename, attributes, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default keyword arguments kwargs.setdefault('parameters',None) kwargs.setdefault('regions',[]) - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) # write data to a pytables HDF5 file df.to_hdf(filename, 'sliderule_segments', format="table", mode="w") # add file attributes @@ -674,7 +809,7 @@ def write_pytables(df, filename, attributes, **kwargs): # set coordinate reference system as attribute fileID.root._v_attrs.crs = kwargs['crs'] # add geospatial attributes - if (kwargs['crs'] == 'EPSG:4326'): + if kwargs['crs'] in ('EPSG:4326', SLIDERULE_EPSG): fileID.root._v_attrs.geospatial_lat_units = \ attributes['geospatial_lat_units'] fileID.root._v_attrs.geospatial_lon_units = \ @@ -716,10 +851,12 @@ def write_pytables(df, filename, attributes, **kwargs): # write pandas dataframe to h5py HDF5 def write_h5py(df, filename, attributes, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default keyword arguments kwargs.setdefault('parameters',None) kwargs.setdefault('regions',[]) - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('units','seconds since 2018-01-01T00:00:00') # open HDF5 file object fileID = h5py.File(filename, mode='w') @@ -755,7 +892,7 @@ def write_h5py(df, filename, attributes, **kwargs): # set coordinate reference system as attribute fileID.attrs['crs'] = kwargs['crs'] # add geospatial attributes - if (kwargs['crs'] == 'EPSG:4326'): + if kwargs['crs'] in ('EPSG:4326', SLIDERULE_EPSG): fileID.attrs['geospatial_lat_units'] = \ attributes['geospatial_lat_units'] fileID.attrs['geospatial_lon_units'] = \ @@ -799,7 +936,7 @@ def write_h5py(df, filename, attributes, **kwargs): def from_hdf(filename, **kwargs): # set default keyword arguments kwargs.setdefault('driver','pytables') - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('return_parameters',False) @@ -815,8 +952,10 @@ def from_hdf(filename, **kwargs): # read pandas dataframe from pytables HDF5 def read_pytables(filename, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default crs - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('return_parameters',False) @@ -883,8 +1022,10 @@ def read_pytables(filename, **kwargs): # read pandas dataframe from h5py HDF5 def read_h5py(filename, **kwargs): + # add warning that function is deprecated + logger.critical(f"Deprecated. Will be removed in a future release") # set default crs - kwargs.setdefault('crs','EPSG:4326') + kwargs.setdefault('crs',SLIDERULE_EPSG) kwargs.setdefault('lon_key','longitude') kwargs.setdefault('lat_key','latitude') kwargs.setdefault('index_key','time') @@ -967,15 +1108,23 @@ def read_h5py(filename, **kwargs): return output # output formats wrapper -def to_file(gdf, filename, format='hdf', **kwargs): +def to_file(gdf, filename, format='parquet', **kwargs): if format.lower() in ('hdf','hdf5','h5'): to_hdf(gdf, filename, **kwargs) elif format.lower() in ('netcdf','nc'): to_nc(gdf, filename, **kwargs) + elif format.lower() in ('geojson','csv','shp'): + gdf.to_file(filename, **kwargs) + elif format.lower() in ('geoparquet','parquet'): + to_parquet(gdf, filename, **kwargs) # input formats wrapper -def from_file(filename, format='hdf', **kwargs): +def from_file(filename, format='parquet', **kwargs): if format.lower() in ('hdf','hdf5','h5'): return from_hdf(filename, **kwargs) - elif format.lower() in ('netcdf','nc'): + elif format.lower() in ('netcdf','netcdf4','nc'): return from_nc(filename, **kwargs) + elif format.lower() in ('geojson','csv','shp'): + return geopandas.from_file(filename, **kwargs) + elif format.lower() in ('geoparquet','parquet'): + return from_parquet(filename, **kwargs) diff --git a/clients/python/sliderule/ipysliderule.py b/clients/python/sliderule/ipysliderule.py index 739f36b37..54c3382d9 100644 --- a/clients/python/sliderule/ipysliderule.py +++ b/clients/python/sliderule/ipysliderule.py @@ -520,8 +520,7 @@ def __init__(self, **kwargs): # dropdown menu for selecting variable to draw on map variable_list = ['h_mean', 'h_sigma', 'dh_fit_dx', 'dh_fit_dy', - 'rms_misfit', 'w_surface_window_final', 'delta_time', - 'cycle', 'rgt'] + 'rms_misfit', 'w_surface_window_final', 'cycle', 'rgt'] self.variable = ipywidgets.Dropdown( options=variable_list, value='h_mean', @@ -659,6 +658,16 @@ def __init__(self, **kwargs): # watch plot kind widgets for changes self.plot_kind.observe(self.set_plot_kind) + # selection for output file format + format_options = ["GeoJSON","csv","geoparquet","netCDF","HDF5"] + self.file_format = ipywidgets.Dropdown( + options=format_options, + value="geoparquet", + description="Format:", + description_tooltip="Format: Output file format", + disabled=False + ) + # button and label for output file selection self.file = copy.copy(self.atl06_filename) self.savebutton = ipywidgets.Button( @@ -1001,12 +1010,16 @@ def saveas_file(self, b): root = tkinter.Tk() root.withdraw() root.call('wm', 'attributes', '.', '-topmost', True) - filetypes = (("HDF5 file", "*.h5"), + filetypes = ( + ("geoparquet file", "*.parquet"), + ("GeoJSON file", "*.geojson"), + ("csv file", "*.csv"), ("netCDF file", "*.nc"), + ("HDF5 file", "*.h5"), ("All Files", "*.*")) b.files = tkinter.filedialog.asksaveasfilename( initialfile=self.file, - defaultextension='h5', + defaultextension='parquet', filetypes=filetypes) self.savelabel.value = b.files self.file = b.files @@ -1024,8 +1037,12 @@ def select_file(self, b): root = tkinter.Tk() root.withdraw() root.call('wm', 'attributes', '.', '-topmost', True) - filetypes = (("HDF5 file", "*.h5"), + filetypes = ( + ("geoparquet file", "*.parquet"), + ("GeoJSON file", "*.geojson"), + ("csv file", "*.csv"), ("netCDF file", "*.nc"), + ("HDF5 file", "*.h5"), ("All Files", "*.*")) b.files = tkinter.filedialog.askopenfilename( defaultextension='h5', @@ -1041,54 +1058,89 @@ def set_loadfile(self, sender): self.file = self.loadlabel.value @property - def atl03_filename(self): + def atl03_filename(self) -> str: """default input and output file string for ATL03 requests """ # get sliderule submission time now = datetime.datetime.now().strftime('%Y%m%d%H%M%S') - args = (now, self.release.value) - return "ATL03-SR_{0}_{1}.h5".format(*args) + args = (now, self.release.value, self.suffix) + return "ATL03-SR_{0}_{1}.{2}".format(*args) @property - def atl06_filename(self): + def atl06_filename(self) -> str: """default input and output file string for ATL06 requests """ # get sliderule submission time now = datetime.datetime.now().strftime('%Y%m%d%H%M%S') - args = (now, self.release.value) - return "ATL06-SR_{0}_{1}.h5".format(*args) + args = (now, self.release.value, self.suffix) + return "ATL06-SR_{0}_{1}.{2}".format(*args) @property - def atl08_filename(self): + def atl08_filename(self) -> str: """default input and output file string for ATL08 requests """ # get sliderule submission time now = datetime.datetime.now().strftime('%Y%m%d%H%M%S') - args = (now, self.release.value) - return "ATL08-SR_{0}_{1}.h5".format(*args) + args = (now, self.release.value, self.suffix) + return "ATL08-SR_{0}_{1}.{2}".format(*args) @property - def format(self): + def format(self) -> str: """return the file format from file string """ - hdf = ('h5','hdf5','hdf') - netcdf = ('nc','netcdf','nc3') - if self.file.endswith(hdf): - return 'hdf' - elif self.file.endswith(netcdf): - return 'netcdf' - else: - return '' + types = {} + types['hdf'] = ('h5','hdf5','hdf') + types['netcdf'] = ('nc','netcdf','nc3') + types['csv'] = ('csv','txt') + types['geojson'] = ('geojson','json') + types['parquet'] = ('geoparquet','parquet','pq') + for key, val in types.items(): + if self.file.endswith(val): + return key + # return empty string if no match + return '' + + @property + def suffix(self) -> str: + """return the file suffix + """ + SUFFIX = { + "GeoJSON":"geojson", + "csv":"csv", + "geoparquet":"parquet", + "netCDF":"nc", + "HDF5":"h5" + } + try: + return SUFFIX[self.file_format.value] + except KeyError: + return "" @property - def _r(self): + def mime_type(self) -> str: + """return the internet media type of the file format + """ + MIME = { + "GeoJSON":"text/json", + "csv":"text/csv", + "geoparquet":"application/vnd.apache.parquet", + "netCDF":"application/x-netcdf", + "HDF5":"application/x-hdf5" + } + try: + return MIME[self.file_format.value] + except KeyError: + return "" + + @property + def _r(self) -> str: """return string for reversed Matplotlib colormaps """ cmap_reverse_flag = '_r' if self.reverse.value else '' return cmap_reverse_flag @property - def colormap(self): + def colormap(self) -> str: """return string for Matplotlib colormaps """ return self.cmap.value + self._r @@ -1885,12 +1937,16 @@ def _load_dict(data): # create traitlets of basemap providers basemaps = _load_dict(providers) +# set default map dimensions +_default_layout = ipywidgets.Layout(width='100%', height='600px') +_default_max_area = 500e3 # draw ipyleaflet map class leaflet: def __init__(self, projection, **kwargs): # set default keyword arguments kwargs.setdefault('map',None) + kwargs.setdefault('layout', _default_layout) kwargs.setdefault('prefer_canvas', False) kwargs.setdefault('attribution', False) kwargs.setdefault('zoom_control', False) @@ -1901,13 +1957,14 @@ def __init__(self, projection, **kwargs): kwargs.setdefault('color', 'green') # create basemap in projection if (projection == 'Global'): - kwargs.setdefault('center', (39,-108)) - kwargs.setdefault('zoom', 9) + kwargs.setdefault('center', (40,-100)) + kwargs.setdefault('zoom', 4) self.map = ipyleaflet.Map(center=kwargs['center'], zoom=kwargs['zoom'], max_zoom=15, world_copy_jump=True, prefer_canvas=kwargs['prefer_canvas'], attribution_control=kwargs['attribution'], - basemap=ipyleaflet.basemaps.Esri.WorldTopoMap) + basemap=ipyleaflet.basemaps.Esri.WorldTopoMap, + layout=kwargs['layout']) self.crs = 'EPSG:3857' elif (projection == 'North'): kwargs.setdefault('center', (90,0)) @@ -1917,20 +1974,22 @@ def __init__(self, projection, **kwargs): prefer_canvas=kwargs['prefer_canvas'], attribution_control=kwargs['attribution'], basemap=basemaps.Esri.ArcticOceanBase, - crs=projections.EPSG5936.ESRIBasemap) + crs=projections.EPSG5936.ESRIBasemap, + layout=kwargs['layout']) # add arctic ocean reference basemap reference = basemaps.Esri.ArcticOceanReference self.map.add(ipyleaflet.basemap_to_tiles(reference)) self.crs = 'EPSG:5936' elif (projection == 'South'): kwargs.setdefault('center', (-90,0)) - kwargs.setdefault('zoom', 2) + kwargs.setdefault('zoom', 3) self.map = ipyleaflet.Map(center=kwargs['center'], zoom=kwargs['zoom'], max_zoom=9, prefer_canvas=kwargs['prefer_canvas'], attribution_control=kwargs['attribution'], basemap=basemaps.Esri.AntarcticBasemap, - crs=projections.EPSG3031.ESRIBasemap) + crs=projections.EPSG3031.ESRIBasemap, + layout=kwargs['layout']) self.crs = 'EPSG:3031' else: # use a predefined ipyleaflet map @@ -1973,7 +2032,7 @@ def __init__(self, projection, **kwargs): draw_control.rectangle = dict(shapeOptions=shapeOptions, metric=['km','m']) draw_control.polygon = dict(shapeOptions=shapeOptions, - allowIntersection=False,showArea=True,metric=['km','m']) + allowIntersection=False, showArea=True, metric=['km','m']) # create regions self.regions = [] draw_control.on_draw(self.handle_draw) @@ -2147,10 +2206,12 @@ def handle_interaction(self, **kwargs): Longitude: {d[1]:8.4f}\u00B0""".format(d=[lat,lon]) # keep track of rectangles and polygons drawn on map - def handle_draw(self, obj, action, geo_json): + def handle_draw(self, obj, action, geo_json, **kwargs): """callback for handling draw events and interactively creating SlideRule region objects """ + kwargs.setdefault('check_valid', True) + kwargs.setdefault('maximum_area', _default_max_area) lon,lat = np.transpose(geo_json['geometry']['coordinates']) lon = sliderule.io.wrap_longitudes(lon) cx,cy = sliderule.io.centroid(lon,lat) @@ -2159,6 +2220,11 @@ def handle_draw(self, obj, action, geo_json): if (wind > 0): lon = lon[::-1] lat = lat[::-1] + # calculate area of region + area = sliderule.io.area(lon,lat) + if kwargs['check_valid'] and (area > kwargs['maximum_area']): + logger.warning(f"Region is too large: {area:0.0f} km^2") + return # create sliderule region from list region = sliderule.io.to_region(lon,lat) # append coordinates to list @@ -2210,7 +2276,7 @@ def GeoData(self, gdf, **kwargs): kwargs.setdefault('stride', None) kwargs.setdefault('max_plot_points', 10000) kwargs.setdefault('tooltip', True) - kwargs.setdefault('fields', self.default_atl06_fields()) + kwargs.setdefault('fields', []) kwargs.setdefault('colorbar', True) kwargs.setdefault('position', 'topright') # add warning that function is deprecated @@ -2377,41 +2443,6 @@ def add_colorbar(self, **kwargs): self.map.add(self.colorbar) plt.close() - @staticmethod - def default_atl03_fields(): - """List of ATL03 tooltip fields - """ - return ['atl03_cnf', 'atl08_class', 'cycle', 'height', - 'pair', 'rgt', 'segment_id', 'track', 'yapc_score'] - - @staticmethod - def default_atl06_fields(): - """List of ATL06-SR tooltip fields - """ - return ['cycle', 'dh_fit_dx', 'gt', 'h_mean', - 'h_sigma', 'rgt', 'rms_misfit', 'w_surface_window_final'] - - @staticmethod - def default_atl08_fields(): - """List of ATL08-SR tooltip fields - """ - return ['canopy_openness', 'cycle', 'gt', 'h_canopy', - 'h_min_canopy', 'h_mean_canopy', - 'h_max_canopy', 'h_te_median', 'rgt'] - - @staticmethod - def default_mosaic_fields(**kwargs): - kwargs.setdefault('with_flags', False) - kwargs.setdefault('zonal_stats', False) - """List of mosaic tooltip fields - """ - columns = ['time','value'] - if kwargs['with_flags']: - columns += ['flags'] - if kwargs['zonal_stats']: - columns += ['count','min','max','mean','median','stdev','mad'] - return [f'mosaic.{c}' for c in columns] - @gpd.pd.api.extensions.register_dataframe_accessor("leaflet") class LeafletMap: """A geopandas GeoDataFrame extension for interactive map plotting, @@ -2425,12 +2456,12 @@ def __init__(self, gdf): # initialize geodataframe self._gdf = gdf # initialize data and colorbars - self.geojson = None - self.tooltip = None - self.tooltip_width = None - self.tooltip_height = None - self.fields = [] - self.colorbar = None + self._geojson = None + self._tooltip = None + self._tooltip_width = None + self._tooltip_height = None + self._fields = [] + self._colorbar = None # initialize hover control self.hover_control = None # initialize selected feature @@ -2478,8 +2509,8 @@ def GeoData(self, m, **kwargs): self.map = m self.crs = m.crs['name'] # remove any prior instances of a data layer - if self.geojson is not None: - self.map.remove(self.geojson) + if self._geojson is not None: + self.map.remove(self._geojson) if kwargs['stride'] is not None: stride = np.copy(kwargs['stride']) elif (self._gdf.shape[0] > kwargs['max_plot_points']): @@ -2487,62 +2518,65 @@ def GeoData(self, m, **kwargs): else: stride = 1 # sliced geodataframe for plotting - geodataframe = self._gdf[slice(None,None,stride)] + self._gdf_selected = self._gdf[slice(None,None,stride)] self.column_name = copy.copy(kwargs['column_name']) - geodataframe['data'] = geodataframe[self.column_name] - # set colorbar limits to 2-98 percentile - # if not using a defined plot range - clim = geodataframe['data'].quantile((0.02, 0.98)).values - if kwargs['vmin'] is None: - vmin = clim[0] - else: - vmin = np.copy(kwargs['vmin']) - if kwargs['vmax'] is None: - vmax = clim[-1] - else: - vmax = np.copy(kwargs['vmax']) + self._gdf_selected['data'] = self._gdf_selected[self.column_name] + # get the normalization bounds + self.get_norm_bounds(**kwargs) # create matplotlib normalization if kwargs['norm'] is None: - norm = colors.Normalize(vmin=vmin, vmax=vmax, clip=True) + self.norm = colors.Normalize(vmin=self.vmin, vmax=self.vmax, clip=True) else: - norm = copy.copy(kwargs['norm']) + self.norm = copy.copy(kwargs['norm']) # normalize data to be within vmin and vmax - normalized = norm(geodataframe['data']) + normalized = self.norm(self._gdf_selected['data']) + # get colormap + self.cmap = copy.copy(cm.get_cmap(kwargs['cmap'])) # create HEX colors for each point in the dataframe - geodataframe["color"] = np.apply_along_axis(colors.to_hex, 1, - cm.get_cmap(kwargs['cmap'], 256)(normalized)) + self._gdf_selected["color"] = np.apply_along_axis(colors.to_hex, 1, + cm.get_cmap(self.cmap.name, 256)(normalized)) # leaflet map point style - point_style = {key:kwargs[key] for key in ['radius','fillOpacity','weight']} + self._point_style = { + key:kwargs[key] for key in ['radius','fillOpacity','weight'] + } # convert to GeoJSON object - self.geojson = ipyleaflet.GeoJSON(data=geodataframe.__geo_interface__, - point_style=point_style, style_callback=self.style_callback) + self._geojson = ipyleaflet.GeoJSON( + data=self._gdf_selected.__geo_interface__, + point_style=self._point_style, + style_callback=self.style_callback + ) # add GeoJSON object to map - self.map.add(self.geojson) + self.map.add(self._geojson) # fields for tooltip views if kwargs['fields'] is None: - self.fields = geodataframe.columns.drop( - [geodataframe.geometry.name, "data", "color"]) + self._fields = self._gdf_selected.columns.drop( + [self._gdf_selected.geometry.name, "data", "color"]) else: - self.fields = copy.copy(kwargs['fields']) + self._fields = copy.copy(kwargs['fields']) # add hover tooltips if kwargs['tooltip']: - self.tooltip = ipywidgets.HTML() - self.tooltip.layout.margin = "0px 20px 20px 20px" - self.tooltip.layout.visibility = 'hidden' - self.tooltip_height = kwargs['tooltip_height'] - self.tooltip_width = kwargs['tooltip_width'] + self._tooltip = ipywidgets.HTML() + self._tooltip.layout.margin = "0px 20px 20px 20px" + self._tooltip.layout.visibility = 'hidden' + self._tooltip_height = kwargs['tooltip_height'] + self._tooltip_width = kwargs['tooltip_width'] # create widget for hover tooltips self.hover_control = ipyleaflet.WidgetControl( - widget=self.tooltip, + widget=self._tooltip, position='bottomright') - self.geojson.on_hover(self.handle_hover) - self.geojson.on_msg(self.handle_mouseout) - self.geojson.on_click(self.handle_click) + self._geojson.on_hover(self.handle_hover) + self._geojson.on_msg(self.handle_mouseout) + self._geojson.on_click(self.handle_click) # add colorbar - if kwargs['colorbar']: - self.add_colorbar(column_name=self.column_name, - cmap=kwargs['cmap'], norm=norm, - position=kwargs['position']) + self.colorbar = kwargs['colorbar'] + self.colorbar_position = kwargs['position'] + if self.colorbar: + self.add_colorbar( + column_name=self.column_name, + cmap=self.cmap, + norm=self.norm, + position=self.colorbar_position + ) # functional call for setting colors of each point def style_callback(self, feature): @@ -2558,12 +2592,12 @@ def handle_hover(self, feature, **kwargs): """callback for creating hover tooltips """ # combine html strings for hover tooltip - self.tooltip.value = '{0}: {1}
'.format('id',feature['id']) - self.tooltip.value += '
'.join(['{0}: {1}'.format(field, - feature["properties"][field]) for field in self.fields]) - self.tooltip.layout.width = self.tooltip_width - self.tooltip.layout.height = self.tooltip_height - self.tooltip.layout.visibility = 'visible' + self._tooltip.value = '{0}: {1}
'.format('id',feature['id']) + self._tooltip.value += '
'.join(['{0}: {1}'.format(field, + feature["properties"][field]) for field in self._fields]) + self._tooltip.layout.width = self._tooltip_width + self._tooltip.layout.height = self._tooltip_height + self._tooltip.layout.visibility = 'visible' self.map.add(self.hover_control) def handle_mouseout(self, _, content, buffers): @@ -2571,10 +2605,10 @@ def handle_mouseout(self, _, content, buffers): """ event_type = content.get('type', '') if event_type == 'mouseout': - self.tooltip.value = '' - self.tooltip.layout.width = "0px" - self.tooltip.layout.height = "0px" - self.tooltip.layout.visibility = 'hidden' + self._tooltip.value = '' + self._tooltip.layout.width = "0px" + self._tooltip.layout.height = "0px" + self._tooltip.layout.visibility = 'hidden' self.map.remove(self.hover_control) # functional calls for click events @@ -2593,13 +2627,13 @@ def handle_region(self, action, **kwargs): """callback for handling region deletions """ # remove any prior instances of a data layer - if (action == 'deleted') and self.geojson is not None: - self.remove(self.geojson) - self.geojson = None + if (action == 'deleted') and self._geojson is not None: + self.remove(self._geojson) + self._geojson = None # remove any prior instances of a colorbar - if (action == 'deleted') and self.colorbar is not None: - self.remove(self.colorbar) - self.colorbar = None + if (action == 'deleted') and self._colorbar is not None: + self.remove(self._colorbar) + self._colorbar = None # remove map layers def remove(self, layer): @@ -2619,6 +2653,121 @@ def remove(self, layer): logger.error(traceback.format_exc()) pass + def get_norm_bounds(self, **kwargs): + # set default keyword arguments + kwargs.setdefault('vmin', None) + kwargs.setdefault('vmax', None) + # set colorbar limits to 2-98 percentile + # if not using a defined plot range + clim = self._gdf_selected['data'].quantile((0.02, 0.98)).values + # set minimum for normalization + fmin = np.finfo(np.float64).min + if (kwargs['vmin'] is None) or np.isclose(kwargs['vmin'], fmin): + self.vmin = clim[0] + self._dynamic = True + else: + self.vmin = np.copy(kwargs['vmin']) + self._dynamic = False + # set maximum for normalization + fmax = np.finfo(np.float64).max + if (kwargs['vmax'] is None) or np.isclose(kwargs['vmax'], fmax): + self.vmax = clim[-1] + self._dynamic = True + else: + self.vmax = np.copy(kwargs['vmax']) + self._dynamic = False + + def redraw(self, *args, **kwargs): + """ + Redraw the GeoJSON on the map + """ + # normalize data to be within vmin and vmax + normalized = self.norm(self._gdf_selected['data']) + # create HEX colors for each point in the dataframe + self._gdf_selected["color"] = np.apply_along_axis(colors.to_hex, 1, + cm.get_cmap(self.cmap.name, 256)(normalized)) + # update data within GeoJSON object + self._geojson.data = self._gdf_selected.__geo_interface__ + + def redraw_colorbar(self, *args, **kwargs): + """ + Redraw the colorbar on the map + """ + try: + if self.colorbar: + self.add_colorbar( + column_name=self.column_name, + cmap=self.cmap, + norm=self.norm, + position=self.colorbar_position + ) + except Exception as exc: + pass + + # observe changes in widget parameters + def set_observables(self, widget, **kwargs): + """observe changes in widget parameters + """ + # set default keyword arguments + # to map widget changes to functions + kwargs.setdefault('variable', [self.set_column_name]) + kwargs.setdefault('cmap', [self.set_colormap]) + kwargs.setdefault('reverse', [self.set_colormap]) + # connect each widget with a set function + for key, val in kwargs.items(): + # try to retrieve the functional + try: + observable = getattr(widget, key) + except AttributeError as exc: + continue + # assert that observable is an ipywidgets object + assert isinstance(observable, ipywidgets.widgets.widget.Widget) + assert hasattr(observable, 'observe') + # for each functional to map + for i, functional in enumerate(val): + # try to connect the widget to the functional + try: + observable.observe(functional) + except (AttributeError, NameError, ValueError) as exc: + pass + + def set_column_name(self, sender): + """update the dataframe variable for a new selected column + """ + # only update variable if a new final + if isinstance(sender['new'], str): + self.column_name = sender['new'] + else: + return + # reduce to variable + self._gdf_selected['data'] = self._gdf_selected[self.column_name] + # check if dynamic normalization is enabled + if self._dynamic: + self.get_norm_bounds() + self.norm.vmin = self.vmin + self.norm.vmax = self.vmax + # try to redraw the selected dataset + self.redraw() + self.redraw_colorbar() + + def set_colormap(self, sender): + """update the colormap for the selected variable + """ + # only update colormap if a new final + if isinstance(sender['new'], str): + cmap_name = self.cmap.name + cmap_reverse_flag = '_r' if cmap_name.endswith('_r') else '' + self.cmap = cm.get_cmap(sender['new'] + cmap_reverse_flag) + elif isinstance(sender['new'], bool): + cmap_name = self.cmap.name.strip('_r') + cmap_reverse_flag = '_r' if sender['new'] else '' + self.cmap = cm.get_cmap(cmap_name + cmap_reverse_flag) + else: + return + # try to redraw the selected dataset + self.redraw() + self.redraw_colorbar() + # add colorbar widget to leaflet map def add_colorbar(self, **kwargs): """Creates colorbars on leaflet maps @@ -2638,21 +2787,22 @@ def add_colorbar(self, **kwargs): kwargs.setdefault('cmap', 'viridis') kwargs.setdefault('norm', None) kwargs.setdefault('alpha', 1.0) - kwargs.setdefault('orientation', 'horizontal') + kwargs.setdefault('orientation', 'vertical') kwargs.setdefault('position', 'topright') - kwargs.setdefault('width', 6.0) - kwargs.setdefault('height', 0.4) + kwargs.setdefault('width', 0.2) + kwargs.setdefault('height', 3.0) # remove any prior instances of a colorbar - if self.colorbar is not None: - self.map.remove(self.colorbar) - # colormap for colorbar - cmap = cm.get_cmap(kwargs['cmap']) + if self._colorbar is not None: + self.map.remove(self._colorbar) # create matplotlib colorbar _, ax = plt.subplots(figsize=(kwargs['width'], kwargs['height'])) - cbar = matplotlib.colorbar.ColorbarBase(ax, cmap=cmap, - norm=kwargs['norm'], alpha=kwargs['alpha'], + cbar = matplotlib.colorbar.ColorbarBase(ax, + cmap=kwargs['cmap'], + norm=kwargs['norm'], + alpha=kwargs['alpha'], orientation=kwargs['orientation'], - label=kwargs['column_name']) + label=kwargs['column_name'] + ) cbar.solids.set_rasterized(True) cbar.ax.tick_params(which='both', width=1, direction='in') # save colorbar to in-memory png object @@ -2661,12 +2811,50 @@ def add_colorbar(self, **kwargs): png.seek(0) # create output widget output = ipywidgets.Image(value=png.getvalue(), format='png') - self.colorbar = ipyleaflet.WidgetControl(widget=output, - transparent_bg=True, position=kwargs['position']) + self._colorbar = ipyleaflet.WidgetControl( + widget=output, + transparent_bg=True, + position=kwargs['position'] + ) # add colorbar - self.map.add(self.colorbar) + self.map.add(self._colorbar) plt.close() + @staticmethod + def default_atl03_fields(): + """List of ATL03 tooltip fields + """ + return ['atl03_cnf', 'atl08_class', 'cycle', 'height', + 'pair', 'rgt', 'segment_id', 'track', 'yapc_score'] + + @staticmethod + def default_atl06_fields(): + """List of ATL06-SR tooltip fields + """ + return ['cycle', 'dh_fit_dx', 'gt', 'h_mean', + 'h_sigma', 'rgt', 'rms_misfit', 'w_surface_window_final'] + + @staticmethod + def default_atl08_fields(): + """List of ATL08-SR tooltip fields + """ + return ['canopy_openness', 'cycle', 'gt', 'h_canopy', + 'h_min_canopy', 'h_mean_canopy', + 'h_max_canopy', 'h_te_median', 'rgt'] + + @staticmethod + def default_mosaic_fields(**kwargs): + kwargs.setdefault('with_flags', False) + kwargs.setdefault('zonal_stats', False) + """List of mosaic tooltip fields + """ + columns = ['time','value'] + if kwargs['with_flags']: + columns += ['flags'] + if kwargs['zonal_stats']: + columns += ['count','min','max','mean','median','stdev','mad'] + return [f'mosaic.{c}' for c in columns] + @gpd.pd.api.extensions.register_dataframe_accessor("icesat2") class ICESat2: """A geopandas GeoDataFrame extension for plotting ICESat-2 transects diff --git a/docs/rtd/source/user_guide/ICESat-2.rst b/docs/rtd/source/user_guide/ICESat-2.rst index 3183cd917..eafb1f9ff 100644 --- a/docs/rtd/source/user_guide/ICESat-2.rst +++ b/docs/rtd/source/user_guide/ICESat-2.rst @@ -135,22 +135,22 @@ For example: } 2.5.1 ATL03 Subsetted Ancillary Data -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +#################################### Ancillary data returned from the ``"atl03s"`` and ``"atl03sp"`` APIs are per-photon values that are read from the ATL03 granules. No processing is performed on the data read out of the ATL03 granule. The fields must come from either a per-photon variable (atl03_ph_fields), or a per-segment variable (atl03_geo_fields). 2.5.2 ATL06-SR Ancillary Data -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +############################# Ancillary data returned from the ``"atl06"`` and ``"atl06p"`` APIs come from the ancillary fields specified for ATL03, but instead of being returned as-is, they are processed using the ATL06 least-squares-fit algorithm, and only the result is returned. In other words, ancillary data points from ATL03 to be included in an ATL06-SR result are treated just like the h_mean, latitude, and longitude variables, and returned as a fitted double-precision floating point value. 2.5.3 ATL06 Subsetted Ancillary Data -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +#################################### Ancillary data returned from the ``"atl06s"`` and ``"atl06sp"`` APIs come from the land_ice_segments group of the ATL06 granules. The data is mostly returned as-is, with one exception. Double-precision and single-precision floating point variables are checked to see if they contain the maximum value of their respective encodings, and if so, a floating point NaN (not-a-number) is returned instead. This check is not performed for integer variables because the maximum value of an encoded integer can sometimes be a valid value (e.g. bit masks). 2.5.4 ATL08-PhoREAL Ancillary Data -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +################################## Ancillary data returned from the ``"atl08"`` and ``"atl08p"`` APIs come from the land_segments group of the ATL08 granules. The data goes through a series of processing steps before being returned back to the user as per-extent (i.e. variable-length segment) result values. @@ -390,7 +390,7 @@ The vegetation GeoDataFrame has the following columns: 4. Callbacks -============= +============ For large processing requests, it is possible that the data returned from the API is too large or impractical to fit in the local memory of the Python interpreter making the request. In these cases, certain APIs in the SlideRule Python client allow the calling application to provide a callback function that is called for every result that is returned by the servers. If a callback is supplied, the API will not return back to the calling application anything associated with the supplied record types, but assumes the callback fully handles processing the data. @@ -420,7 +420,7 @@ Here is an example of a callback being used for the ``atl03sp`` function: 5. Endpoints -============= +============ atl06 -----