diff --git a/pointCollection/ATL06/data.py b/pointCollection/ATL06/data.py index 5932a58..d0562bb 100644 --- a/pointCollection/ATL06/data.py +++ b/pointCollection/ATL06/data.py @@ -33,6 +33,7 @@ def __init__(self, fields=None, SRS_proj4=None, field_dict=None, beam_pair=2, co self.fields=list(fields) self.fields=fields self.SRS_proj4=SRS_proj4 + self.coordinates=['longitude','latitude'] self.columns=columns self.shape=(0,2) self.size=0 diff --git a/pointCollection/ATL11_prerelease/__init__.py b/pointCollection/ATL11_prerelease/__init__.py deleted file mode 100644 index 22a3805..0000000 --- a/pointCollection/ATL11_prerelease/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .data import data -from .crossover_data import crossover_data - diff --git a/pointCollection/ATL11_prerelease/crossover_data.py b/pointCollection/ATL11_prerelease/crossover_data.py deleted file mode 100755 index 1a7e517..0000000 --- a/pointCollection/ATL11_prerelease/crossover_data.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed May 20 13:27:26 2020 - -@author: ben -""" - -import pointCollection as pc -import numpy as np -import h5py - -class crossover_data(pc.data): - def __init__(self, pair=2, **kwargs): - self.pair=pair - self.pair_name = f'pt{int(self.pair)}' - self.columns=2 - super().__init__(**kwargs) - - - def __default_XO_field_dict__(self): - - field_dict= {'crossing_track_data':\ - ['along_track_rss', 'atl06_quality_summary', 'cycle_number', - 'delta_time', 'h_corr', 'h_corr_sigma', 'h_corr_sigma_systematic', - 'latitude', 'longitude', 'ref_pt', 'rgt']} - return {self.pair_name+'/'+key:field_dict[key] for key in field_dict} - - - def from_h5(self, filename, pair=None): - if pair is None: - pair=self.pair - else: - self.pair_name = f'pt{int(self.pair)}' - D_at=pc.ATL11.data().from_h5(filename, pair=pair) - D_xo=pc.data().from_h5(filename, group=self.pair_name+'/crossing_track_data', field_dict=self.__default_XO_field_dict__()) - D_xo.index(np.isfinite(D_xo.ref_pt) & np.isfinite(D_xo.cycle_number)) - with h5py.File(filename,'r') as h5f: - rgt=int(h5f[self.pair_name].attrs['ReferenceGroundTrack']) - - D_at.assign({'rgt':np.zeros_like(D_at.delta_time)+rgt}) - n_cycles=D_at.cycle_number[-1, -1] - - u_pt_xo=np.unique(D_xo.ref_pt) - D_at.index(np.in1d(D_at.ref_pt[:,0], u_pt_xo)) - - theshape=(u_pt_xo.size, n_cycles,2) - self.assign({field:np.zeros(theshape)+np.nan for field in D_xo.fields}) - - row=np.searchsorted( u_pt_xo, D_at.ref_pt.ravel()) - col=D_at.cycle_number.ravel().astype(int)-1 - for field in self.fields: - if field not in D_at.fields: - continue - getattr(self, field).flat[np.ravel_multi_index((row, col, np.zeros_like(row, dtype=int)), theshape)]=getattr(D_at, field).ravel() - - row=np.searchsorted(u_pt_xo, D_xo.ref_pt) - col=D_xo.cycle_number.astype(int)-1 - for field in self.fields: - getattr(self, field).flat[np.ravel_multi_index((row, col, 1+np.zeros_like(row, dtype=int)), theshape)]=getattr(D_xo, field) - - self.__update_size_and_shape__() - return self - diff --git a/pointCollection/ATL11_prerelease/data.py b/pointCollection/ATL11_prerelease/data.py deleted file mode 100755 index 36e7675..0000000 --- a/pointCollection/ATL11_prerelease/data.py +++ /dev/null @@ -1,84 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed May 6 12:07:54 2020 - -@author: ben -""" - -import numpy as np -import pointCollection as pc -import os -import h5py - - -class data(pc.data): - def __init__(self, pair=2, **kwargs): - self.pair=pair - self.pair_name = f'pt{int(self.pair)}' - self.cycle_number=np.array([]) - super().__init__(**kwargs) - - - def __default_field_dict__(self, field_weight='light'): - if field_weight=='light': - field_dict={'corrected_h':['latitude','longitude','h_corr',\ - 'h_corr_sigma', 'h_corr_sigma_systematic',\ - 'delta_time','quality_summary', 'ref_pt'], \ - 'ref_surf': ['dem_h', 'x_atc']} - return self.__convert_field_dict__(field_dict) - - def __convert_field_dict__(self, field_dict): - - temp={} - for key in field_dict: - if key is None: - temp[self.pair_name]=field_dict[key] - if key in ['__calc_internal__']: - # skip appending the pair name to the __calc_internal__ field - temp[key]=field_dict[key] - else: - temp[self.pair_name+'/'+key]=field_dict[key] - return temp - - - def __tile_fields__(self): - rows=self.shape[0] - cols=self.shape[1] - for field in self.fields: - temp=getattr(self, field) - if temp.size==cols: - if temp.ndim==1: - temp.shape=(1,cols) - setattr(self, field, np.tile(temp, [rows,1])) - elif temp.size==rows: - if temp.ndim==1: - temp.shape=(rows,1) - setattr(self, field, np.tile(temp, [1, cols])) - self.__update_size_and_shape__() - return self - - def __internal_field_calc__(self, field_dict): - if 'rgt' in field_dict['__calc_internal__']: - self.__update_size_and_shape__() - with h5py.File(self.filename,'r') as h5f: - self.rgt=h5f['/ancillary_data/start_rgt'][0]+np.zeros(self.shape) - - def from_h5(self, filename, pair=None, field_weight='light', tile_fields=True, **kwargs): - if pair is not None: - self.pair_name=f'pt{int(pair)}' - self.field_dict=self.__default_field_dict__(field_weight=field_weight) - with h5py.File(filename,'r') as h5f: - cycle_number = np.array(h5f[self.pair_name]['corrected_h']['cycle_number']) - - if 'field_dict' in kwargs: - kwargs['field_dict']=self.__convert_field_dict__(kwargs['field_dict'].copy()) - - super().from_h5(filename, **kwargs) - self.cycle_number=cycle_number - self.columns=len(self.cycle_number) - self.shape=(self.latitude.size, self.columns) - if tile_fields: - self.fields += ['cycle_number'] - self.__tile_fields__() - return self diff --git a/pointCollection/ATL11_prerelease/demo.py b/pointCollection/ATL11_prerelease/demo.py deleted file mode 100755 index ab901b1..0000000 --- a/pointCollection/ATL11_prerelease/demo.py +++ /dev/null @@ -1,12 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed May 6 14:10:10 2020 - -@author: ben -""" - -import pointCollection as pc - -SRS_proj4='+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs' -gI=pc.geoIndex(SRS_proj4=SRS_proj4).for_file('/Volumes/ice2/ben/scf/AA_11/U07/ATL11_000110_0306_02_vU07.h5', file_type='ATL11') diff --git a/pointCollection/CS2/data.py b/pointCollection/CS2/data.py index 0165abc..878417a 100755 --- a/pointCollection/CS2/data.py +++ b/pointCollection/CS2/data.py @@ -35,7 +35,6 @@ import numpy as np import pointCollection as pc -import netCDF4 import re import os @@ -918,6 +917,7 @@ def cryosat_baseline_D(self, full_filename, unpack=False): """ Read L2 MDS variables for CryoSat Baseline D (netCDF4) """ + import netCDF4 # open netCDF4 file for reading fid = netCDF4.Dataset(os.path.expanduser(full_filename),'r') # use original unscaled units unless unpack=True diff --git a/pointCollection/CS2_wfm/data.py b/pointCollection/CS2_wfm/data.py index 4d8dc99..468aca9 100755 --- a/pointCollection/CS2_wfm/data.py +++ b/pointCollection/CS2_wfm/data.py @@ -32,7 +32,6 @@ import numpy as np import pointCollection as pc -import netCDF4 import re import os @@ -1461,6 +1460,7 @@ def cryosat_baseline_D(self, full_filename, unpack=False): """ Read L1b MDS variables for CryoSat Baseline D (netCDF4) """ + import netCDF4 # open netCDF4 file for reading fid = netCDF4.Dataset(os.path.expanduser(full_filename),'r') # use original unscaled units unless unpack=True diff --git a/pointCollection/__init__.py b/pointCollection/__init__.py index cb3bf21..e99d7ce 100644 --- a/pointCollection/__init__.py +++ b/pointCollection/__init__.py @@ -10,21 +10,21 @@ from .tile import tile from .tilingSchema import tilingSchema from .geoIndex import geoIndex -# data classes -from pointCollection import ATM_Qfit -from pointCollection import ATM_ICESSN -from pointCollection import CS2_wfm -from pointCollection import CS2 -from pointCollection import CS2_retracked_POCA -from pointCollection import CS2_retracked_SW -from pointCollection import ATL06 -from pointCollection import indexedH5 -from pointCollection import ATM_WF -from pointCollection import glah12 -from pointCollection import glah06 -from pointCollection import ATL11 -from pointCollection import ATL11_prerelease from pointCollection import grid +# data classes — loaded lazily on first access (e.g. pc.ATL06, pc.ATM_Qfit) +_SUBPACKAGES = frozenset({ + 'ATM_Qfit', 'ATM_ICESSN', 'CS2_wfm', 'CS2', 'CS2_retracked_POCA', + 'CS2_retracked_SW', 'ATL06', 'indexedH5', 'ATM_WF', 'glah12', + 'glah06', 'ATL11', 'grid', +}) + +def __getattr__(name): + if name in _SUBPACKAGES: + import importlib + mod = importlib.import_module(f'pointCollection.{name}') + globals()[name] = mod + return mod + raise AttributeError(f"module 'pointCollection' has no attribute {name!r}") # utilities from .tools.bin_rows import bin_rows @@ -44,4 +44,3 @@ from .xover_search import cross_tracks from pointCollection.ps_scale_for_lat import ps_scale_for_lat from pointCollection.reconstruct_ATL06_tracks import reconstruct_ATL06_tracks -#from pointCollection.dataPicker import dataPicker diff --git a/pointCollection/data.py b/pointCollection/data.py index 9676807..eca29cc 100644 --- a/pointCollection/data.py +++ b/pointCollection/data.py @@ -4,10 +4,10 @@ @author: ben """ -import h5py import numpy as np import pointCollection as pc -import pyproj + +#TBD: read/write additional formats (parquet?) class data(object): """ @@ -30,7 +30,7 @@ class data(object): -a copy of the object sliced by the index otherwise """ np.seterr(invalid='ignore') - def __init__(self, data_dict=None, fields=None, SRS_proj4=None, EPSG=None, field_dict=None, columns=0, filename=None): + def __init__(self, data_dict=None, fields=None, SRS_proj4=None, EPSG=None, field_dict=None, columns=0, filename=None, coordinates=['x', 'y']): """ Initialize a new pointCollection.data object @@ -51,6 +51,9 @@ def __init__(self, data_dict=None, fields=None, SRS_proj4=None, EPSG=None, field of shape [Ndata, columns]. The default is 0, specifying a shape of [Ndata] filename : str, optional a filename indicating the source of the data. The default is None. + coordinates : list, optional + two-element list giving the names of the x and y coordinate fields. + The default is ['x', 'y']. Returns ------- @@ -62,6 +65,8 @@ def __init__(self, data_dict=None, fields=None, SRS_proj4=None, EPSG=None, field else: self.field_dict=field_dict + self.shape=None + self.size=None if fields is None: fields=list() if field_dict is not None: @@ -69,14 +74,15 @@ def __init__(self, data_dict=None, fields=None, SRS_proj4=None, EPSG=None, field for field in self.field_dict[group]: fields.append(field) if isinstance(fields, dict): - self.fields=list(fields) - self.fields=fields + self.fields=list() + self.assign(fields) + else: + self.fields=fields self.SRS_proj4=SRS_proj4 self.EPSG=EPSG self.columns=columns - self.shape=None - self.size=None self.filename=filename + self.coordinates=list(coordinates) if data_dict is not None: self.assign(data_dict) self.attrs={} @@ -97,13 +103,21 @@ def __default_field_dict__(self): def __copy__(self): other=self.copy_attrs() for field in self.fields: - setattr(other, field, getattr(self, field).copy()) + try: + setattr(other, field, getattr(self, field).copy()) + except AttributeError: + if isinstance(getattr(self, field), (int, float)): + setattr(other, field, getattr(self, field)) return other - def __update_size_and_shape__(self): + def __update_size_and_shape__(self, shape=None): """ When data size and shape may have changed, update the size and shape atttributes """ + if shape is not None: + self.shape=shape + self.size=int(np.prod(shape)) + return self self.size=0 self.shape=[0] for field in self.fields: @@ -117,6 +131,14 @@ def __update_size_and_shape__(self): pass return self + @property + def _x_coord(self): + return self.coordinates[0] + + @property + def _y_coord(self): + return self.coordinates[1] + def __getitem__(self, *args, **kwargs): """ return a field or a subset of self @@ -169,7 +191,7 @@ def choose_crs(self, *args, proj4_string=None, EPSG=None, SRS_proj4=None, SRS_EP elif EPSG is not None: crs=EPSG else: - if self.SRS_proj4 is None: + if self.SRS_proj4 is not None: crs=self.SRS_proj4 else: crs=self.EPSG @@ -186,7 +208,7 @@ def copy_attrs(self): copy attributes to a new data object """ out = type(self)() - for field in ['fields', 'SRS_proj4', 'EPSG', 'columns']: + for field in ['fields', 'SRS_proj4', 'EPSG', 'columns', 'coordinates']: try: temp=getattr(self, field) except AttributeError: @@ -273,12 +295,13 @@ def from_h5(self, filename, group=None, fields=None, field=None, field_dict=None called to fill in the missing fields. """ + import h5py if filename is None: filename=self.filename else: self.filename=filename - if fields is None and fields is not None: + if field is not None and fields is None: fields = [field] if fields is not None: field_dict = {group:fields} @@ -293,9 +316,10 @@ def from_h5(self, filename, group=None, fields=None, field=None, field_dict=None # build the field dict from the group if not isinstance(group, (list, tuple)): group=[group] + field_dict={} for this_group in group: - field_dict={this_group: [key for key in h5_f[this_group].keys() \ - if isinstance(h5_f[this_group][key], h5py.Dataset)]} + field_dict[this_group] = [key for key in h5_f[this_group].keys() \ + if isinstance(h5_f[this_group][key], h5py.Dataset)] else: if self.field_dict is not None: field_dict=self.field_dict @@ -353,6 +377,27 @@ def from_h5(self, filename, group=None, fields=None, field=None, field_dict=None if self.shape is not None: for field in nan_fields: setattr(self, field, np.zeros(self.shape)+np.nan) + # infer coordinates from CF 'coordinates' attribute + coord_names = [] + for grp_name in field_dict.keys(): + if grp_name == '__calc_internal__': + continue + for ds_name in field_dict[grp_name]: + try: + ds = h5_f[ds_name] if grp_name is None else h5_f[grp_name][ds_name] + except KeyError: + continue + coords_attr = ds.attrs.get('coordinates', None) + if coords_attr is not None: + if isinstance(coords_attr, bytes): + coords_attr = coords_attr.decode('utf-8') + for name in coords_attr.split(): + if name not in coord_names: + coord_names.append(name) + break + inferred = [c for c in coord_names if c in self.fields] + if inferred: + self.coordinates = inferred if '__calc_internal__' in field_dict: try: self.__internal_field_calc__(field_dict) @@ -363,17 +408,20 @@ def from_h5(self, filename, group=None, fields=None, field=None, field_dict=None self.__update_size_and_shape__() return self - def get_xy(self, *args, **kwargs): + def get_xy(self, *args, src_coords=None, dst_coords=None, **kwargs): ''' Calculate projected coordinates from latitude and longitude fields x and y fields are calculated based on an object's latitude and - longitude fields. + longitude fields. The object's coordinates will be set to ['x','y'] Parametes *args: iterable list of arguments passed to self.choose_crs() + dst_coords: iterable of strings + these will be used as the names of the output coordinates. + default is 'x','y' **kwargs: dict keyword=pair arguments passed to self.choose_crs() Returns @@ -381,9 +429,9 @@ def get_xy(self, *args, **kwargs): self: pointCollection.data Current object with updated x and y fields ''' + import pyproj crs=self.choose_crs(*args,**kwargs) - try: # this is compatible with pyproj 3.7.0 # Note that EPSG:4326 takes latitude first @@ -399,10 +447,13 @@ def get_xy(self, *args, **kwargs): xy=np.array(pyproj.Proj("+init=epsg:"+str(crs))(self.longitude, self.latitude)) else: xy=np.array(pyproj.Proj(crs)(self.longitude, self.latitude)) - self.x=xy[0,:].reshape(self.shape) - self.y=xy[1,:].reshape(self.shape) - if 'x' not in self.fields: - self.fields += ['x','y'] + if dst_coords is None: + dst_coords = ['x','y'] + + for coord_ind in [0, 1]: + self.assign({ dst_coords[coord_ind] : xy[coord_ind,:].reshape(self.shape)}) + # this may be a problem for data that have one or three coordinates. + self.coordinates = dst_coords return self def get_latlon(self, *args, **kwargs): @@ -425,22 +476,25 @@ def get_latlon(self, *args, **kwargs): Current object with updated latitude and longitude fields ''' + import pyproj crs=self.choose_crs(*args,**kwargs) + _x=getattr(self, self._x_coord) + _y=getattr(self, self._y_coord) try: # Compatible with pyproj 3.7.0 # Note that EPSG:4326 takes latitude first latlon = np.array(pyproj.Transformer.from_crs(pyproj.CRS(crs), - pyproj.CRS(4326)).transform(self.x, self.y)) + pyproj.CRS(4326)).transform(_x, _y)) self.assign({"longitude":latlon[1,:].reshape(self.shape),\ "latitude":latlon[0,:].reshape(self.shape)}) except Exception: if hasattr(pyproj, 'proj'): - lonlat=np.array(pyproj.proj.Proj(crs)(self.x, self.y, inverse=True)) + lonlat=np.array(pyproj.proj.Proj(crs)(_x, _y, inverse=True)) else: if isinstance(crs, int): - lonlat=np.array(pyproj.Proj("+init=epsg:"+str(crs))(self.x, self.y)) + lonlat=np.array(pyproj.Proj("+init=epsg:"+str(crs))(_x, _y, inverse=True)) else: - lonlat=np.array(pyproj.Proj(crs)(self.x, self.y)) + lonlat=np.array(pyproj.Proj(crs)(_x, _y, inverse=True)) self.assign({"longitude":lonlat[0,:].reshape(self.shape), \ "latitude":lonlat[1,:].reshape(self.shape)}) return self @@ -464,7 +518,6 @@ def from_dict(self, dd, fields=None): """ - if fields is not None: self.fields=fields else: @@ -473,7 +526,7 @@ def from_dict(self, dd, fields=None): try: default_shape=dd[next(iter(dd))].shape except StopIteration: - print("HERE!") + raise ValueError("pointCollection.data.from_dict: empty dict") for field in self.fields: if field in dd: setattr(self, field, dd[field]) @@ -605,31 +658,38 @@ def blockmedian(self, scale, field='z'): number of elements, the average of the elements spanning the median is returned. + The median index is computed from a 1-D field. Other fields may be + 1-D or multi-column: for a field of shape (N, M) the index selects + rows, returning shape (K, M). + Parameters ---------- scale : float Scale over which the blockmedian is calculated. field : str, optional - Points within the object are selected based on the values in this - field. The default is 'z'. + 1-D field used to compute the median index. The default is 'z'. + A ValueError is raised if this field is multi-column. Returns ------- pointCollection.data - Object containing the meidan values + Object containing the median values """ - if self.size<2: return self - ind = pc.pt_blockmedian(self.x, self.y, np.float64(getattr(self, field)), scale, return_index=True)[3] - try: - for field in self.fields: - temp_field=getattr(self, field) - setattr(self, field, 0.5*temp_field[ind[:,0]] + 0.5*temp_field[ind[:,1]]) - except IndexError: - pass + ref = getattr(self, field) + if ref.ndim > 1: + raise ValueError( + f"blockmedian: field '{field}' must be 1-D (shape {ref.shape} given)") + _x = getattr(self, self._x_coord) + _y = getattr(self, self._y_coord) if len(self.coordinates) >= 2 else np.zeros_like(_x) + ind = pc.pt_blockmedian(_x, _y, np.float64(ref), scale, return_index=True)[3] + new_fields = {name: 0.5*getattr(self, name)[ind[:,0]] + 0.5*getattr(self, name)[ind[:,1]] + for name in self.fields} + for name, val in new_fields.items(): + setattr(self, name, val) self.__update_size_and_shape__() return self @@ -649,9 +709,13 @@ def complete_fields(self, fields): """ + if self.shape is not None: + default_shape=self.shape + else: + default_shape=[0] for field in fields: if field not in self.fields: - self.assign({field:np.zeros(self.shape)+np.nan}) + self.assign({field:np.zeros(default_shape)+np.nan}) def copy_subset(self, index, by_row=False, datasets=None, fields=None): """ @@ -675,7 +739,7 @@ def copy_subset(self, index, by_row=False, datasets=None, fields=None): """ dd=dict() - if self.columns is not None and self.columns >=1 and by_row is not None or isinstance(index, slice): + if (self.columns is not None and self.columns >= 1) or isinstance(index, slice): by_row=True if fields is not None and datasets is None: datasets=fields @@ -704,20 +768,48 @@ def copy_subset(self, index, by_row=False, datasets=None, fields=None): print("IndexError") return self.copy_attrs().from_dict(dd, fields=datasets) - def cropped(self, bounds, return_index=False, **kwargs): + def _crop_index(self, bounds): + """ + Build a boolean index selecting points within per-dimension bounds. + + Bounds may be given either as one 2-iterable of (min, max) per + dimension, in order matching self.coordinates (e.g. (XR, YR)), or + packed into a single iterable of iterables (e.g. [XR, YR]). + Trailing dimensions may be omitted, and any dimension's bounds may + be None, to leave that dimension uncropped. + """ + if len(bounds) == 1 and isinstance(bounds[0][0], (list, tuple, np.ndarray)): + bounds = bounds[0] + if len(bounds) > len(self.coordinates): + raise ValueError(f"crop: got bounds for {len(bounds)} dimensions, " + f"but self.coordinates has only {len(self.coordinates)}") + bounds = list(bounds) + [None] * (len(self.coordinates) - len(bounds)) + + ind = np.ones(getattr(self, self.coordinates[0]).shape, dtype=bool) + for coord, b in zip(self.coordinates, bounds): + if b is None: + continue + val = getattr(self, coord) + ind &= (val >= b[0]) & (val <= b[1]) + return ind + + def cropped(self, *bounds, return_index=False, **kwargs): """ Return a cropped copy of an object - Returns a copy self for which - bounds[0][0] <= self.x <= bounds[0][1] - and - bounds[1][0] <= self.y <= bounds[1][1] + Returns a copy of self restricted to points for which each + coordinate in self.coordinates falls within the corresponding + range in bounds (e.g. bounds[0][0] <= self.x <= bounds[0][1] and + bounds[1][0] <= self.y <= bounds[1][1]). Parameters ---------- - bounds : iterable - A list of two iterables containing the range of x values and y values - to include in the output. + *bounds : iterable + Bounds for each dimension in self.coordinates, given either as + separate arguments in order (e.g. cropped(XR, YR)) or packed + into a single iterable of iterables (e.g. cropped([XR, YR])). + Trailing dimensions may be omitted, and any dimension's bounds + may be None, to leave that dimension uncropped. return_index : bool, optional If True, return only the indices of the points within the bounds. The default is False. @@ -730,34 +822,27 @@ def cropped(self, bounds, return_index=False, **kwargs): cropped version of the inptu data """ - - ind = (self.x >= bounds[0][0]) & (self.x <= bounds[0][1]) &\ - (self.y >= bounds[1][0]) & (self.y <= bounds[1][1]) + ind = self._crop_index(bounds) if return_index: return ind return self.copy_subset(ind, **kwargs) - def crop(self, bounds): + def crop(self, *bounds): """ Crop current object in place. Parameters ---------- - bounds : iterable - A list of two iterables containing the range of x values and y values - to include in the output. - return_index : bool, optional - If True, return only the indices of the points within the bounds. - The default is False. + *bounds : iterable + Bounds for each dimension in self.coordinates. See `cropped` + for the accepted forms. Returns ------- None """ - - self.index((self.x >= bounds[0][0]) & (self.x <= bounds[0][1]) &\ - (self.y >= bounds[1][0]) & (self.y <= bounds[1][1])) + self.index(self._crop_index(bounds)) def to_h5(self, fileOut=None, h5f_out=None, @@ -806,6 +891,7 @@ def to_h5(self, fileOut=None, None. """ + import h5py # check whether overwriting existing files # append to existing files as default mode = 'w' if replace else 'a' @@ -815,117 +901,121 @@ def to_h5(self, fileOut=None, else: close_file = False - if group is not None: - if not group in h5f_out: - h5f_out.create_group(group.encode('ascii')) - field_dict={} - if meta_dict is None: - field_dict = {field:field for field in self.fields} - else: - field_dict = {} - for out_field, this_md in meta_dict.items(): - if 'source_field' in this_md and this_md['source_field'] is not None and this_md['source_field'] in self.fields: - field_dict[out_field] = this_md['source_field'] - elif out_field in self.fields: - field_dict[out_field] = out_field - - if meta_dict is None: - meta_dict = {out_field:{} for out_field in field_dict.keys()} - - # establish the coordinate fields first - dimension_fields=[] - non_dimension_fields=[] - for out_field in field_dict.keys(): - if 'dimension' in meta_dict[out_field] and meta_dict[out_field]['dimension']:#(out_field==meta_dict[out_field]['dimensions'] or out_field in meta_dict[out_field]['dimensions']): - dimension_fields += [out_field] + try: + if group is not None: + if not group in h5f_out: + h5f_out.create_group(group) + field_dict={} + if meta_dict is None: + field_dict = {field:field for field in self.fields} else: - non_dimension_fields += [out_field] - - for out_field in dimension_fields + non_dimension_fields: - this_data=getattr(self, field_dict[out_field]) - maxshape=this_data.shape - if extensible: - maxshape=list(maxshape) - maxshape[0]=None - # try prepending the 'group' entry of meta_dict to the output field. - # catch the exception thrown if meta_dict is None or if the 'group' - # entry is not there - out_field_name=out_field - try: - out_field_name = meta_dict[out_field]['group'] + out_field - except (TypeError, KeyError) as e: - if DEBUG: - print(e) - out_field_name = group + '/' +out_field - kwargs = dict( compression=compression, - maxshape=tuple(maxshape)) - if meta_dict is not None and 'precision' in meta_dict[out_field] and meta_dict[out_field]['precision'] is not None: - kwargs['scaleoffset']=int(meta_dict[out_field]['precision']) - if 'datatype' in meta_dict[out_field]: - dtype = meta_dict[out_field]['datatype'].lower() - kwargs['dtype']=dtype - if 'int' in dtype: - kwargs['fillvalue'] = np.iinfo(np.dtype(dtype)).max + field_dict = {} + for out_field, this_md in meta_dict.items(): + if 'source_field' in this_md and this_md['source_field'] is not None and this_md['source_field'] in self.fields: + field_dict[out_field] = this_md['source_field'] + elif out_field in self.fields: + field_dict[out_field] = out_field + + if meta_dict is None: + meta_dict = {out_field:{} for out_field in field_dict.keys()} + + # establish the coordinate fields first + dimension_fields=[] + non_dimension_fields=[] + for out_field in field_dict.keys(): + if 'dimension' in meta_dict[out_field] and meta_dict[out_field]['dimension']:#(out_field==meta_dict[out_field]['dimensions'] or out_field in meta_dict[out_field]['dimensions']): + dimension_fields += [out_field] else: - kwargs['fillvalue'] = np.finfo(np.dtype(dtype)).max - # use an explicit '_FillValue' for preference over an overall fillvalue - if '_FillValue' in meta_dict[out_field]: - this_data = np.nan_to_num(this_data,nan=meta_dict[out_field]['_FillValue']) - elif 'fillvalue' in kwargs: - this_data = np.nan_to_num(this_data,nan=kwargs['fillvalue']) - if 'dtype' in kwargs: - this_data = this_data.astype(kwargs['dtype']) - - # Create the dataset - dset = h5f_out.create_dataset(out_field_name.encode('ASCII'), - data=this_data, **kwargs) - if 'dimension' in meta_dict[out_field] and meta_dict[out_field]['dimension']: - dset.make_scale() - if out_field in meta_dict: - for key, val in meta_dict[out_field].items(): - if key.lower() not in ['group','source_field','precision','dimensions']: - if isinstance(val, str): + non_dimension_fields += [out_field] + + for out_field in dimension_fields + non_dimension_fields: + this_data=getattr(self, field_dict[out_field]) + maxshape=this_data.shape + if extensible: + maxshape=list(maxshape) + maxshape[0]=None + # try prepending the 'group' entry of meta_dict to the output field. + # catch the exception thrown if meta_dict is None or if the 'group' + # entry is not there + out_field_name=out_field + try: + out_field_name = meta_dict[out_field]['group'] + '/' + out_field + except (TypeError, KeyError) as e: + if DEBUG: + print(e) + out_field_name = group + '/' + out_field + kwargs = dict( compression=compression, + maxshape=tuple(maxshape)) + if meta_dict is not None and 'precision' in meta_dict[out_field] and meta_dict[out_field]['precision'] is not None: + kwargs['scaleoffset']=int(meta_dict[out_field]['precision']) + if 'datatype' in meta_dict[out_field]: + dtype = meta_dict[out_field]['datatype'].lower() + kwargs['dtype']=dtype + if 'int' in dtype: + kwargs['fillvalue'] = np.iinfo(np.dtype(dtype)).max + else: + kwargs['fillvalue'] = np.finfo(np.dtype(dtype)).max + # use an explicit '_FillValue' for preference over an overall fillvalue + if '_FillValue' in meta_dict[out_field]: + this_data = np.nan_to_num(this_data,nan=meta_dict[out_field]['_FillValue']) + elif 'fillvalue' in kwargs: + this_data = np.nan_to_num(this_data,nan=kwargs['fillvalue']) + if 'dtype' in kwargs: + this_data = this_data.astype(kwargs['dtype']) + + # Create the dataset + dset = h5f_out.create_dataset(out_field_name, + data=this_data, **kwargs) + if 'dimension' in meta_dict[out_field] and meta_dict[out_field]['dimension']: + dset.make_scale() + if self.coordinates and out_field not in self.coordinates: + dset.attrs['coordinates'] = ' '.join(self.coordinates) + if out_field in meta_dict: + for key, val in meta_dict[out_field].items(): + if key.lower() not in ['group','source_field','precision','dimensions']: + if isinstance(val, str): # h5f_out[out_field_name.encode('ASCII')].attrs[key] = str(val).encode('utf-8') - h5f_out[out_field_name.encode('ASCII')].attrs.create(str(key), str(val).encode('utf-8'), None, dtype=' 0: - newdata |= kwargs + newdata = newdata | kwargs + if self.shape is None: + shape=next((v.shape for v in newdata.values() if hasattr(v, 'shape')), None) + if shape is None: + raise ValueError("pointCollection.data.assign: cannot determine shape; no array-valued fields in input") + self.__update_size_and_shape__(shape=shape) for field in newdata.keys(): if isinstance(newdata[field], (int, float)): setattr(self, field, np.zeros(self.shape) + newdata[field]) @@ -995,29 +1091,26 @@ def assign(self, *args, **kwargs): self.fields.append(field) return self - def coords(self): + def coords(self, order=None): """ - Return the coordinates of an object + Return the coordinate arrays of an object, in the order given by + self.coordinates. - Returns the object's coordinates, in (y, x, time) order, or (y, x) - if a 't' or 'time' field is not present. Note that the coordinate - order for this method is not consistent with that in the bounds() - and crop() functions + Parameters + ----------- + order: iterable, optional + order of arrays to return (defaults to self.coordinates) Returns ------- tuple - tuple of object coordinates + tuple of coordinate arrays """ + if order is None: + order = self.coordinates.copy() - - if 'time' in self.fields: - return (self.y, self.x, self.time) - elif 't' in self.fields: - return (self.y, self.x, self.t) - else: - return self.y, self.x + return tuple(getattr(self, c) for c in order) def bounds(self, pad=0): """ @@ -1033,11 +1126,14 @@ def bounds(self, pad=0): XR, YR: minimum and maximum of x and y """ - if len(self.x)==0: - return None, None - - return np.array([np.nanmin(self.x)-pad, np.nanmax(self.x)+pad]), \ - np.array([np.nanmin(self.y)-pad, np.nanmax(self.y)+pad]) + _x = getattr(self, self._x_coord) + if len(_x) == 0: + return (None,) * len(self.coordinates) + result = [] + for c in self.coordinates: + v = getattr(self, c) + result.append(np.array([np.nanmin(v)-pad, np.nanmax(v)+pad])) + return tuple(result) def ravel_fields(self): """ diff --git a/pointCollection/geoIndex.py b/pointCollection/geoIndex.py index 8bfaff5..0f2290d 100644 --- a/pointCollection/geoIndex.py +++ b/pointCollection/geoIndex.py @@ -11,7 +11,6 @@ """ import numpy as np import re -import h5py #from osgeo import osr #import matplotlib.pyplot as plt import pointCollection as pc @@ -199,6 +198,7 @@ def from_file(self, index_file, read_file=False, group='index'): faster than reading the whole file. """ + import h5py h5_f = h5py.File(os.path.expanduser(index_file),'r') h5_i = h5_f[group] if read_file: @@ -232,6 +232,7 @@ def to_file(self, filename): """ write the current geoindex to h5 file 'filename' """ + import h5py indexF = h5py.File(os.path.expanduser(filename),'a', libver='latest') if 'index' in indexF: del indexF['index'] @@ -294,7 +295,7 @@ def for_file(self, filename, file_type, number=0, dir_root='', group=None): if self.DEBUG: raise(e) pass - self.from_list(temp) + self.from_list(temp) if file_type in ['h5']: D=pc.data().from_h5(filename, field_dict={group:['x','y']}) if D.x.size > 0: @@ -332,6 +333,7 @@ def for_file(self, filename, file_type, number=0, dir_root='', group=None): self.attrs['n_files']=1 self.from_xy(xy_bin, filename=filename_out, file_type=file_type, number=number, fake_offset_val=-1) if file_type in ['indexed_h5']: + import h5py h5f=h5py.File(filename,'r') if 'INDEX' in h5f: xy=[np.array(h5f['INDEX']['bin_x']), np.array(h5f['INDEX']['bin_y'])] @@ -362,6 +364,7 @@ def for_file(self, filename, file_type, number=0, dir_root='', group=None): self.attrs['dir_root']=dir_root h5f.close() if file_type in ['indexed_h5_from_matlab']: + import h5py h5f=h5py.File(filename,'r') xy=[np.array(h5f['INDEX']['bin_x']), np.array(h5f['INDEX']['bin_y'])] first_last=None @@ -393,12 +396,13 @@ def query_latlon(self, lat, lon, get_data=True, fields=None, error_action='warn' temp=pc.data().from_dict({'latitude':lat,'longitude':lon}).get_xy(SRS_proj4=self.attrs['SRS_proj4']) x=temp.x y=temp.y - delta=self.attribs['delta'] + delta=self.attrs['delta'] xb=np.round(x/delta[0])*delta[0] yb=np.round(y/delta[1])*delta[1] - return self.query_xy(xb, yb, get_data=get_data, fields=fields, error_action=error_action) + return self.query_xy([xb, yb], get_data=get_data, fields=fields, error_action=error_action) - def query_xy_box(self, xr, yr, get_data=True, fields=None, dir_root='', error_action='warn'): + def query_xy_box(self, xr, yr, get_data=True, fields=None, dir_root='', + full_path=False, error_action='warn'): """ query the current geoIndex for all bins in the box specified by box [xr,yr] """ @@ -406,7 +410,11 @@ def query_xy_box(self, xr, yr, get_data=True, fields=None, dir_root='', error_ac these=(xy_bin[0] >= xr[0]) & (xy_bin[0] <= xr[1]) &\ (xy_bin[1] >= yr[0]) & (xy_bin[1] <= yr[1]) return self.query_xy([xy_bin[0][these], xy_bin[1][these]], get_data=get_data, \ - fields=fields, dir_root=dir_root, bounds=[xr, yr], error_action=error_action) + fields=fields, + dir_root=dir_root, + bounds=[xr, yr], + full_path = full_path, + error_action = error_action) def intersect(self, other, pad=[0, 0]): """ @@ -420,8 +428,15 @@ def intersect(self, other, pad=[0, 0]): other_sub=other.copy_subset(xyBin=[xyB[:,0], xyB[:,1]], pad=pad[1]) return self_sub, other_sub - def query_xy(self, xyb, cleanup=True, get_data=True, fields=None, pad=None, - dir_root='', strict=False, bounds=None, error_action='warn'): + def query_xy(self, xyb, cleanup = True, + get_data = True, + full_path = True, + fields=None, + pad=None, + dir_root='', + strict=False, + bounds=None, + error_action='warn'): """ check if data exist within the current geo index for bins in lists/arrays xb and yb. @@ -489,12 +504,14 @@ def query_xy(self, xyb, cleanup=True, get_data=True, fields=None, pad=None, i1=i1[keep] xy=xy[keep,:] # if the file_N attribute begins with ':', it's a group in the current file, so add the current filename - this_query_file=self.attrs['file_%d' % out_file_num] + this_query_file = self.attrs['file_%d' % out_file_num] if this_query_file is not None and this_query_file[0] == ':': file_base=self.filename.replace(dir_root,'') if 'dir_root' in self.attrs: file_base=file_base.replace(self.attrs['dir_root'],'') this_query_file = file_base + this_query_file + if full_path: + this_query_file = self.resolve_path(this_query_file, dir_root) query_results[this_query_file]={ 'type':self.attrs['type_%d' % out_file_num], 'offset_start':i0, @@ -594,7 +611,7 @@ def get_data(self, query_results, fields=None, data=None, dir_root='', [np.min(all_y)-delta[1]/2, np.max(all_y)+delta[1]/2]] for file_key, result in query_results.items(): - this_file=self.resolve_path(file_key, dir_root) + this_file = self.resolve_path(file_key, dir_root) try: if not (os.path.isfile(this_file) or os.path.isfile(this_file.split(':')[0])): print(f'geoIndex.get_data(): missing file {this_file}') @@ -607,23 +624,21 @@ def get_data(self, query_results, fields=None, data=None, dir_root='', elif result['type'] == 'h5_geoindex': D=geoIndex().from_file(this_file).query_xy((result['x'], result['y']), fields=fields, get_data=True, dir_root=dir_root, error_action=error_action) elif result['type'] == 'ATL06': + this_file, pair = this_file.split(':pair') if fields is None: fields={None:(u'latitude',u'longitude',u'h_li',u'delta_time')} - D6_file, pair=this_file.split(':pair') D=[pc.ATL06.data(beam_pair=int(pair), fields=field_list, field_dict=field_dict).from_h5(\ - filename=D6_file, index_range=np.array(temp)) \ + filename=this_file, index_range=np.array(temp)) \ for temp in zip(result['offset_start'], result['offset_end'])] elif result['type'] == 'ATL11': - D11_file, pair = this_file.split(':pair') + this_file, pair = this_file.split(':pair') try: - if not os.path.isfile(D11_file): - print(D11_file) D=[pc.ATL11.data().from_h5(\ - filename=D11_file, index_range=np.array(temp), \ + filename=this_file, index_range=np.array(temp), \ pair=int(pair), field_dict=field_dict) \ for temp in zip(result['offset_start'], result['offset_end'])] except Exception as e: - print(f"pointCollection.geoIndex: problem with ATL11 file:{D11_file} for beam pair {pair}.") + print(f"pointCollection.geoIndex: problem with ATL11 file:{this_file} for beam pair {pair}.") print(" Indexing information:") print(result) print(" Exception:") @@ -711,7 +726,7 @@ def bins_as_array(self): def bin_latlon(self): xy_bin=self.bins_as_array() - temp=pc.data().from_dict({'x':xy_bin[:,0],'y':xy_bin[:,1]}).get_latlon(SRS_proj4=self.attrs['SRS_proj4']) + temp=pc.data().from_dict({'x':xy_bin[0],'y':xy_bin[1]}).get_latlon(SRS_proj4=self.attrs['SRS_proj4']) #internal_srs=osr.SpatialReference() #internal_srs.ImportFromProj4(self.attrs['SRS_proj4']) diff --git a/pointCollection/grid/data.py b/pointCollection/grid/data.py index 9e2bc46..90fd6fc 100755 --- a/pointCollection/grid/data.py +++ b/pointCollection/grid/data.py @@ -6,7 +6,8 @@ @author: ben """ -from osgeo import gdal, gdalconst, osr, ogr + + import numpy as np import re @@ -15,9 +16,6 @@ import bz2 import gzip import uuid -import h5py -import pyproj -import netCDF4 import posixpath #from scipy.interpolate import RegularGridInterpolator #from scipy.stats import scoreatpercentile @@ -43,7 +41,11 @@ class data(object): as 'upside down' in some plotting routines. """ - def __init__(self, fields=None, fill_value=np.nan, t_axis=2): + def __init__(self, + fields=None, + fill_value=np.nan, + t_axis=2, + coordinates=None): self.x=None self.y=None self.projection=None @@ -55,10 +57,18 @@ def __init__(self, fields=None, fill_value=np.nan, t_axis=2): self.fill_value=fill_value self.size=None self.shape=None - self.t_axis=t_axis + self._t_axis=t_axis self.xform=None self.srs_epsg=None self.srs_proj4=None + if coordinates is not None: + self.coordinates = coordinates + elif self._t_axis==2: + self.coordinates = ['y','x','time'] + elif self._t_axis==0: + self.coordinates = ['time', 'y', 'x'] + else: + self.coordinates = ['y','x'] self._time=None self._spacing=[None, None] if fields is None: @@ -71,19 +81,23 @@ def __init__(self, fields=None, fill_value=np.nan, t_axis=2): @property def spacing(self): try: - return [self.x[1]-self.x[0], self.y[1]-self.y[0]] + col = getattr(self, self._col_coord) + row = getattr(self, self._row_coord) + return [col[1]-col[0], row[1]-row[0]] except Exception: return self._spacing @spacing.setter def spacing(self, val): self._spacing=val return + @property def t(self): return self._time @t.setter def t(self, value): self._time = value + @property def time(self): return self._time @@ -97,14 +111,60 @@ def im_kws(self): return {'extent':self.img_extent, 'origin':'lower'} + @property + def t_axis(self): + return self._t_axis + + @t_axis.setter + def t_axis(self, val): + self._t_axis = val + + @property + def _row_coord(self): + """Name of the row (first spatial) coordinate; default 'y'.""" + if self.coordinates: + if len(self.coordinates) == 3 and self._t_axis == 0: + return self.coordinates[1] + return self.coordinates[0] + return 'y' + + @property + def _col_coord(self): + """Name of the column (second spatial) coordinate; default 'x'.""" + if self.coordinates: + if len(self.coordinates) == 3 and self._t_axis == 0: + return self.coordinates[2] + if len(self.coordinates) >= 2: + return self.coordinates[1] + return 'x' + + @property + def _band_coord(self): + """Name of the band (time-like) coordinate; None for 2-D data.""" + if self.coordinates and len(self.coordinates) == 3: + return self.coordinates[self._t_axis] + return None + def __copy__(self, fields=None): """Return a copy of a grid, optionally with a subset of fields.""" if fields is None: fields=self.fields temp=pc.grid.data() - for field in ['x','y','projection','filename','extent','time', 't', 't_axis']: + for field in ['projection','filename','extent','img_extent', + 'coordinates','_t_axis','xform','srs_epsg','srs_proj4']: if hasattr(self, field): setattr(temp, field, getattr(self, field)) + for coord_name in (self.coordinates or []): + val = getattr(self, coord_name, None) + if val is not None: + setattr(temp, coord_name, val) + for field in ['x', 'y']: + if field not in (self.coordinates or []): + val = getattr(self, field, None) + if val is not None: + setattr(temp, field, val) + if self._time is not None: + temp._time = self._time for field in fields: if hasattr(self, field): setattr(temp, field, getattr(self, field).copy()) @@ -122,66 +182,80 @@ def __repr__(self): def copy(self, fields=None): """Return a copy of the current dataset.""" - return self.__copy__() + return self.__copy__(fields=fields) def copy_meta(self): """Return an empty dataset matching the current dataset.""" temp=pc.grid.data() - for field in ['x','y','projection','filename','extent','time', 't', 't_axis']: + for field in ['projection','filename','extent','img_extent', + 'coordinates','_t_axis','xform','srs_epsg','srs_proj4']: if hasattr(self, field): setattr(temp, field, getattr(self, field)) + for coord_name in (self.coordinates or []): + val = getattr(self, coord_name, None) + if val is not None: + setattr(temp, coord_name, val) + for field in ['x', 'y']: + if field not in (self.coordinates or []): + val = getattr(self, field, None) + if val is not None: + setattr(temp, field, val) + if self._time is not None: + temp._time = self._time temp.__update_size_and_shape__() temp.__update_extent__() return temp - def __getitem__(self, *args, **kwargs): - """Return a subset of the data.""" - return self.copy_subset(*args, **kwargs) + def __getitem__(self, key, **kwargs): + """Return a field by name or a spatial subset of the data.""" + if isinstance(key, str): + if key not in self.fields and not hasattr(self, key): + raise KeyError(f"{key!r} is not a field in this grid") + return getattr(self, key) + return self.copy_subset(key, **kwargs) def __update_extent__(self): """Update the extent of the data to match its x and y fields.""" try: - self.extent=[np.min(self.x), np.max(self.x), np.min(self.y), np.max(self.y)] - if len(self.x)>1: - hdx=np.abs(self.x[1]-self.x[0])/2 + col = getattr(self, self._col_coord) + row = getattr(self, self._row_coord) + self.extent=[np.min(col), np.max(col), np.min(row), np.max(row)] + if len(col)>1: + hdx=np.abs(col[1]-col[0])/2 else: hdx=0 - if len(self.y)>1: - hdy=np.abs(self.y[1]-self.y[0])/2 + if len(row)>1: + hdy=np.abs(row[1]-row[0])/2 else: hdy=0 - self.img_extent=[np.min(self.x)-hdx, np.max(self.x)+hdx, np.min(self.y)-hdy, np.max(self.y)+hdy] - except ValueError: - # usually happens when self.x or self.y is empty + self.img_extent=[np.min(col)-hdx, np.max(col)+hdx, np.min(row)-hdy, np.max(row)+hdy] + except (ValueError, TypeError): + # usually happens when coordinates are empty or None self.extent=[None, None, None, None] def __update_size_and_shape__(self): """Update the size and shape parameters of the object to match that of its data fields.""" - if hasattr(self, 'dimensions') and self.dimensions is not None and self.dimensions[0] is not None: - self.shape=self.dimensions.copy() - self.size=np.prod(self.shape) - return - elif hasattr(self, 'x') and self.x is not None and hasattr(self,'y') and self.y is not None: - self.shape=[len(self.y), len(self.x)] - if hasattr(self, 't') and self.t is not None: - try: - self.shape += [len(self.t)] - except TypeError: - pass - elif hasattr(self, 'time') and self.time is not None: - try: - self.shape += [len(self.time)] - except TypeError: - pass - self.size =np.prod(self.shape) - return - for field in ['z']+self.fields: + + if self.coordinates is not None and self.coordinates[0] is not None: + coords = self.coordinates.copy() + else: + coords = ['y', 'x', 'time'] + self.shape = [] + for coord in coords: try: - self.size=getattr(self, field).size - self.shape=getattr(self, field).shape - break - except Exception: + self.shape.append(len(getattr(self, coord))) + except (AttributeError, TypeError): pass + if not self.shape: + # haven't found any coordinates. Try to get shape from the fields + for field in ['z']+self.fields: + try: + self.size=getattr(self, field).size + self.shape=getattr(self, field).shape + break + except Exception: + pass + self.size=np.prod(self.shape) def summary(self, return_table=True, return_dict=False): """ @@ -233,12 +307,13 @@ def from_dict(self, thedict): Dictionary of spatial grid variables, should contain fields 'x','y', and optionally 'time' or 't' """ + coord_names = set(self.coordinates or []) | {'x', 'y', 'time', 't'} for field in thedict: - if field in ['t','time']: + if field in ('t', 'time'): self._time = thedict[field] else: setattr(self, field, thedict[field]) - if field not in self.fields and field not in ['x','y','time', 't']: + if field not in self.fields and field not in coord_names: self.fields.append(field) self.__update_extent__() self.__update_size_and_shape__() @@ -252,16 +327,19 @@ def assign(self, *args, **kwargs): newdata=dict() if len(kwargs) > 0: newdata |= kwargs - for field in newdata.keys(): - if field in ['t','time']: + coord_names = set(self.coordinates or []) | {'x', 'y', 'time', 't'} + for field, val in newdata.items(): + if isinstance(val, (float, int)): + # allow assign(x=0) to produce x = np.zeros(self.shape) + val = np.zeros(self.shape) + val + if field in ('t', 'time'): self._time=newdata[field] else: - setattr(self, field, newdata[field]) - if field not in self.fields and field not in ['x','y','time','t']: + setattr(self, field, val) + if field not in self.fields and field not in coord_names: self.fields.append(field) return self - def bounds_and_spacing_for_list(self, D_list, spacing_from='last'): xmin,xmax,ymin,ymax = [np.inf,-np.inf,np.inf,-np.inf] spacing = [None]*2 @@ -282,7 +360,7 @@ def bounds_and_spacing_for_list(self, D_list, spacing_from='last'): if spacing[0] is None or spacing_from=='last': spacing =this_spacing elif spacing_from=='min': - spacing = np.minumum(spacing, this_spacing) + spacing = np.minimum(spacing, this_spacing) elif spacing_from=='max': spacing = np.maximum(spacing, this_spacing) @@ -323,12 +401,14 @@ def from_list(self, D_list, t_axis=None, sort=False): # calculate x and y dimensions with new extents nx = np.int64((xmax - xmin)/spacing[0]) + 1 ny = np.int64((ymax - ymin)/spacing[1]) + 1 - # calculate x and y arrays - self.x = np.arange(xmin, xmax+0.1*spacing[0], spacing[0]) - self.y = np.arange(ymin, ymax+0.1*spacing[1], spacing[1]) + # calculate coordinate arrays + setattr(self, self._col_coord, np.arange(xmin, xmax+0.1*spacing[0], spacing[0])) + setattr(self, self._row_coord, np.arange(ymin, ymax+0.1*spacing[1], spacing[1])) # try to extract times time = np.zeros((nt)) i = 0 + # default time field name is 'time' + time_field='time' for D in D_list: try: ntime = D.shape[D.t_axis] @@ -338,7 +418,7 @@ def from_list(self, D_list, t_axis=None, sort=False): if hasattr(D,'time') and D.time is not None: time[i:i+ntime] = D.time time_field='time' - elif hasattr(D,'t'): + elif hasattr(D,'t') and D.t is not None: time[i:i+ntime] = D.t time_field='t' i += ntime @@ -367,8 +447,8 @@ def from_list(self, D_list, t_axis=None, sort=False): try: assert(np.all(np.array(self.spacing)==np.array(D.spacing))) # calculate grid coordinates for merging fields - iy = np.array((D.y[:,None]-ymin)/spacing[1],dtype=int) - ix = np.array((D.x[None,:]-xmin)/spacing[0],dtype=int) + iy = np.array((getattr(D, D._row_coord)[:,None]-ymin)/spacing[1],dtype=int) + ix = np.array((getattr(D, D._col_coord)[None,:]-xmin)/spacing[0],dtype=int) if (D.t_axis == 0): temp = np.zeros((ntime,ny,nx)) temp[:,iy,ix] = getattr(D,field) @@ -378,9 +458,9 @@ def from_list(self, D_list, t_axis=None, sort=False): except AssertionError: if D.t is not None: temp = np.atleast_3d( - D.interp(self.x, self.y, np.atleast_1d(D.t), field=field, gridded=True)) + D.interp(getattr(self, self._col_coord), getattr(self, self._row_coord), np.atleast_1d(D.t), field=field, gridded=True)) else: - temp = np.atleast_3d(D.interp(self.x, self.y, field=field, gridded=True)) + temp = np.atleast_3d(D.interp(getattr(self, self._col_coord), getattr(self, self._row_coord), field=field, gridded=True)) # merge fields if (t_axis == 0) and (D.t_axis == 0): data_field[i:i+ntime,:,:] = temp[:] @@ -448,6 +528,7 @@ def from_geotif(self, file, date_format=None, **kwargs): - ``'year'``: WorldView year from filename - ``'matlab'``: WorldView matlab date from filename """ + from osgeo import gdal, gdalconst self.filename=file if date_format is not None: try: @@ -470,7 +551,8 @@ def from_geotif(self, file, date_format=None, **kwargs): def from_gdal(self, ds, field='z', bands=None, bounds=None, extent=None, skip=1, fill_value=np.nan, min_res=None, meta_only=False, t_range=None, - verbose=False): + group=None, + verbose=False, **kwargs): """ Make a pointCollection.grid.data from a gdal dataset. @@ -497,6 +579,7 @@ def from_gdal(self, ds, field='z', bands=None, bounds=None, extent=None, read bands that have metadata 'time' values between t_range[0] and t_range[1] meta_only : return raster extent without reading data verbose: if true, report errors, etc + group: keyword will be ignored Raises ------ AttributeError @@ -610,6 +693,8 @@ def h5_open(self, h5_file, mode='r', compression=None): - ``'bzip'`` - ``'gzip'`` """ + # lazy import of h5py + import h5py if (compression is None): return h5py.File(h5_file,mode=mode) elif (compression == 'bzip'): @@ -625,37 +710,40 @@ def h5_open(self, h5_file, mode='r', compression=None): fid.seek(0) return h5py.File(fid, 'r') - def select_slices(self, bounds, x, y, bands): + def select_slices(self, bounds, x, y, bands, skip=1): # get orientation of y-axis yorient = np.sign(y[1] - y[0]) + col_key = self._col_coord + row_key = self._row_coord + band_key = self._band_coord or 'time' # reduce raster to bounds and orient to lower slices={} if (bounds is not None) and (yorient > 0): # indices to read xind, = np.nonzero((x >= bounds[0][0]) & (x <= bounds[0][1])) - slices['x'] = slice(xind[0],xind[-1],1) + slices[col_key] = slice(xind[0],xind[-1]+1,skip) yind, = np.nonzero((y >= bounds[1][0]) & (y <= bounds[1][1])) - slices['y'] = slice(yind[0],yind[-1],1) + slices[row_key] = slice(yind[0],yind[-1]+1,skip) elif (bounds is not None) and (yorient < 0): # indices to read with reversed y xind, = np.nonzero((x >= bounds[0][0]) & (x <= bounds[0][1])) - slices['x'] = slice(xind[0],xind[-1],1) + slices[col_key] = slice(xind[0],xind[-1]+1,skip) yind, = np.nonzero((y >= bounds[1][0]) & (y <= bounds[1][1])) - slices['y'] = slice(yind[-1],yind[0],-1) + slices[row_key] = slice(yind[-1],yind[0]-1 if yind[0] > 0 else None,-skip) elif (yorient < 0): # indices to read (all) with reversed y - slices['x'] = slice(None,None,1) - slices['y'] = slice(None,None,-1) + slices[col_key] = slice(None,None,skip) + slices[row_key] = slice(None,None,-skip) else: # indices to read (all) - slices['x'] = slice(None,None,1) - slices['y'] = slice(None,None,1) + slices[col_key] = slice(None,None,skip) + slices[row_key] = slice(None,None,skip) if bands is None: - slices['time'] = slice(None) + slices[band_key] = slice(None) else: - slices['time'] = list(bands) + slices[band_key] = list(bands) return slices, yorient def choose_bands_by_time(self, t=None, bounds=None, t_range=None): @@ -680,6 +768,10 @@ def choose_bands_by_time(self, t=None, bounds=None, t_range=None): range of time values to be read. """ + + if t is None: + t=self.time + bands=None if bounds is not None and len(bounds)==3: if self.t_axis==0: @@ -689,11 +781,35 @@ def choose_bands_by_time(self, t=None, bounds=None, t_range=None): t_range=bounds[2] bounds=bounds[:2] - if t_range is not None: + if t_range is not None and t is not None: bands = np.flatnonzero((t>=t_range[0]) & (t<=t_range[1])) return bands, t_range + def select_fields(self, fields=None): + """ + retain a list of fields in the current object + + Parameters + ---------- + fields : TYPE, iterable + Fields to retain. Fields not in the list will be dropped. The default is None. + + Returns + ------- + None. + + """ + if fields is None: + return + self_fields = self.fields.copy() + for field in set(self_fields): + if field not in fields: + self.fields.remove(field) + try: + delattr(self, field) + except AttributeError: + pass def read_data(self, src, i0, i1, bands): """ @@ -740,7 +856,8 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ meta_only=False, bounds=None, skip=1, fill_value=None, t_axis=None, t_range=None, bands=None, - compression=None, swap_xy=False, source_fillvalue=None): + compression=None, swap_xy=False, source_fillvalue=None, + coord_mapping=None): """ Read a raster from an HDF5 file. @@ -789,7 +906,7 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ pc.grid.data object containing the map data. """ if t_axis is not None: - self.t_axis=t_axis + self.t_axis = t_axis if fill_value is not None: self.fill_value = fill_value if fields is None and field is not None: @@ -797,7 +914,11 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ if field_mapping is None: field_mapping={} self.filename=h5_file - dims=[xname,yname,'t','time'] + if coord_mapping is None: + coord_mapping = {xname: self._col_coord, + yname: self._row_coord, + timename: self._band_coord or 'time'} + dims = [xname, yname, 't', 'time'] + list(self.coordinates or []) if not group.startswith('/'): group='/'+group t=None @@ -808,8 +929,8 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ #default yorient=1 with self.h5_open(h5_file, mode='r', compression=compression) as h5f: - x=np.array(h5f[group][xname]).ravel() - y=np.array(h5f[group][yname]).ravel() + x = np.array(h5f[group][xname]).ravel() + y = np.array(h5f[group][yname]).ravel() for time_var_name in set(['time','t', timename]): if time_var_name in h5f[group].keys(): try: @@ -848,9 +969,9 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ if bounds is not None: # indices to read xind, = np.nonzero((x >= bounds[0][0]) & (x <= bounds[0][1])) - cols = slice(xind[0],xind[-1], skip) + cols = slice(xind[0],xind[-1]+1, skip) yind, = np.nonzero((y >= bounds[1][0]) & (y <= bounds[1][1])) - rows = slice(yind[0],yind[-1], skip) + rows = slice(yind[0],yind[-1]+1, skip) else: # indices to read (all) rows = slice(None,None, skip) @@ -873,7 +994,7 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ nT=len(src_t) except TypeError: nT=1 - if t_axis==0: + if self.t_axis == 0: default_shape_3d = [nT] + default_shape_2d else: default_shape_3d = default_shape_2d + [nT] @@ -889,9 +1010,13 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ if f_field.shape in [tuple(default_shape_3d), tuple(default_shape_2d) ]: z=self.read_data(f_field, i0, i1, bands) elif src_t is not None and len(f_field)==np.prod(default_shape_3d): + # special case: 3d data have been raveled. + # Read the whole dataset and reshape to 3d array z=self.read_data(np.array(f_field).reshape(default_shape_3d), i0, i1, bands) elif len(f_field)==np.prod(default_shape_2d): - z=self.read_data(np.array(f_field).reshape(default_shape_3d), i0, i1, bands) + # special case: 2d data have been raveled. + # Read the whole dataset and reshape to 2d array + z=self.read_data(np.array(f_field).reshape(default_shape_2d), i0, i1, bands) else: raise(IndexError(f'from filename {h5_file}, field {f_field_name} has shape:{f_field.shape} incompatible with data shape:{default_shape_3d}.')) @@ -933,13 +1058,18 @@ def from_h5(self, h5_file, field_mapping=None, group='/', \ if self_field not in self.fields: self.fields.append(self_field) - self.x=x[cols] - self.y=y[rows] + setattr(self, coord_mapping.get(xname, self._col_coord), x[cols]) + row_arr = y[rows] if yorient==-1: - y=y[::-1] + row_arr = row_arr[::-1] + setattr(self, coord_mapping.get(yname, self._row_coord), row_arr) if t is not None: - self.t=t + t_dest = coord_mapping.get(timename, self._band_coord or 'time') + if t_dest in ('t', 'time'): + self._time = t + else: + setattr(self, t_dest, t) # try to retrieve grid mapping and add to projection if grid_mapping_name is not None: self.projection = {} @@ -970,6 +1100,7 @@ def nc_open(self, nc_file, mode='r', compression=None): - ``'bzip'`` - ``'gzip'`` """ + import netCDF4 if (compression is None): return netCDF4.Dataset(nc_file, mode=mode) elif (compression == 'bzip'): @@ -1036,7 +1167,7 @@ def from_nc(self, nc_file, field_mapping=None, group='', self pc.grid.data object containing the map data. """ - dim_names={xname:'x',yname:'y', timename:'time'} + dim_names={xname: self._col_coord, yname: self._row_coord, timename: self._band_coord or 'time'} if t_axis is not None: self.t_axis=t_axis @@ -1052,7 +1183,7 @@ def from_nc(self, nc_file, field_mapping=None, group='', if field_mapping is None: field_mapping={} self.filename=nc_file - dims=[xname,yname,'t','time'] + dims=[xname, yname, 't', 'time'] + list(self.coordinates or []) t=None grid_mapping_name = None with self.nc_open(nc_file,mode='r',compression=compression) as fileID: @@ -1068,7 +1199,7 @@ def from_nc(self, nc_file, field_mapping=None, group='', t=ncf.variables[this_time_var_name][:].copy() timename=this_time_var_name break - dim_names={xname:'x',yname:'y', timename:'time'} + dim_names={xname: self._col_coord, yname: self._row_coord, timename: self._band_coord or 'time'} if bands is None: bands, t_range = self.choose_bands_by_time(t=t, bounds=bounds, t_range=t_range) @@ -1092,9 +1223,9 @@ def from_nc(self, nc_file, field_mapping=None, group='', if hasattr(var,'shape') and var.shape: field_mapping.update({key:key}) - slices=self.select_slices(bounds, x, y, bands)[0] + slices=self.select_slices(bounds, x, y, bands, skip=skip)[0] # check that raster can be sliced - if len(x[slices['x']]) == 0 or len(y[slices['y']]) == 0: + if len(x[slices[self._col_coord]]) == 0 or len(y[slices[self._row_coord]]) == 0: self.__update_extent__() self.__update_size_and_shape__() return self @@ -1108,13 +1239,12 @@ def from_nc(self, nc_file, field_mapping=None, group='', # establish the order of output dimensions for this field # based on the time-axis dimension and the number of dimensions if len(f_field.shape)==2: - out_dims=['y','x'] + out_dims=[self._row_coord, self._col_coord] else: - # N.B. In the next lines, just changed timename to 'time' if self.t_axis==2: - out_dims=['y', 'x', 'time'] + out_dims=[self._row_coord, self._col_coord, self._band_coord or 'time'] else: - out_dims=['time','y','x'] + out_dims=[self._band_coord or 'time', self._row_coord, self._col_coord] if f_field.dimensions is None: f_dims = out_dims else: @@ -1146,10 +1276,14 @@ def from_nc(self, nc_file, field_mapping=None, group='', self.assign({self_field:z}) # reduce x and y to bounds - self.x=x[slices['x']] - self.y=y[slices['y']] + setattr(self, self._col_coord, x[slices[self._col_coord]]) + setattr(self, self._row_coord, y[slices[self._row_coord]]) if t is not None: - self.t=t + band_name = self._band_coord or 'time' + if band_name in ('t', 'time'): + self._time = t + else: + setattr(self, band_name, t) # try to retrieve grid mapping and add to projection if grid_mapping_name is not None: self.projection = {} @@ -1162,13 +1296,23 @@ def from_nc(self, nc_file, field_mapping=None, group='', self.__update_size_and_shape__() return self - def to_h5(self, out_file, fields=None, group='/', replace=False, nocompression=False, attributes={}, fill_value=None, overwrite_coords=False, **kwargs): + def to_h5(self, out_file, fields=None, group='/', replace=False, + nocompression=False, + attributes=None, + fill_value=None, + overwrite_coords=False, + VERBOSE=False, + **kwargs): + # lazy import of h5py + import h5py """Write a grid data object to an hdf5 file.""" kwargs.setdefault('srs_proj4', None) kwargs.setdefault('srs_wkt', None) kwargs.setdefault('srs_epsg', None) kwargs.setdefault('dimensions', {}) kwargs.setdefault('grid_mapping_name', 'crs') + if attributes is None: + attributes={} # check whether overwriting existing files # append to existing files as default mode = 'w' if replace else 'a' @@ -1205,11 +1349,12 @@ def to_h5(self, out_file, fields=None, group='/', replace=False, nocompression=F except Exception: pass - for field in ['x','y','time', 't'] + fields: + coord_fields = list(self.coordinates) if self.coordinates else ['x','y','time'] + for field in coord_fields + fields: f_field_name = posixpath.join(group,field) # if field exists, and overwrite_coords is True, overwrite it if field in h5f[group]: - if field in ['x','y','time', 't']: + if field in coord_fields: if overwrite_coords: if hasattr(self, field): temp=getattr(self, field) @@ -1220,7 +1365,7 @@ def to_h5(self, out_file, fields=None, group='/', replace=False, nocompression=F else: #Otherwise, try to create the dataset try: - if nocompression or field in ['x','y','time']: + if nocompression or field in coord_fields: h5f.create_dataset(f_field_name, data=getattr(self, field)) else: h5f.create_dataset(f_field_name, data=getattr(self, field), @@ -1231,22 +1376,27 @@ def to_h5(self, out_file, fields=None, group='/', replace=False, nocompression=F try: for att_name,att_val in attributes[field].items(): h5f[f_field_name].attrs[att_name] = att_val - except Exception: + except Exception as exc: + if VERBOSE: + print(f"failed to sef field attributes, with exception {exc}") pass # try adding dimensions try: - if field in ['x','y','time', 't'] and any(kwargs['dimensions']): + if field in coord_fields and any(kwargs['dimensions']): h5f[f_field_name].make_scale(field) elif any(kwargs['dimensions']): for i,dim in enumerate(kwargs['dimensions'][f_field_name]): h5f[f_field_name].dims[i].attach_scale(h5f[dim]) except Exception as exc: + if VERBOSE: + print(f"failed to set dimensions with exception {exc}") pass - # add crs attributes if applicable + # add crs attributes if applicable + if self.crs: + # add grid mapping attribute + if field not in coord_fields: + h5f[f_field_name].attrs['grid_mapping'] = kwargs['grid_mapping_name'] if self.crs: - # add grid mapping attribute to each grid field - for field in fields: - h5f[f_field_name].attrs['grid_mapping'] = kwargs['grid_mapping_name'] # add grid mapping variable with projection attributes h5crs = h5f.create_dataset(kwargs['grid_mapping_name'], (), dtype=np.byte) for att_name,att_val in self.crs.items(): @@ -1257,8 +1407,10 @@ def to_h5(self, out_file, fields=None, group='/', replace=False, nocompression=F def to_nc(self, out_file, fields=None, group='', replace=False, time_name='time', - nocompression=False, - attributes={}, fill_value=None, **kwargs): + nocompression = False, + attributes = None, + fill_value =None, + **kwargs): """Write a grid data object to a netCDF4 file.""" kwargs.setdefault('srs_proj4', None) kwargs.setdefault('srs_wkt', None) @@ -1267,6 +1419,9 @@ def to_nc(self, out_file, fields=None, group='', replace=False, time_name='time' kwargs.setdefault('xy_units', 'meter') kwargs.setdefault('t_units', None) + if attributes is None: + attributes = {} + # check whether overwriting existing files # append to existing files as default mode = 'w' if replace else 'a' @@ -1275,8 +1430,12 @@ def to_nc(self, out_file, fields=None, group='', replace=False, time_name='time' fields=self.fields # try getting the time variable name try: - t_name = [field for field in ('t','time') if hasattr(self, field) - and np.any(getattr(self,field))].pop() + band_coord = self._band_coord + if band_coord is not None and getattr(self, band_coord, None) is not None: + t_name = band_coord + else: + t_name = [field for field in ('t', 'time') + if getattr(self, field, None) is not None].pop() except Exception: t_name = None @@ -1284,6 +1443,7 @@ def to_nc(self, out_file, fields=None, group='', replace=False, time_name='time' if fill_value is not None: self.replace_invalid(fields=fields, fill_value=fill_value) + import netCDF4 # get crs attributes self.crs_attributes(**kwargs) with netCDF4.Dataset(out_file,mode) as fileID: @@ -1308,29 +1468,32 @@ def to_nc(self, out_file, fields=None, group='', replace=False, time_name='time' pass # for each dimension variable - for field in ['x','y', time_name]: - # if field exists, overwrite it - if field in ncf.variables.keys() and hasattr(self, field): - var = getattr(self, field) - if var is not None: - ncf.variables[field][:] = var - elif hasattr(self, field): - var = getattr(self, field) - if var is not None: - ncf.createDimension(field, len(np.atleast_1d(var))) - this_var = ncf.createVariable(field, var.dtype, (field,)) - if field in ['x','y']: - this_var.units = kwargs['xy_units'] - if field in ['t','time'] and kwargs['t_units'] is not None: - this_var.units = kwargs['t_units'] - ncf.variables[field][:] = var + row_name = self._row_coord + col_name = self._col_coord + spatial_coord_names = {row_name, col_name, 'x', 'y'} + time_coord_names = {'t', 'time'} + write_coords = [col_name, row_name] + ([t_name] if t_name else []) + for field in write_coords: + var = getattr(self, field, None) + if var is None: + continue + if field in ncf.variables.keys(): + ncf.variables[field][:] = var + else: + ncf.createDimension(field, len(np.atleast_1d(var))) + this_var = ncf.createVariable(field, var.dtype, (field,)) + if field in spatial_coord_names: + this_var.units = kwargs['xy_units'] + if field in time_coord_names and kwargs['t_units'] is not None: + this_var.units = kwargs['t_units'] + ncf.variables[field][:] = var for field in fields: if (self.t_axis == 2) and t_name is not None: - nc_dim=('y','x',t_name,) + nc_dim=(row_name, col_name, t_name,) elif (self.t_axis == 0) and t_name is not None: - nc_dim=(t_name,'y','x',) + nc_dim=(t_name, row_name, col_name,) else: - nc_dim=('y','x',) + nc_dim=(row_name, col_name,) # if field exists, overwrite it if field in ncf.variables.keys() and hasattr(self, field): var = getattr(self, field) @@ -1387,7 +1550,7 @@ def to_geotif(self, out_file, **kwargs): return out_ds def to_gdal(self, driver='MEM', out_file='', field='z', srs_proj4=None, - srs_wkt=None, srs_epsg=None, EPSG=None,dtype=gdal.GDT_Float32, + srs_wkt=None, srs_epsg=None, EPSG=None, dtype=None, options=["compress=LZW"]): """ Write a grid object to a gdal memory object. @@ -1414,10 +1577,15 @@ def to_gdal(self, driver='MEM', out_file='', field='z', srs_proj4=None, GDAL creation options for output raster """ + from osgeo import gdal, gdalconst, osr, ogr + if dtype is None: + dtype = gdal.GDT_Float32 z=np.atleast_3d(getattr(self, field)) ny,nx,nband = [*map(int, z.shape)] - dx=np.abs(np.diff(self.x[0:2]))[0] - dy=np.abs(np.diff(self.y[0:2]))[0] + col = getattr(self, self._col_coord) + row = getattr(self, self._row_coord) + dx=np.abs(np.diff(col[0:2]))[0] + dy=np.abs(np.diff(row[0:2]))[0] # no supported creation options with in memory rasters if driver=='MEM': @@ -1428,7 +1596,7 @@ def to_gdal(self, driver='MEM', out_file='', field='z', srs_proj4=None, # top left x, w-e pixel resolution, rotation # top left y, rotation, n-s pixel resolution - out_ds.SetGeoTransform((self.x.min()-dx/2, dx, 0, self.y.max()+dy/2, 0., -dy)) + out_ds.SetGeoTransform((col.min()-dx/2, dx, 0, row.max()+dy/2, 0., -dy)) if EPSG is not None: srs_epsg=EPSG @@ -1481,43 +1649,40 @@ def as_points(self, fields=None, keep_all=False): fields=self.fields + row_arr = getattr(self, self._row_coord) + col_arr = getattr(self, self._col_coord) + band_name = self._band_coord + band_arr = getattr(self, band_name) if band_name and hasattr(self, band_name) else self.time or self.t if len(self.shape) == 2: - x,y=np.meshgrid(self.x, self.y) + x,y=np.meshgrid(col_arr, row_arr) else: if self.t_axis==0: - if self.time is not None: - t, y, x = np.meshgrid(self.time, self.y, self.x, indexing='ij') - else: - t, y, x = np.meshgrid(self.t, self.y, self.x, indexing='ij') + t, y, x = np.meshgrid(band_arr, row_arr, col_arr, indexing='ij') elif self.t_axis==2: - if self.time is not None: - y, x, t = np.meshgrid(self.y, self.x, self.time, indexing='ij') - else: - y, x, t = np.meshgrid(self.y, self.x, self.t, indexing='ij') + y, x, t = np.meshgrid(row_arr, col_arr, band_arr, indexing='ij') + row_key = self._row_coord + col_key = self._col_coord + band_key = self._band_coord or ('time' if self.time is not None else 't') if keep_all: result = pc.data(filename=self.filename).\ - from_dict({'x':x.ravel(),'y':y.ravel()}) + from_dict({col_key:x.ravel(), row_key:y.ravel()}) for field in fields: result.assign({field:getattr(self, field).ravel()}) else: good=np.isfinite(getattr(self, fields[0])).ravel() result = pc.data(filename=self.filename).\ - from_dict({'x':x.ravel()[good],'y':y.ravel()[good]}) + from_dict({col_key:x.ravel()[good], row_key:y.ravel()[good]}) for field in fields: result.assign({field:getattr(self, field).ravel()[good]}) if len(self.shape)==2: if self.time is not None: result.assign({'time':self.time+np.zeros_like(getattr(result, field))}) else: - if self.time is not None: - time_var='time' - else: - time_var='t' if keep_all: - result.assign({time_var:t.ravel()}) + result.assign({band_key:t.ravel()}) else: - result.assign({time_var:t.ravel()[good]}) + result.assign({band_key:t.ravel()[good]}) return result def get_latlon(self, srs_proj4=None, srs_wkt=None, srs_epsg=None): @@ -1540,6 +1705,7 @@ def get_latlon(self, srs_proj4=None, srs_wkt=None, srs_epsg=None): latitude: float latitude coordinates of grid cells """ + import pyproj # set the spatial projection reference information if srs_proj4 is not None: source = pyproj.CRS.from_proj4(srs_proj4) @@ -1564,7 +1730,7 @@ def add_alpha_band(self, alpha=None, field='z', nodata_vals=None): if alpha is None: if nodata_vals is not None: alpha=np.ones_like(getattr(self, field)[:,:,0]) - if hasattr(nodata_vals, 'len') and len(nodata_vals)==3: + if hasattr(nodata_vals, '__len__') and len(nodata_vals)==3: for ii in range(3): alpha[~np.isfinite(getattr(self, field)[:,:,ii]) | (getattr(self, field)[:,:,ii]==nodata_vals[ii])]=0 elif nodata_vals is not None: @@ -1641,7 +1807,7 @@ def calc_gradient(self, field='z', band=0): else: zz=getattr(self, field) - gy, gx=np.gradient(zz, self.y, self.x) + gy, gx=np.gradient(zz, getattr(self, self._row_coord), getattr(self, self._col_coord)) self.assign({field+'_x':gx, field+'_y':gy}) def toRGB(self, cmap=None, field='z', bands=None, caxis=None, alpha=None): @@ -1685,15 +1851,21 @@ def index(self, row_ind, col_ind, band_ind=None, fields=None ): """Slice a grid by row or column.""" if fields is None: fields=self.fields - self.x=self.x[col_ind] - self.y=self.y[row_ind] + setattr(self, self._col_coord, getattr(self, self._col_coord)[col_ind]) + setattr(self, self._row_coord, getattr(self, self._row_coord)[row_ind]) if band_ind is not None: - for field in ['t','time']: - try: - setattr(self, field, getattr(self, field)[band_ind]) - break - except: - pass + band_name = self._band_coord + if band_name and band_name not in ('t', 'time'): + band_arr = getattr(self, band_name, None) + if band_arr is not None: + setattr(self, band_name, band_arr[band_ind]) + else: + for field in ['t','time']: + try: + setattr(self, field, getattr(self, field)[band_ind]) + break + except: + pass for field in fields: if len(getattr(self, field).shape) == 2: setattr(self, field, getattr(self, field)[row_ind,:][:, col_ind]) @@ -1708,6 +1880,8 @@ def index(self, row_ind, col_ind, band_ind=None, fields=None ): setattr(self, field, getattr(self, field)[:, row_ind,:][:, :, col_ind]) else: setattr(self, field, getattr(self, field)[band_ind,:,:][:, row_ind, :][:, :, col_ind]) + + self.select_fields(fields) self.__update_extent__() self.__update_size_and_shape__() return self @@ -1741,17 +1915,21 @@ def copy_subset(self, rc_ind, band_ind=None, fields=None): return self.copy(fields=fields).index(rc_ind[1], rc_ind[2], band_ind=rc_ind[0]) return self.copy(fields=fields).index(rc_ind[0], rc_ind[1], band_ind=band_ind) - def crop(self, XR, YR, TR=None, fields=None): + def crop(self, XR=None, YR=None, TR=None, fields=None): ''' - Crop self to specified bounds + Crop self to specified bounds inplace Parameters ---------- - XR, YR, TR : iterables - two-element iterables specifying the range in each dimension. TR is optional + XR, iterable + a 2-iterable of numerics specifies the range of x values, + or a 2- or 3-iterable of 2-iterables specifies the range of values + in each dimension + YR, TR : iterables, optional + two-element iterables specifying the range in each dimension. fields : iterable, optional - strings specifying fields to include in the output. The default is None. + strings specifying fields to include in the output. The default is to include all fields. Raises ------ @@ -1759,19 +1937,32 @@ def crop(self, XR, YR, TR=None, fields=None): If TR is specified and neither self.time nor self.t is defined, the subset is not defined + The current object will be cropped in place. To return a cropped copy + of the current object, use the 'cropped' method + Returns ------- pointCollection.grid.data cropped version of self ''' + # Bounds for all dimensions may be packed into a single argument instead + # of passed separately, e.g. crop([XR, YR]) or crop([XR, YR, TR]). + if XR is not None and YR is None and isinstance(XR[0], (list, tuple, np.ndarray)): + if len(XR) == 2: + XR, YR = XR + elif len(XR) == 3: + XR, YR, TR = XR + else: + raise ValueError('packed bounds must contain 2 (XR, YR) or 3 (XR, YR, TR) entries') + if XR is not None: - col_ind = np.flatnonzero((self.x >= XR[0]) & (self.x <= XR[1])) + col_ind = np.flatnonzero((getattr(self, self._col_coord) >= XR[0]) & (getattr(self, self._col_coord) <= XR[1])) else: col_ind=slice(None) if YR is not None: - row_ind = np.flatnonzero((self.y >= YR[0]) & (self.y <= YR[1])) + row_ind = np.flatnonzero((getattr(self, self._row_coord) >= YR[0]) & (getattr(self, self._row_coord) <= YR[1])) else: row_ind=slice(None) @@ -1827,7 +2018,7 @@ def show(self, field='z', band=None, ax=None, xy_scale=1, gradient=False, ddt=No zz /= (t[ddt[1]]-t[ddt[0]]) if gradient: - zz=np.gradient(zz.squeeze(), self.x[1]-self.x[0], self.y[1]-self.y[0])[0] + zz=np.gradient(zz.squeeze(), getattr(self, self._col_coord)[1]-getattr(self, self._col_coord)[0], getattr(self, self._row_coord)[1]-getattr(self, self._row_coord)[0])[0] if 'stretch_pct' not in kwargs and 'clim' not in kwargs: stretch_pct=[5, 95] if 'cmap' not in kwargs: @@ -1881,13 +2072,17 @@ def interp(self, x, y, t=None, gridded=False, band=None, field='z', replace=Fals z0 = getattr(self, field).copy() else: z0 = getattr(self, field).copy() + row = getattr(self, self._row_coord) + col = getattr(self, self._col_coord) + band_name = self._band_coord + band_arr = getattr(self, band_name) if band_name else self.t if len(getattr(self, field).shape)==2 or band is not None: - grid_vars=(self.y, self.x) + grid_vars=(row, col) else: if self.t_axis==0: - grid_vars=(self.t, self.y, self.x) + grid_vars=(band_arr, row, col) else: - grid_vars=(self.y, self.x, self.t) + grid_vars=(row, col, band_arr) self.interpolator[field] = scipy.interpolate.RegularGridInterpolator( grid_vars, z0, bounds_error=False) @@ -1927,7 +2122,9 @@ def bounds(self, pad=0): """ - return [[np.min(self.x)-pad, np.max(self.x)+pad], [np.min(self.y)-pad, np.max(self.y)+pad]] + col = getattr(self, self._col_coord) + row = getattr(self, self._row_coord) + return [[np.min(col)-pad, np.max(col)+pad], [np.min(row)-pad, np.max(row)+pad]] def boundary(self, type='image'): """ @@ -1978,12 +2175,12 @@ def fill_smoothed(self, w_smooth=1, fields=None): def rasterize_poly(self, in_geom, field='z', burn_value=1, epsg=1024, raster_epsg=None, poly_epsg=None): """Rasterize a shapely polygon""" - + from osgeo import gdal, osr, ogr # use whatever epsgs are specified, if none, use epsg 1024 (generic xy cartesian) if raster_epsg is None: if self.srs_epsg is not None: raster_epsg = self.srs_epsg - if poly_epsg is not None: + elif poly_epsg is not None: raster_epsg = poly_epsg else: raster_epsg = epsg @@ -2032,8 +2229,12 @@ def crs_attributes(self, srs_proj4=None, srs_wkt=None, srs_epsg=None, **kwargs): srs_epsg: int or NoneType, default None EPSG projection code """ + from osgeo import osr # output projection attributes dictionary self.crs = {} + # return early if no projection was specified (avoids touching osr) + if srs_proj4 is None and srs_wkt is None and srs_epsg is None: + return # set the spatial projection reference information sr = osr.SpatialReference() if srs_proj4 is not None: @@ -2042,8 +2243,6 @@ def crs_attributes(self, srs_proj4=None, srs_wkt=None, srs_epsg=None, **kwargs): sr.ImportFromWkt(srs_wkt) elif srs_epsg is not None: sr.ImportFromEPSG(srs_epsg) - else: - return # convert proj4 string to dictionary proj4_dict = {p[0]:p[1] for p in re.findall(r'\+(.*?)\=([^ ]+)',sr.ExportToProj4())} # get projection attributes diff --git a/pointCollection/grid/mosaic.py b/pointCollection/grid/mosaic.py index a6c5d86..ee43807 100644 --- a/pointCollection/grid/mosaic.py +++ b/pointCollection/grid/mosaic.py @@ -122,13 +122,14 @@ def image_coordinates(self, temp, validate=False): def setup_bounds_from_list(self, in_list, group=None, fields=None, - bounds=None): + bounds=None, + bands=None): for item in in_list.copy(): if isinstance(item, str): # read tile grid from file try: - temp=pc.grid.data().from_file(item, group=group, meta_only=True) + temp=pc.grid.data().from_file(item, group=group, meta_only=True, bands=bands) if bounds is not None: temp=temp.cropped(*bounds) if temp is not None and (len(temp.x)>0) and (len(temp.y) > 0): @@ -142,6 +143,8 @@ def setup_bounds_from_list(self, in_list, print(f"failed to read group {group} "+ str(item)) in_list.remove(item) else: + if bands is not None: + item=item[:,:,bands] if bounds is not None: item=item.cropped(*bounds) if item is not None and (len(item.x)>0) and (len(item.y) > 0): @@ -157,7 +160,7 @@ def setup_bounds_from_list(self, in_list, self.__update_size_and_shape__() - def setup_fields(self, item, group=None, fields=None): + def setup_fields(self, item, group=None, fields=None, bands=None): ''' Set up fields based on an input data item ''' @@ -166,8 +169,10 @@ def setup_fields(self, item, group=None, fields=None): # read data grid from the first tile HDF5, use it to set the field dimensions if isinstance(item, str): - prototype=pc.grid.mosaic().from_file(item, group=group, fields=fields) + prototype=pc.grid.mosaic().from_file(item, group=group, fields=fields, bands=bands) else: + if bands is not None: + item=item[:,:,bands] prototype=item if len(prototype.fields) == 0: message = f"pointCollection.grid.mosaic.py: did not find fields {fields} in file {prototype.filename}" @@ -190,7 +195,7 @@ def setup_fields(self, item, group=None, fields=None): self.__update_size_and_shape__() - def replace(self, item, group=None, fields=None): + def replace(self, item, group=None, fields=None, bands=None): """ Overwrite a section of the mosaic with an input file or mosaic """ @@ -199,8 +204,10 @@ def replace(self, item, group=None, fields=None): fields=self.fields.copy() # read data grid from HDF5 if isinstance(item, str): - temp=pc.grid.mosaic().from_file(item, group=group, fields=fields) + temp=pc.grid.mosaic().from_file(item, group=group, fields=fields, bands=bands) else: + if bands is not None: + item=item[:,:,bands] temp=item if not temp.overlaps(self): return @@ -230,7 +237,7 @@ def replace(self, item, group=None, fields=None): def add(self, item, fields, group=None, use_time=False, - pad=0, feather=0): + pad=0, feather=0, bands=None): """ Add all bands from an item to a mosaic. @@ -251,6 +258,8 @@ def add(self, item, fields, group=None, use_time=False, pad the weights by this distance. The default is 0. feather : float optional feather weights by this distance. The default is 0. + bands : iterable, optional + Bands to read from the item. The default is None (read all bands). Returns ------- @@ -264,8 +273,10 @@ def add(self, item, fields, group=None, use_time=False, self.normalized=False # read data grid from file if isinstance(item, str): - temp=pc.grid.mosaic().from_file(item, group=group, fields=fields) + temp=pc.grid.mosaic().from_file(item, group=group, fields=fields, bands=bands) else: + if bands is not None: + item=item[:,:,bands] if isinstance(item, self.__class__): temp=item else: @@ -535,6 +546,7 @@ def from_list(self, in_list, by_band=True, verbose=False, spacing=[None, None], + bands=None, ): """ Generate a mosaic from a list of inputs. @@ -560,6 +572,9 @@ def from_list(self, in_list, Smooth blending length for inputs. The default is 0. group : str, optional group in hdf5 or netcdf4 files to read. The default is '/'. + bands : iterable, optional + Bands (e.g. time slices) to read from each input, in order. If + not specified, all bands in each input are read. The default is None. Returns ------- @@ -569,8 +584,8 @@ def from_list(self, in_list, """ weight = (pad is not None and pad > 0) or (feather is not None and feather>0) - self.setup_bounds_from_list(in_list, group=group, fields=fields, bounds=bounds) - message = self.setup_fields(in_list[0], group=group, fields=fields) + self.setup_bounds_from_list(in_list, group=group, fields=fields, bounds=bounds, bands=bands) + message = self.setup_fields(in_list[0], group=group, fields=fields, bands=bands) if message is not None: return message # check if using a weighted summation scheme for calculating mosaic @@ -583,8 +598,11 @@ def from_list(self, in_list, for band in band_list: self.invalid = np.ones(self.dimensions[0:2],dtype=bool) self.weight = np.zeros((self.dimensions[0],self.dimensions[1])) + # if specific input bands were requested, map the output + # band index back to the corresponding input band + in_band = bands[band] if (bands is not None and band is not None) else band for item in in_list: - self.add_to_band(item, group=group, fields=fields, pad=pad, feather=feather, band=band) + self.add_to_band(item, group=group, fields=fields, pad=pad, feather=feather, in_band=in_band, out_band=band) self.normalize(band=band) else: self.invalid = np.ones(self.dimensions[0:2],dtype=bool) @@ -592,7 +610,7 @@ def from_list(self, in_list, # for each file in the list for item in in_list: try: - self.add(item, group=group, fields=fields, pad=pad, feather=feather) + self.add(item, group=group, fields=fields, pad=pad, feather=feather, bands=bands) except Exception as e: print(f"mosaic.from_list : problem with {item} for group={group} and fields={fields}") print(e) @@ -602,7 +620,7 @@ def from_list(self, in_list, # for each file in the list self.invalid = np.ones(self.dimensions[0:2],dtype=bool) for item in in_list: - self.replace(item, group=group, fields=fields) + self.replace(item, group=group, fields=fields, bands=bands) self.normalize(by_weight=False) return self diff --git a/pointCollection/dataPicker.py b/pointCollection/old/dataPicker.py similarity index 100% rename from pointCollection/dataPicker.py rename to pointCollection/old/dataPicker.py diff --git a/pointCollection/points_to_grid.py b/pointCollection/points_to_grid.py index bb3cc26..b59bae6 100755 --- a/pointCollection/points_to_grid.py +++ b/pointCollection/points_to_grid.py @@ -57,7 +57,9 @@ def apply_bin_fn(D_pt, res, fn=None, fields=['z'], xy0=[0, 0]): return result -def points_to_grid(D_pt, res=None, grid=None, field='z', background=np.nan): +def points_to_grid(D_pt, res=None, grid=None, field='z', + out_field='z', + background=np.nan): if grid is not None: res=np.abs(grid.x[1]-grid.x[0]) x=np.round(D_pt.x/res)*res @@ -75,13 +77,16 @@ def points_to_grid(D_pt, res=None, grid=None, field='z', background=np.nan): else: XR=grid.extent[0:2] YR=grid.extent[2:4] - zg=grid.z.copy() + if 'z' in grid.fields: + zg = grid.z.copy() + else: + zg = np.zeros(grid.shape) + background c=((x-XR[0])/res).astype(int) r=((y-YR[0])/res).astype(int) ii = (c>=0) & (c=0) & (r= 0.5*self.tile_spacing - tol) xy[dim] = np.append(xy[dim], ctrs[bdry_ind] + sgn*(self.tile_spacing/2 + tol), axis=0) xy[other_dim] = np.append(xy[other_dim], xy[other_dim][bdry_ind]) - tile_xys = self.mapping_function( np.c_[xy[0], xy[1]] / self.tile_spacing ) * self.tile_spacing + tile_xys = self.mapping_function( + (np.c_[xy[0], xy[1]]-xy0) / self.tile_spacing ) * self.tile_spacing + xy0 if unique: return np.unique(tile_xys, axis=0) else: return tile_xys def tile_filename(self, xy_t): - return os.path.join(self.directory, self.format_str % tuple([xy_t[0]/self.scale, xy_t[1]/self.scale]))+self.extension + if set(['xmin','xmax','ymin','ymax']) == set(self.format_variables): + var_val = {var:val for var, val in zip(['xmin','xmax','ymin','ymax'], + np.concatenate(self.tile_bounds(xy_t)))} + vals_sorted = [(var_val[var]/self.scale) + for var in self.format_variables] + return os.path.join(self.directory, + self.format_str % tuple(vals_sorted))\ + +self.extension + elif set(['x','y']) == set(self.format_variables): + xy0 = [xy_t[0]/self.scale, xy_t[1]/self.scale] + return os.path.join(self.directory, + self.format_str + % tuple(xy0))+self.extension def filenames_for_xy(self, xy0): if np.isscalar(xy0[0]): @@ -187,5 +213,13 @@ def file_xy(self, filenames=None): for filename in filenames: m = self.format_re.search(os.path.basename(filename)) if m is not None: - xy.append(np.array([*map(float, m.groups())])*self.scale) + if set(['x','y']) == set(self.format_variables): + xy.append(np.array([*map(float, m.groups())])*self.scale) + elif ['xmin','xmax','ymin','ymax'] == self.format_variables: + this_xy =np.array([*map(float, m.groups())])*self.scale + xy.append(np.array(np.mean(this_xy[0:2]), + np.mean(this_xy[2:4]))) return xy + +# example: +# ts = tilingSchema(format_str='z0%d_%d_%d_%d', format_variables=['xmin','xmax','ymin','ymax'], tile_spacing=2.e5, EPSG=3031, tile_offset=[1.e5, 1.e5]) diff --git a/pointCollection/tools/convert_ITRF.py b/pointCollection/tools/convert_ITRF.py index 0752505..c9a080f 100644 --- a/pointCollection/tools/convert_ITRF.py +++ b/pointCollection/tools/convert_ITRF.py @@ -1,6 +1,3 @@ -import pyproj - - """ History 2/2024 : Written by Tyler Sutterley @@ -9,7 +6,6 @@ """ -# PURPOSE: convert heights between reference frames def convert_ITRF(lon, lat, z, tdec, epochs=[2014, 2020], direction='forward'): """Converts height data between ITRF2014 and ITRF2020 @@ -29,25 +25,24 @@ def convert_ITRF(lon, lat, z, tdec, epochs=[2014, 2020], direction='forward'): - ``'forward'``: early to late - ``'inverse'``: late to early """ - # get the transform for converting from ITRF2014 to ITRF2020 - if tuple(epochs)==(2014, 2020): - transform = wgs84_itrf2014_to_wgs84_itrf2020() - elif tuple(epochs)==(2005, 2014): - transform = wgs84_itrf2005_to_wgs84_itrf2014() - elif tuple(epochs)==(2008, 2014): - transform = wgs84_itrf2008_to_wgs84_itrf2014() - elif tuple(epochs)==(2005, 2020): - transform = wgs84_itrf2005_to_wgs84_itrf2020() - elif tuple(epochs)==(2008, 2020): - transform = wgs84_itrf2008_to_wgs84_itrf2020() - - # transform the data to ITRF + import pyproj + if tuple(epochs) == (2014, 2020): + pipeline = wgs84_itrf2014_to_wgs84_itrf2020() + elif tuple(epochs) == (2005, 2014): + pipeline = wgs84_itrf2005_to_wgs84_itrf2014() + elif tuple(epochs) == (2008, 2014): + pipeline = wgs84_itrf2008_to_wgs84_itrf2014() + elif tuple(epochs) == (2005, 2020): + pipeline = wgs84_itrf2005_to_wgs84_itrf2020() + elif tuple(epochs) == (2008, 2020): + pipeline = wgs84_itrf2008_to_wgs84_itrf2020() + transform = pyproj.Transformer.from_pipeline(pipeline) return transform.transform(lon, lat, z, tdec, direction=direction) -# WGS84 Ellipsoid in ITRF2014 to WGS84 Ellipsoid in ITRF2020 + def wgs84_itrf2014_to_wgs84_itrf2020(): - """``pyproj`` transform for WGS84 Ellipsoid in ITRF2014 to WGS84 Ellipsoid in ITRF2020""" - pipeline = """+proj=pipeline + """Pipeline string for WGS84 Ellipsoid in ITRF2014 to WGS84 Ellipsoid in ITRF2020""" + return """+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=0.0014 +y=0.0009 +z=-0.0014 +rx=0 +ry=0 +rz=0 +s=0.00042 @@ -55,13 +50,11 @@ def wgs84_itrf2014_to_wgs84_itrf2020(): +t_epoch=2015 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m""" - return pyproj.Transformer.from_pipeline(pipeline) -# WGS84 Ellipsoid in ITRF2008 to WGS84 Ellipsoid in ITRF2014 def wgs84_itrf2008_to_wgs84_itrf2014(): - """``pyproj`` transform for WGS84 Ellipsoid in ITRF2008 to WGS84 Ellipsoid in ITRF2014""" - pipeline = """+proj=pipeline + """Pipeline string for WGS84 Ellipsoid in ITRF2008 to WGS84 Ellipsoid in ITRF2014""" + return """+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=-0.0016 +y=-0.0019 +z=-0.0024 +rx=0 +ry=0 +rz=0 +s=0.00002 @@ -69,13 +62,11 @@ def wgs84_itrf2008_to_wgs84_itrf2014(): +t_epoch=2015 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m""" - return pyproj.Transformer.from_pipeline(pipeline) -# WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2014 def wgs84_itrf2005_to_wgs84_itrf2014(): - """``pyproj`` transform for WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2014""" - pipeline = """+proj=pipeline + """Pipeline string for WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2014""" + return """+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=-0.0026 +y=-0.001 +z=0.0023 +rx=0 +ry=0 +rz=0 +s=-0.00092 @@ -83,13 +74,11 @@ def wgs84_itrf2005_to_wgs84_itrf2014(): +t_epoch=2010 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m""" - return pyproj.Transformer.from_pipeline(pipeline) -# WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2020 def wgs84_itrf2005_to_wgs84_itrf2020(): - """``pyproj`` transform for WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2020""" - pipeline = """+proj=pipeline + """Pipeline string for WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2020""" + return """+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=-0.0027 +y=-0.0001 +z=0.0014 +rx=0 +ry=0 +rz=0 +s=-0.00065 @@ -97,12 +86,11 @@ def wgs84_itrf2005_to_wgs84_itrf2020(): +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m""" - return pyproj.Transformer.from_pipeline(pipeline) -# WGS84 Ellipsoid in ITRF2008 to WGS84 Ellipsoid in ITRF2020 + def wgs84_itrf2008_to_wgs84_itrf2020(): - """``pyproj`` transform for WGS84 Ellipsoid in ITRF2005 to WGS84 Ellipsoid in ITRF2014""" - pipeline = """+proj=pipeline + """Pipeline string for WGS84 Ellipsoid in ITRF2008 to WGS84 Ellipsoid in ITRF2020""" + return """+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=WGS84 +step +proj=helmert +x=-0.0002 +y=-0.001 +z=-0.0033 +rx=0 +ry=0 +rz=0 +s=0.00029 @@ -110,4 +98,3 @@ def wgs84_itrf2008_to_wgs84_itrf2020(): +t_epoch=2015 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m""" - return pyproj.Transformer.from_pipeline(pipeline) diff --git a/test/test_geoindex.py b/test/test_geoindex.py deleted file mode 100755 index 2f87a54..0000000 --- a/test/test_geoindex.py +++ /dev/null @@ -1,25 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Sat Apr 4 08:32:55 2020 - -@author: ben -""" - -import os -import warnings -import pytest -import pointCollection as pc -import numpy as np - -def test_geoindex(): - xy0=(np.array(-360000.), np.array(-1140000.)) - geoIndex_file = '/Volumes/ice2/ben/CS2_data/GL/retrack_BC/2019/index_POCA.h5' - # only run test if file exists - try: - gi=pc.geoIndex().from_file(geoIndex_file) - except OSError: - print('{0} not found on file system'.format(geoIndex_file)) - else: - D=gi.query_xy(xy0, fields=['x','y', 'h']) - gi=None diff --git a/test/test_point.py b/test/test_point.py deleted file mode 100644 index 89a053d..0000000 --- a/test/test_point.py +++ /dev/null @@ -1,39 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Sat Apr 4 08:32:55 2020 - -@author: ben -""" -import os -import inspect -import warnings -import pytest -import numpy as np -import pointCollection as pc - -def test_point(): - x=np.arange(10) - y=x**2 - z=y - D1 = pc.data().from_dict({'x':x,'y':y,'z':z}) - D2 = pc.data().from_list([D1]) - assert all((getattr(D1,key) == getattr(D2,key)).all() for key in D1.fields) - -def test_read_ATL06(): - filename = inspect.getframeinfo(inspect.currentframe()).filename - filepath = os.path.dirname(os.path.abspath(filename)) - ATL06_file = os.path.join(filepath,'..','test_data', - 'ATL06_20190205041106_05910210_209_01.h5') - # general list of ATL06 fields with projected coordinate fields - ATL06_fields = ['delta_time','h_li','h_li_sigma','latitude','longitude', - 'atl06_quality_summary','segment_id','sigma_geo_h','x_atc','y_atc', - 'seg_azimuth','sigma_geo_at','sigma_geo_xt','dh_fit_dx','dh_fit_dx_sigma', - 'h_mean','dh_fit_dy','h_rms_misfit','h_robust_sprd','n_fit_photons', - 'signal_selection_source','snr_significance','w_surface_window_final', - 'bsnow_conf','bsnow_h','cloud_flg_asr','cloud_flg_atm','r_eff', - 'tide_ocean','valid','BP','LR','cycle_number','rgt','x','y'] - # read ATL06 file and get projected coordinates - D_06 = pc.ATL06.data().from_h5(os.path.relpath(ATL06_file)) - D_06.get_xy(EPSG=3031) - assert (D_06.fields == ATL06_fields) diff --git a/tests/test_data.py b/tests/test_data.py new file mode 100644 index 0000000..33bdfc5 --- /dev/null +++ b/tests/test_data.py @@ -0,0 +1,288 @@ +""" +Tests for pointCollection.data +""" +import os +import numpy as np +import pytest +import pointCollection as pc + +TEST_H5 = os.path.join(os.path.dirname(__file__), '..', 'test_data', + 'ATL06_20190205041106_05910210_209_01.h5') + + +# --------------------------------------------------------------------------- +# __init__ / assign +# --------------------------------------------------------------------------- + +def test_fields_dict_assigns_data(): + x = np.array([1.0, 2.0, 3.0]) + y = np.array([4.0, 5.0, 6.0]) + D = pc.data(fields={'x': x, 'y': y}) + assert isinstance(D.fields, list) + assert set(D.fields) == {'x', 'y'} + np.testing.assert_array_equal(D.x, x) + np.testing.assert_array_equal(D.y, y) + + +def test_fields_dict_is_list_not_dict(): + D = pc.data(fields={'x': np.zeros(3), 'y': np.zeros(3)}) + _ = D.fields[0] # would fail if fields were a dict + assert len(D.fields) == 2 + + +def test_fields_list_unchanged(): + fields = ['x', 'y', 'z'] + D = pc.data(fields=fields) + assert D.fields == fields + + +def test_fields_none_default(): + D = pc.data() + assert D.fields == [] + + +def test_assign_sets_shape_from_array(): + D = pc.data() + D.assign({'x': np.arange(5.0), 'y': np.arange(5.0)}) + assert D.shape == (5,) + assert D.size == 5 + + +def test_assign_scalar_with_array_broadcasts(): + D = pc.data() + D.assign({'x': np.arange(4.0), 'scale': 2.0}) + assert D.shape == (4,) + np.testing.assert_array_equal(D.scale, np.full(4, 2.0)) + + +def test_assign_all_scalars_no_shape_raises(): + D = pc.data() + with pytest.raises(ValueError, match="cannot determine shape"): + D.assign({'scale': 1.0}) + + +def test_assign_kwargs_does_not_mutate_caller_dict(): + # |= on a dict is in-place; using | creates a new dict so the caller's dict is unchanged + D = pc.data(fields={'x': np.arange(4.0)}) + original = {'y': np.arange(4.0)} + D.assign(original, z=np.ones(4)) + assert 'z' not in original + + +# --------------------------------------------------------------------------- +# __update_size_and_shape__ +# --------------------------------------------------------------------------- + +def test_update_size_and_shape_explicit(): + D = pc.data() + D.__update_size_and_shape__(shape=(10, 3)) + assert D.shape == (10, 3) + assert D.size == 30 + + +# --------------------------------------------------------------------------- +# from_h5 with multiple groups +# --------------------------------------------------------------------------- + +def test_from_h5_multiple_groups_reads_all(): + groups = ['gt1l/land_ice_segments/dem', 'gt1l/land_ice_segments/fit_statistics'] + D = pc.data().from_h5(TEST_H5, group=groups) + assert 'dem_h' in D.fields + assert 'h_mean' in D.fields + assert hasattr(D, 'dem_h') + assert hasattr(D, 'h_mean') + assert D.dem_h.shape == D.h_mean.shape + + +# --------------------------------------------------------------------------- +# copy_subset by_row behaviour +# --------------------------------------------------------------------------- + +def test_copy_subset_by_row_true_slices_rows(): + D = pc.data(fields={'x': np.arange(6.0).reshape(3, 2)}, columns=2) + sub = D.copy_subset(np.array([0, 2]), by_row=True) + np.testing.assert_array_equal(sub.x, np.array([[0., 1.], [4., 5.]])) + + +def test_copy_subset_columns_auto_sets_by_row(): + # columns >= 1 forces row-wise slicing regardless of the by_row default + D = pc.data(fields={'x': np.arange(6.0).reshape(3, 2)}, columns=2) + sub = D.copy_subset(np.array([0, 2])) + np.testing.assert_array_equal(sub.x, np.array([[0., 1.], [4., 5.]])) + + +def test_copy_subset_1d_ravels(): + # 1-D fields use element-wise indexing (no by_row auto-override) + D = pc.data(fields={'x': np.arange(6.0)}) + sub = D.copy_subset(np.array([0, 2, 4])) + np.testing.assert_array_equal(sub.x, np.array([0., 2., 4.])) + + +# --------------------------------------------------------------------------- +# get_xy / get_latlon roundtrip +# --------------------------------------------------------------------------- + +def test_blockmedian_raises_on_multicolumn_field(): + D = pc.data(fields={ + 'x': np.array([0., 0., 1., 1.]), + 'y': np.array([0., 0., 0., 0.]), + 'z': np.ones((4, 2)), # multi-column — must not be used as median field + }) + with pytest.raises(ValueError, match="must be 1-D"): + D.blockmedian(1.0, field='z') + + +def test_blockmedian_applies_index_to_multicolumn_fields(): + # Four points in two spatial bins; each bin contains two identical rows. + # After blockmedian the averaged rows should equal the original rows. + x = np.array([0.1, 0.2, 1.1, 1.2]) + y = np.zeros(4) + z = np.array([1.0, 2.0, 3.0, 4.0]) + v = np.array([[1., 10.], [2., 20.], [3., 30.], [4., 40.]]) # multi-column + D = pc.data(fields={'x': x, 'y': y, 'z': z, 'v': v}, columns=2) + D.blockmedian(1.0, field='z') + assert D.v.ndim == 2 + assert D.v.shape[1] == 2 + + +def test_to_h5_attrs_written_to_group(tmp_path): + import h5py + D = pc.data(fields={'x': np.arange(5.0), 'y': np.arange(5.0)}) + D.attrs = {'source': 'test', 'version': 1} + outfile = str(tmp_path / 'out.h5') + D.to_h5(outfile, group='/data') + with h5py.File(outfile, 'r') as f: + assert 'source' in f['/data'].attrs + assert 'version' in f['/data'].attrs + # attrs must NOT appear on a dataset + assert 'source' not in f['/data/x'].attrs + + +def test_coordinates_roundtrip_h5(tmp_path): + import h5py + x = np.arange(5.0) + y = np.arange(5.0) * 2 + z = np.ones(5) + D = pc.data(fields={'x': x, 'y': y, 'z': z}) + outfile = str(tmp_path / 'coords.h5') + D.to_h5(outfile, group='/data') + # verify coordinates attr written on non-coordinate field, not on coordinate fields + with h5py.File(outfile, 'r') as f: + assert f['/data/z'].attrs.get('coordinates') == 'x y' + assert 'coordinates' not in f['/data/x'].attrs + assert 'coordinates' not in f['/data/y'].attrs + # verify coordinates recovered on read + D2 = pc.data().from_h5(outfile, group='/data') + assert D2.coordinates == ['x', 'y'] + + +def test_coordinates_roundtrip_h5_custom(tmp_path): + lon = np.linspace(-10.0, 10.0, 6) + lat = np.linspace(60.0, 70.0, 6) + z = np.ones(6) + D = pc.data(fields={'lon': lon, 'lat': lat, 'z': z}, coordinates=['lon', 'lat']) + outfile = str(tmp_path / 'custom_coords.h5') + D.to_h5(outfile, group='/data') + D2 = pc.data().from_h5(outfile, group='/data') + assert D2.coordinates == ['lon', 'lat'] + + +def test_get_xy_get_latlon_roundtrip(): + lat = np.array([-75.0, -80.0]) + lon = np.array([0.0, 45.0]) + D = pc.data(fields={'latitude': lat, 'longitude': lon}) + D.get_xy(EPSG=3031) + assert 'x' in D.fields and 'y' in D.fields + D.get_latlon(EPSG=3031) + np.testing.assert_allclose(D.latitude, lat, atol=1e-6) + np.testing.assert_allclose(D.longitude, lon, atol=1e-6) + + +# --------------------------------------------------------------------------- +# configurable coordinates +# --------------------------------------------------------------------------- + +def test_custom_coordinates_bounds(): + cx = np.array([1.0, 2.0, 3.0]) + cy = np.array([10.0, 20.0, 30.0]) + D = pc.data(fields={'cx': cx, 'cy': cy}, coordinates=['cx', 'cy']) + xr, yr = D.bounds() + np.testing.assert_array_equal(xr, np.array([1.0, 3.0])) + np.testing.assert_array_equal(yr, np.array([10.0, 30.0])) + + +def test_custom_coordinates_coords(): + cx = np.array([1.0, 2.0]) + cy = np.array([3.0, 4.0]) + D = pc.data(fields={'cx': cx, 'cy': cy}, coordinates=['cx', 'cy']) + result = D.coords() + # coords() returns fields in self.coordinates order + np.testing.assert_array_equal(result[0], cx) + np.testing.assert_array_equal(result[1], cy) + + +def test_coords_custom_order(): + cx = np.array([1.0, 2.0]) + cy = np.array([3.0, 4.0]) + D = pc.data(fields={'cx': cx, 'cy': cy}, coordinates=['cx', 'cy']) + result = D.coords(order=['cy', 'cx']) + np.testing.assert_array_equal(result[0], cy) + np.testing.assert_array_equal(result[1], cx) + + +def test_single_coordinate_coords(): + d = np.array([0.0, 1.0, 2.0]) + D = pc.data(fields={'dist': d}, coordinates=['dist']) + result = D.coords() + assert len(result) == 1 + np.testing.assert_array_equal(result[0], d) + + +def test_single_coordinate_bounds(): + d = np.array([1.0, 3.0, 5.0]) + D = pc.data(fields={'dist': d}, coordinates=['dist']) + result = D.bounds() + assert len(result) == 1 + np.testing.assert_array_equal(result[0], np.array([1.0, 5.0])) + + +def test_single_coordinate_crop(): + d = np.array([0.0, 1.0, 2.0, 3.0]) + z = np.array([10.0, 20.0, 30.0, 40.0]) + D = pc.data(fields={'dist': d, 'z': z}, coordinates=['dist']) + D.crop([[0.5, 2.5]]) + np.testing.assert_array_equal(D.dist, np.array([1.0, 2.0])) + np.testing.assert_array_equal(D.z, np.array([20.0, 30.0])) + + +def test_single_coordinate_blockmedian(): + d = np.array([0.1, 0.2, 1.1, 1.2]) + z = np.array([1.0, 3.0, 5.0, 7.0]) + D = pc.data(fields={'dist': d, 'z': z}, coordinates=['dist']) + D.blockmedian(1.0, field='z') + assert D.size == 2 + + +def test_custom_coordinates_crop(): + cx = np.array([0.0, 1.0, 2.0, 3.0]) + cy = np.array([0.0, 0.0, 0.0, 0.0]) + D = pc.data(fields={'cx': cx, 'cy': cy}, coordinates=['cx', 'cy']) + D.crop([[0.5, 2.5], [-1.0, 1.0]]) + np.testing.assert_array_equal(D.cx, np.array([1.0, 2.0])) + + +def test_coordinates_copied_in_copy_subset(): + cx = np.array([0.0, 1.0, 2.0]) + cy = np.array([5.0, 6.0, 7.0]) + D = pc.data(fields={'cx': cx, 'cy': cy}, coordinates=['cx', 'cy']) + sub = D.copy_subset(np.array([0, 2])) + assert sub.coordinates == ['cx', 'cy'] + xr, yr = sub.bounds() + np.testing.assert_array_equal(xr, np.array([0.0, 2.0])) + + +def test_default_coordinates_are_xy(): + D = pc.data() + assert D.coordinates == ['x', 'y'] + assert D._x_coord == 'x' + assert D._y_coord == 'y'