diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..006cd363ab6dc798ae2aef34d1596f8d3cee1c65 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.pyc +sandbox/ +floenavi/CFGDIR diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000000000000000000000000000000000000..35b4677d51ac3cea12cd06ec0fda50cc0243134a --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2019 Stefan Hendricks + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000000000000000000000000000000000000..1605f1dcdf69c54df1cbf713eef502ac16ecdc34 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,4 @@ +include LICENSE +include floenavi/VERSION +include requirements.txt +recursive-include bin *.py \ No newline at end of file diff --git a/README.rst b/README.rst new file mode 100644 index 0000000000000000000000000000000000000000..97cd3fa4edb97bbb090bb31916b44abdef174a0f --- /dev/null +++ b/README.rst @@ -0,0 +1,3 @@ +FloeNavi +======== + diff --git a/bin/floenavi-set-cfgdir.py b/bin/floenavi-set-cfgdir.py new file mode 100644 index 0000000000000000000000000000000000000000..0d627b1d8a698707ed518402a2e0404df9eb1ad4 --- /dev/null +++ b/bin/floenavi-set-cfgdir.py @@ -0,0 +1,13 @@ +# -*- coding: utf-8 -*- + +""" +Short script to set the config directory for the floenavi package +""" + +import os +from floenavi import PACKAGE_ROOT_DIR + +config_dir_def = os.path.abspath(os.path.join(PACKAGE_ROOT_DIR, "CFGDIR")) +CFG_DIR = input("Enter local directory with floenavi config files: ") +with open(config_dir_def, "w") as f: + f.write(CFG_DIR) diff --git a/bin/floenavi_gpstrack2xy.py b/bin/floenavi_gpstrack2xy.py new file mode 100644 index 0000000000000000000000000000000000000000..8dada736f3fc84e4801d00c67f39881bdc8b1d53 --- /dev/null +++ b/bin/floenavi_gpstrack2xy.py @@ -0,0 +1,85 @@ +# -*- coding: utf-8 -*- + +""" +Objective of the script: Create a gps track with FloeNavi coordinates from a gps track file. + +Usage: + + python floenavi_gpstrack2xy.py <file_path> | <file_search_pattern> + +""" + +# --- Imports --- + +import re +import sys +import glob + +from loguru import logger + +from floenavi import FloeNavi +from floenavi.dataset import GPSTrack + + +# --- Globals --- + +USAGE_STR = "python floenavi_gpstrack2xy.py <file_path> | <file_search_pattern>" + + +def floenavi_gpstrack2xy(): + """ + Workflow script for computing FloeNavi (x,y) coordinates based on a a FloeNavi standard GPS track csv file. + Output will be written in file with additional floenavi tag. + :return: + """ + + # Get the command line argument + try: + file_search = sys.argv[1] + except IndexError: + logger.error("Error: missing command line argument -> Terminating") + logger.info("Usage: %s" % USAGE_STR) + sys.exit(1) + + # Get the list of files + file_paths = sorted(glob.glob(file_search)) + + # Check number of files to process + if len(file_paths) == 0: + logger.error("Error: No such file(s) [%s] -> Terminating" % str(file_search)) + sys.exit(1) + logger.info("Found %d file(s)" % len(file_paths)) + + # Process files + floenavi_app = FloeNavi() + + # Loop over all files + counter = 0 + for file_path in file_paths: + + logger.info("Converting: %s" % file_path) + + # Check if file has the correct naming + if not re.search(r"track\.csv$", file_path): + logger.error("File %s does not meet naming requirements (*track.csv) -> skipping") + continue + + # Read the input file + track = GPSTrack.from_csv(file_path) + + # Compute coordinate transformation + floenavi_app.get_xy_coordinates(track) + + # Write the output file + output_filepath = file_path.replace(".csv", "-floenavi.csv") + track.export_csv(output_filepath) + + # Count successfully converted files + counter += 1 + + # Report and close + logger.info("Converted %d files" % counter) + + +if __name__ == "__main__": + floenavi_gpstrack2xy() diff --git a/floenavi/VERSION b/floenavi/VERSION new file mode 100644 index 0000000000000000000000000000000000000000..e6d5cb833c634c4e2972335825d1100492e7c96b --- /dev/null +++ b/floenavi/VERSION @@ -0,0 +1 @@ +1.0.2 \ No newline at end of file diff --git a/floenavi/__init__.py b/floenavi/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..d2cd46dabe5af7b72223c6f220dd308e1e96011f --- /dev/null +++ b/floenavi/__init__.py @@ -0,0 +1,1183 @@ +# -*- coding: utf-8 -*- + +import os +import sys +import glob +import yaml +import parse +import pyproj +import numpy as np +from datetime import datetime, timedelta +from loguru import logger + +from floenavi.ais import DshipAISExtract + +# The algorithm for sea ice coordinate system is implemented in the icedrift module +from icedrift import GeoReferenceStation, GeoPositionData, IceCoordinateSystem + + +# Get version from VERSION in package root +PACKAGE_ROOT_DIR = os.path.dirname(os.path.abspath(__file__)) +version_file = open(os.path.abspath(os.path.join(PACKAGE_ROOT_DIR, "VERSION"))) +try: + with version_file as f: + version = f.read().strip() +except IOError: + sys.exit("Cannot find VERSION file in package (expected: %s" % version_file) + +# Get the location of the config dir +config_dir_def = os.path.abspath(os.path.join(PACKAGE_ROOT_DIR, "CFGDIR")) +try: + with open(config_dir_def, "r") as f: + CFG_DIR = f.read().strip() +except IOError: + print("Cannot find CFGDIR file in package (expected: `%s`)" % config_dir_def) + CFG_DIR = input("Enter local directory to floenavi config files: ") + with open(config_dir_def, "w") as f: + f.write(CFG_DIR) + print("Written to CFGDIR") + +# Get default config files +DEFAULT_CONFIG_FILE = os.path.abspath(os.path.join(CFG_DIR, "floenavi_cfg.yaml")) +DEFAULT_PATHDEF_FILE = os.path.abspath(os.path.join(CFG_DIR, "floenavi_local_path_def.yaml")) + +# Package Metadata + +__all__ = ["ais", "dataset", # modules + "FloeNavi", "FloeNaviConfig", # main classes + PACKAGE_ROOT_DIR] + +__version__ = version +__author__ = "Stefan Hendricks" +__author_email__ = "stefan.hendricks@awi.de" + + +class FloeNavi(object): + """ + The main FloeNavi class and interface for the ice coordinate transformation for floenavi.datasets instances + + Usage: + + floenavi = FloeNavi() + + dataset = dataset.Waypoints(...) + or dataset = dataset.GPSTrack(...) + + floenavi.get_xy_coordinates(dataset) + + dataset.export_csv(filename) + """ + + def __init__(self, cfg=None): + """ + Main FloeNavi app with optional FloeNavi configuration instance. + :param cfg: (Optional) FloeNaviConfig instance that overrides the default FloeNavi Configuration + """ + + # --- Store arguments --- + + # Use default Floenavi configuration if no specific configuration is specified + if cfg is None: + cfg = FloeNaviConfig() + self.cfg = cfg + + def get_xy_coordinates(self, dset): + """ + Compute the (x, y) of the ice coordinate system for a given data set. The data set will be changed in-place + :param dset: a dataset.WayPoints or dataset.GPSTrack instance + :return: + """ + + # Check temporal coverage of data set + msg = "Dataset time range: [%s - %s]" % tuple([time.isoformat() for time in dset.time_bounds]) + logger.info(msg) + + # Temporal coverage of base station data + # These depend on the temporal coverage of the data set + floenavi_cfg = self.cfg.get_floenavi_cfg(dset.time_bounds) + + if floenavi_cfg is None: + logger.error("No FloeNavi configuration found for selected time range in config file -> terminating") + sys.exit(1) + + logger.info("Requested AIS reference stations: %s" % (", ".join([str(mmsi) for mmsi in floenavi_cfg.mmsis]))) + + # Get the (icedrift) reference station data for the selected base stations and data period + refstat = self.get_virtual_reference_stations(dset.time_bounds, floenavi_cfg) + + # --- Init the ice coordinate system transformation --- + icecs_trns = IceCoordinateSystem(refstat, **floenavi_cfg.icecs_trns_kwargs) + logger.info("Initialized Ice Coordinate System (icecs) transformation") + logger.debug("Settings: %s" % str(floenavi_cfg.icecs_trns_kwargs)) + + # --- Convert the coordinates into ice coordinates --- + pos_data = GeoPositionData(dset.time, dset.longitude, dset.latitude) + ice_pos, icecs = icecs_trns.get_xy_coordinates(pos_data) + logger.info("Computed (x,y) coordinates") + logger.debug("projection parameters: %s" % icecs.proj4str) + + # --- Update the data set in place ---- + dset.update_floenavi_coords(ice_pos.xc, ice_pos.yc) + + def get_virtual_reference_stations(self, dataset_time_bounds, basestation_setup): + """ + Get two virtual reference stations (origin ref and heading ref) from merging + the information of all available FloeNavi reference stations. The step is + necessary as the convention of 1 basestation at grid origin and 1 basestation + at a distance defining a axis of the grid could not be maintained during + MOSAiC (had to be moved where there was power available). Since the icedrift + package requires an origin and heading reference station, this two are computed + as virtual basestation from the actual FloeNavi basestations that have a known + position on the the FloeNavi grid. + + :param dataset_time_bounds: (datetime list) time range of the dataset + :param basestation_setup: (FloeNaviBaseStationSetup instance) + :return: + """ + + # Load the collection of FloeNavi basestation files + logger.info("--- Establish FloeNavi grid ---") + floenavi_grid = FloeNaviGrid(self.cfg, basestation_setup) + floenavi_grid.set_time_bounds(dataset_time_bounds) + + # Get Reference station data from two positions + pos_ref, pos_hdg = floenavi_grid.get_virtual_reference_station_positions() + refstat_ais = GeoReferenceStation.from_two_basestations(pos_ref, pos_hdg) + return refstat_ais + + +class FloeNaviConfig(object): + """ + The configuration of the FloeNavi post processing system, including local paths + to data and the list of FloeNavi setups (list of basestations and their reference + positions) + """ + + def __init__(self, yaml_config_file=None, local_path_file=None): + """ + Loads the FloeNavi configuration (base station setup and local path) + :param yaml_config_file: (filepath): Optional: link to floenavi config file (base station lists, + floe navi setups, ...). Only required if a special config file needs to be used. If None, + the module will use the default config file in the root dir of the floenavi package + :param local_path_file: (filepath): Optional: link to floenavi input data folder. If None, + the module will use the default config file in the root dir of the floenavi package + """ + + # --- Save and digest arguments --- + + # Config file with list of base stations and the setups for different periods + if yaml_config_file is None: + self.yaml_config_file = DEFAULT_CONFIG_FILE + else: + self.yaml_config_file = yaml_config_file + logger.info("Custom floenavi config file: %s" % self.yaml_config_file) + if not os.path.isfile(self.yaml_config_file): + logger.error("Invalid floenavi config file: %s -> terminating" % str(self.yaml_config_file)) + sys.exit(1) + + if local_path_file is None: + self.local_path_file = DEFAULT_PATHDEF_FILE + else: + self.local_path_file = local_path_file + if not os.path.isfile(self.local_path_file): + logger.error("Invalid floenavi path definition file: %s -> terminating" % str(self.local_path_file)) + sys.exit(1) + + # Init properties + self._ref_ctlg = None + + # Init methods + self._parse_yaml() + self._create_valid_mmsi_catalog() + self._create_reference_data_catalog() + self._create_period_catalog() + + def get_station_ids(self, dataset_time_bounds): + """ + Searches for the mmsi number of the reference and heading station + :param dataset_time_bounds: (datetime list) The start and end time of the dataset + :return: + """ + + # Check time coverage for each FloeNavi configuration period + # (period section in the floenavi yaml config file) + mmsis = [] + for period_cfg in self.period_ctlg: + period_time_bounds = period_cfg.time_bounds + start_time_matches = dataset_time_bounds[0] >= period_time_bounds[0] + end_time_matches = dataset_time_bounds[0] <= period_time_bounds[1] + if start_time_matches and end_time_matches: + mmsis = period_cfg.mmsis + return mmsis + + def get_floenavi_cfg(self, dataset_time_bounds): + """ + Searches for the mmsi number of the reference and heading station + :param dataset_time_bounds: (datetime list) The start and end time of the dataset + :return: + """ + + # Check time coverage for each FloeNavi configuration period + # (period section in the floenavi yaml config file) + floenavi_cfg = None + for period_cfg in self.period_ctlg: + period_time_bounds = period_cfg.time_bounds + start_time_matches = dataset_time_bounds[0] >= period_time_bounds[0] + end_time_matches = dataset_time_bounds[0] <= period_time_bounds[1] + if start_time_matches and end_time_matches: + floenavi_cfg = period_cfg + + return floenavi_cfg + + def get_dship_id(self, mmsi): + """ + Return the dship id + :param mmsi: + :return: + """ + return self.mmsi_ctlg.get_basestation_dship_id(mmsi) + + def get_basestations_ctlg_entries(self, mmsis=None, period=None): + """ + Return a list of all catalog entries fitting optional search criteria + :param mmsis: (int) optional list of mmsis. If None, all mmsi's will be used + :param period: (list of datetime.date) optional period, If None, full period will be used + :return: list of FloeNaviBaseStation instances + """ + + # Use full mmsi range, if no search criteria for mmsi has been specified + if mmsis is None: + mmsis = self.mmsi_ctlg.base_mmsis + + # Use full mmsi range, if no search criteria for period has been specified + if period is None: + period = self.refstations_ctlg.full_date_range + + # Get the list of days from the time bounds of the data set + # The list of days is already in format contained in the filename + tdelta = period[1]-period[0] + n_days = tdelta.days + date = period[0] + date_strs = [date.strftime("%Y%m%d")] + for day_num in np.arange(n_days+1): + next_date = date + timedelta(days=int(day_num)) + date_strs.append(next_date.strftime("%Y%m%d")) + logger.debug("Period requires base station data from %g day(s)" % len(date_strs)) + + # Return file(s) as a list as it could be more than one + stations = [] + for date_str in date_strs: + station_list = self.refstations_ctlg.get_station_entries(date_str, mmsi_filter=mmsis) + stations.extend(station_list) + + # All done, return + return stations + + def get_basestation_files(self, mmsis=None, period=None): + """ + Get the list of file path + :param mmsis: (int) optional list of mmsis. If None, all mmsi's will be used + :param period: (list of datetime.date) optional period, If None, full period will be used + :return: list of all filepaths + """ + + # Use full mmsi range, if no search criteria for mmsi has been specified + if mmsis is None: + mmsis = self.mmsi_ctlg.base_mmsis + + # Use full mmsi range, if no search criteria for period has been specified + if period is None: + period = self.refstations_ctlg.full_date_range + + # Get the list of days from the time bounds of the data set + # The list of days is already in format contained in the filename + tdelta = period[1]-period[0] + n_days = tdelta.days + date = period[0] + date_strs = [date.strftime("%Y%m%d")] + for day_num in np.arange(n_days): + next_date = date + timedelta(days=int(day_num)) + date_strs.append(next_date.strftime("%Y%m%d")) + logger.debug("Period requires base station data from %g day(s)" % len(date_strs)) + + # Return file(s) as a list as it could be more than one + station_filepaths = [] + for date_str in date_strs: + file_paths = self.refstations_ctlg.get_station_filepath(date_str, mmsi_filter=mmsis) + station_filepaths.extend(file_paths) + + # All done, return + return station_filepaths + + def _parse_yaml(self): + """ + Read the FloeNavi config file and map contents to a dictionary + :return: + """ + logger.debug("Read FloeNavi config file: %s" % os.path.abspath(self.yaml_config_file)) + with open(self.yaml_config_file) as fh: + self._cfg = yaml.load(fh, Loader=yaml.FullLoader) + + logger.debug("Read FloeNavi path definition file: %s" % os.path.abspath(self.local_path_file)) + with open(self.local_path_file) as fh: + self._pathdef = yaml.load(fh, Loader=yaml.FullLoader) + + def _create_valid_mmsi_catalog(self): + """ + Get the valid mmsi's for base stations and mobile stations + :return: + """ + self._mmsi_ctlg = FloeNaviMMSICtlg(self._cfg["mmsis"]) + + def _create_reference_data_catalog(self): + """ + Create a catalog of data files (for each mmsi and day) + :return: + """ + + # Find reference station files + filename_lookup = self._pathdef["filename_template"].format(mmsi="*", yyyymmdd="*") + filepath_lookup = os.path.join(self._pathdef["lookup_dir"], filename_lookup) + refstat_files = sorted(glob.glob(filepath_lookup)) + msg = "Found %g %s files in %s" + msg = msg % (len(refstat_files), self._pathdef["data_source_id"], self._pathdef["lookup_dir"]) + logger.info(msg) + + # Add file to catalog + self._ref_ctlg = RefDataFileCatalog(refstat_files, self._pathdef["filename_template"]) + + def _create_period_catalog(self): + """ + Create a catalog of periods with defined FloeNavi configurations + (mmsi of deployed basestations and their location on the FloeNavi grid + :return: + """ + self._prd_ctlg = [] + for period_cfg in self._cfg["periods"]: + self._prd_ctlg.append(FloeNaviBaseStationSetup(period_cfg)) + + @property + def mmsi_ctlg(self): + return self._mmsi_ctlg + + @property + def period_ctlg(self): + return self._prd_ctlg + + @property + def refstations_ctlg(self): + return self._ref_ctlg + + @property + def basestation_ids(self): + return self.mmsi_ctlg.base_mmsis + + @property + def mobiletation_ids(self): + return self.mmsi_ctlg.mobile_mmsis + + @property + def pathdef(self): + return dict(self._pathdef) + + @property + def date_coverage_start(self): + return self.refstations_ctlg.full_date_range[0].date() + + @property + def date_coverage_end(self): + return self.refstations_ctlg.full_date_range[1].date() + + +class FloeNaviMMSICtlg(object): + """ + Container for the list of valid MMSI for basestations and mobile stations + """ + + def __init__(self, cfg_dict): + """ + + :param cfg_dict: + """ + + # save args + self._cfg_dict = cfg_dict + + # Get dicts for base stations and mobile stations: {mmsi: labels} + self._basestations = {} + for item in self._cfg_dict["basestations"]: + self._basestations[item["mmsi"]] = BaseStationDef(item) + + self._mobilestations = {} + for item in self._cfg_dict["mobilestations"]: + self._mobilestations[item["mmsi"]] = BaseStationDef(item) + + def get_type(self, mmsi): + """ + Return the type (base or mobile) of a given mmsi, None if invalid mmsi + :param mmsi: + :return: + """ + if mmsi in self.base_mmsis: + station_type = "base" + elif mmsi in self.mobile_mmsis: + station_type = "mobile" + else: + station_type = None + return station_type + + def get_station_name(self, mmsi): + """ + Return the station name for a given mmsi, None if invalid mmsi + :param mmsi: + :return: + """ + if mmsi in self.base_mmsis: + station_name = self._basestations[mmsi].label + elif mmsi in self.mobile_mmsis: + station_name = self._mobilestations[mmsi].label + else: + station_name = None + return station_name + + def get_basestation_dship_id(self, mmsi): + """ + Return the DSHIP id ( + :param mmsi: + :return: + """ + # DShips ids only available for base stations + if mmsi not in self.base_mmsis: + logger.error("MMSI %d not a registered basestation -> no DSHIP id exists") + return None + return self._basestations[mmsi].dship_id + + @property + def base_mmsis(self): + return sorted(self._basestations.keys()) + + @property + def mobile_mmsis(self): + return sorted(self._mobilestations.keys()) + + +class BaseStationDef(object): + """ + Small data class for the mmsi station definition in the config file + Attributes: + - mmsi + - label + - dship_id (optional) + """ + + def __init__(self, cfg_dict): + """ + Create a basestation definition data class + :param cfg_dict: + """ + self.mmsi = cfg_dict["mmsi"] + self.label = cfg_dict["label"] + self.dship_id = cfg_dict.get("dship_id", None) + + +class FloeNaviBaseStationSetup(object): + """ + The configuration of the basestations (mmsi and their positions) + in the ice coordinate system (FloeNavi Grid) + """ + + def __init__(self, period_cfg_dict): + """ + + :param period_cfg_dict: + """ + self.period_cfg_dict = period_cfg_dict + + def get_grid_location(self, mmsi): + """ + Get the reference position of a basestation within the + FloeNavi grid + :param mmsi: + :return: + """ + for cfg in self.period_cfg_dict["basestations"]: + if int(cfg["mmsi"]) == mmsi: + return cfg["pos"] + return None, None + + @property + def time_bounds(self): + tcs = self.period_cfg_dict["valid_since"] + tce = self.period_cfg_dict["valid_until"] + if tce is None: + tce = datetime.now() + return [tcs, tce] + + @property + def mmsis(self): + return [entry["mmsi"] for entry in self.period_cfg_dict["basestations"]] + + @property + def icecs_trns_kwargs(self): + return self.period_cfg_dict["icecs_transformation"] + + +class RefDataFileCatalog(object): + """ + Collection of metadata for all files in the basestation data folder + """ + + def __init__(self, filepaths, parser): + + self.filepaths = filepaths + self._catalog = [] + for filepath in self.filepaths: + try: + catalog_item = RefDataCatalogItem.from_filepath(filepath, parser) + except ValueError: + logger.error("Unable to parse file: %s -> skipping file" % filepath) + continue + self._catalog.append(catalog_item) + + # Check results + # logger.info("MMSI's in catalog:") + # for mmsi, start_time, stop_time in self.mssi_coverage: + # msg = "mmsi: %d [%s - %s]" + # msg = msg % (mmsi, start_time.isoformat(), stop_time.isoformat()) + # logger.info(msg) + + def get_items(self, mmsi): + return [item for item in self._catalog if item.mmsi == mmsi] + + def get_station_entries(self, yyyymmdd, mmsi_filter=None): + """ + Get filepaths for files for a given date and mmsi (potentially more than one) + :param yyyymmdd: + :param mmsi_filter: + :return: + """ + if mmsi_filter is None: + mmsi_filter = self.mmsis + if isinstance(mmsi_filter, int): + mmsi_filter = [mmsi_filter] + filepaths = [item for item in self._catalog if item.yyyymmdd == yyyymmdd and item.mmsi in mmsi_filter] + return filepaths + + def get_station_filepath(self, yyyymmdd, mmsi_filter=None): + """ + Get filepaths for files for a given date and mmsi (potentially more than one) + :param yyyymmdd: + :param mmsi_filter: + :return: + """ + if mmsi_filter is None: + mmsi_filter = self.mmsis + if isinstance(mmsi_filter, int): + mmsi_filter = [mmsi_filter] + filepaths = [item.filepath for item in self._catalog if item.yyyymmdd == yyyymmdd and item.mmsi in mmsi_filter] + return filepaths + + @property + def mmsis(self): + return set(sorted([item.mmsi for item in self._catalog])) + + @property + def mssi_coverage(self): + # TODO: Refactor to a list that also contains the gaps between the first and last day + mssi_coverage = [] + for mmsi in self.mmsis: + items = self.get_items(mmsi) + start_times = [item.time_bounds[0] for item in items] + stop_times = [item.time_bounds[1] for item in items] + mssi_coverage.append([mmsi, np.amin(start_times), np.amax(stop_times)]) + return mssi_coverage + + @property + def full_date_range(self): + start_times = [] + stop_times = [] + for mmsi in self.mmsis: + items = self.get_items(mmsi) + start_times.extend([item.time_bounds[0] for item in items]) + stop_times.extend([item.time_bounds[1] for item in items]) + return [np.amin(start_times), np.amax(stop_times)] + + +class RefDataCatalogItem(object): + """ + Container for metadata for one AIS stations (basestation of mobilestation) file + """ + + def __init__(self, filepath, mmsi=None, yyyymmdd=None, time_coverage="daily"): + self.filepath = filepath + self.mmsi = int(mmsi) + self.time_coverage = time_coverage + self.yyyymmdd = yyyymmdd + self.ref_time = datetime.strptime(yyyymmdd, "%Y%m%d") + self.time_bounds = [ + datetime(self.ref_time.year, self.ref_time.month, self.ref_time.day, 0, 0, 0), + datetime(self.ref_time.year, self.ref_time.month, self.ref_time.day, 23, 59, 59)] + + @classmethod + def from_filepath(cls, filepath, parser): + + # Get parameter from the filename + filename = os.path.split(filepath)[-1] + result = parse.parse(parser, filename) + attrs = dict() + if result: + for key in result.named.keys(): + attrs[key] = result.named[key] + return cls(filepath, **attrs) + + @property + def filename(self): + return os.path.split(self.filepath)[-1] + + @property + def date(self): + return self.ref_time.date() + + +class FloeNaviGrid(object): + """ + A container for all basestation data in a defined period + and functionality to compute two virtual basestation + that define heading and origin for the ice coordindate system + based on two or more FloeNavi reference stations + """ + + def __init__(self, cfg, setup): + """ + + :param cfg: + :param setup: + """ + + # Store the variables + self.cfg = cfg + self.setup = setup + + # Attributes + self.time_bounds = None + self.basestations = None + self.basestation_dict = {} + + def set_time_bounds(self, time_bounds): + """ + Specify the time bounds of the data set, which will trigger parsing the appropriate + basestation data. + :param time_bounds: + :return: + """ + + # Store the current time_bounds + self.time_bounds = time_bounds + + # TODO: Narrow list of mmsi to one specified in config file (not all available in repo) + + # Get the list of base station files (yet unsorted) + basestations_ctlg_entries = self.cfg.get_basestations_ctlg_entries( + mmsis=self.setup.mmsis, period=time_bounds) + + # Load the base station data and append data for specific mmsi's + # This step is there to group data by mmsi in case there more than + # one file per mmsi + self.basestation_dict = {} + for basestation in basestations_ctlg_entries: + mmsi = basestation.mmsi + logger.info("Read Reference station file(s): %s" % basestation.filepath) + ais_data = DshipAISExtract(basestation.filepath) + pos = GeoPositionData(ais_data.time, ais_data.longitude, ais_data.latitude) + if mmsi not in self.basestation_dict: + self.basestation_dict[mmsi] = pos + else: + self.basestation_dict[mmsi].append(pos) + + def get_virtual_reference_station_positions(self, time_resolution_sec=10.0, **refstat_kwargs): + """ + Compute the position objects of two virtual reference stations, based on the number of available + basestations. This is primary input for the icedrift correction in the icedrift module. + :param time_resolution_sec: (float) The temporal resolution of the virtual base station positions + :param refstat_kwargs: (dict) keywords for the reference station computation + :return: + """ + + # Compute common timestamp + logger.info("Compute timestamp of virtual base stations [resolution: %.2f sec]" % time_resolution_sec) + duration_seconds = (self.time_bounds[1]-self.time_bounds[0]).seconds + seconds_increments = np.arange(0, duration_seconds+time_resolution_sec, time_resolution_sec) + time = np.array([self.time_bounds[0]+timedelta(seconds=sec) for sec in seconds_increments]) + + # Create a collection of all base stations + # NOTE: Making sure all base stations have the same temporal resolutions + logger.info("Create base station collection with harmonized timestamp") + self.basestations = FloeNaviBaseStationCollection() + for mmsi in self.basestation_dict: + basestation_pos = self.basestation_dict[mmsi] + basestation_pos.resample_to_time(time) + ref_pos = self.setup.get_grid_location(mmsi) + self.basestations.add_basestation(FloeNaviBaseStation(mmsi, ref_pos, basestation_pos)) + + # Pass information to grid solver and compute origin and heading from + # all base stations + logger.info("Compute FloeNavi grid solution from basestation collection") + fngrid = FloeNaviGridSolver(self.basestations) + fngrid.solve() + + # Get the virtual origin and heading position data objects + logger.info("Compute origin and heading station positions") + orig, hdg = fngrid.get_reference_stations(**refstat_kwargs) + return orig, hdg + + +class FloeNaviGridSolver(object): + """ + Class for managing and obtaining the origin of the FloeNavi grid in geographical + coordinates and heading of the x-axis from an arbitrary set of deployed AIS + basestations throughout the floe. + + This algorithm + is implemented in the FloeNaviGridSolution class. + + This class will repeat the computation over all unique permutations of basestations + and provides a merged solution with optional smoothing and filtering. + + """ + + def __init__(self, basestations, minimum_distance_m=100., smoothing_window=None): + """ + Init the FloeNavi Grid solver + :param basestations: + :param minimum_distance_m: + :param smoothing_window: + """ + + # Input arguments + self.basestations = basestations + self.minimum_distance_m = minimum_distance_m + self.smoothing_window = smoothing_window + + # Attributes + self.solutions = [] + self.origin_longitude = None + self.origin_latitude = None + self.origin_uncertainty_m = None + self.axes_heading = None + + def solve(self): + """ + Compute the multi base station solution for the origin and axes heading of the + FloeNavi grid from all possible permutations of the base station combinations + """ + + # Get permutations + mmsi_permutations = self.mmsi_permutations + n_permutations = len(mmsi_permutations) + logger.info("Base stations combinations available for FloeNavi grid solution: %d" % n_permutations) + + # Input validation + if n_permutations == 0: + logger.error("To few basestations available -> Terminating") + sys.exit(1) + + # Loop over all permutations + for mmsi1, mmsi2 in self.mmsi_permutations: + + # Get base station data objects + logger.info("Using combination: [%d -> %d]" % (mmsi1, mmsi2)) + bs1 = self.basestations.get_by_mmsi(mmsi1) + bs2 = self.basestations.get_by_mmsi(mmsi2) + + # Check distance between the two base stations + bs_distance_m = self.basestations_grid_distance(bs1, bs2) + if bs_distance_m >= self.minimum_distance_m: + msg = "Distance between stations: %.0fm [min: %.0fm -> baseline ok]" + msg = msg % (bs_distance_m, self.minimum_distance_m) + logger.info(msg) + if bs_distance_m < self.minimum_distance_m: + msg = "Distance between stations: %.0fm [min: %.0fm -> baseline to short]" + msg = msg % (bs_distance_m, self.minimum_distance_m) + logger.warning(msg) + continue + + # Compute the solutions and add to the list + self.solutions.append(FloeNaviGridSolution(bs1, bs2)) + + # --- Check result --- + + # No solution -> error + if self.n_solutions == 0: + logger.error("No FloeNavi grid solution available -> terminating") + sys.exit(1) + + # One solution -> return as is + if self.n_solutions == 1: + logger.info("Single solution, skipping merging") + solution = self.solutions[0] + self.origin_longitude = solution.origin_longitude + self.origin_latitude = solution.origin_latitude + self.origin_uncertainty_m = solution.origin_uncertainty_m + self.axes_heading = solution.axes_heading + return + + # Multiple solutions -> merge + logger.info("Merge %d FloeNavi grid solutions" % self.n_solutions) + # NOTE: Currently, this is a simple point-wise mean of the parameters + # TODO: Add quality check (e.g. the min-max distance) + parameter_names = ["origin_longitude", "origin_latitude", "axes_heading"] + parameter_units = dict(origin_longitude="deg", origin_latitude="deg", origin_uncertainty_m="m", + axes_heading="deg") + n_records = self.solutions[0].n_records + for parameter_name in parameter_names: + + # Merge solution for parameters in to a multi-dimensional array + # -> easier to calculate mean and statistics with numpy + var_array = np.ndarray(shape=(n_records, self.n_solutions)) + for solution_num in np.arange(self.n_solutions): + var = getattr(self.solutions[solution_num], parameter_name) + if var is None: + var = np.full((n_records), np.nan) + var_array[:, solution_num] = var + + # Report statistics based on spread of solutions + unit = parameter_units[parameter_name] + min_vals = np.nanmin(var_array, axis=1) + max_vals = np.nanmax(var_array, axis=1) + spread = max_vals - min_vals + msg = "Spread of solutions for %s: max: %f%s, sdev:%f%s" + msg = msg % (parameter_name, np.nanmax(spread), unit, np.nanstd(spread), unit) + logger.debug(msg) + + # Assign the mean solution + mean_solution = np.nanmean(var_array, axis=1) + setattr(self, parameter_name, mean_solution) + + @staticmethod + def basestations_grid_distance(bs1, bs2): + return np.sqrt((bs2.xc_grid-bs1.xc_grid)**2.0 + (bs2.yc_grid-bs1.yc_grid)**2.0) + + def get_reference_stations(self, heading_station_distance_m=500): + """ + Returns two position object (one for origin and one for specifying the heading of the primary + axis) that can be used as input for the icedrift package. + :return: origin, hdg: Position objects for origin and heading of the FloeNavi x-axis + """ + + # Origin stations + origin = GeoPositionData(self.basestations.time, self.origin_longitude, self.origin_latitude) + + # Compute heading station positions + hdg = origin.clone() + az = self.axes_heading + r = heading_station_distance_m + g = pyproj.Geod(ellps='WGS84') + for i in np.arange(hdg.n_records): + hdg.longitude[i], hdg.latitude[i], _ = g.fwd(hdg.longitude[i], hdg.latitude[i], az[i], r) + + return origin, hdg + + @property + def mmsis(self): + return self.basestations.mmsis + + @property + def mmsi_permutations(self): + """ + Return all permutations of base stations that can be used to solve for grid origin + and location + :return: + """ + mmsi_permutations = [] + values = self.mmsis + all_values = list(self.mmsis) + for entry in all_values: + values.pop(values.index(entry)) + for other_entry in values: + mmsi_permutations.append((entry, other_entry)) + return mmsi_permutations + + @property + def n_solutions(self): + return len(self.solutions) + + +class FloeNaviGridSolution(object): + """ + Class containing the algorithm for computing origin and heading from FloeNavi + basestations at arbitrary (but known) positions in the FloeNavi Grid. + + The main idea of the algorithm is to compute the azimuth of the FloeNavi x-axis by + comparing the measured azimuth between two base stations against the defined azimuth + in the FloeNavi grid coordinate system. The origin of the FloeNavi grid in + geographical coordinates than can be computed for the known positions of the + base station in the FloeNavi grid and the heading of the x-axis. + """ + + def __init__(self, bs1, bs2, distance_deviation_warning_m=20.0): + """ + Compute origin and heading of the FloeNavi grid based on the geographical + position of two base station with known positions in the grid + :param bs1: Basestation 1 (floenavi.FloeNaviBaseStation) + :param bs2: Basestation 2 (floenavi.FloeNaviBaseStation) + :param distance_deviation_warning_m: (float) Warning threshold for discrepancies between observed + and specified base station distances + """ + + # Store arguments + self.bs1 = bs1 + self.bs2 = bs2 + self.distance_deviation_warning_m = distance_deviation_warning_m + + # Init attributes + self.azimuth_stations_geo = None + self.azimuth_stations_grid = None + self.origin_longitude = None + self.origin_latitude = None + self.origin_uncertainty_m = None + self.axes_heading = None + + # Compute origin and heading + self._compute() + + def _compute(self): + """ + Get origin in geographical coordinates and heading of x-axis of the + FloeNavi grid + :return: + """ + + # TODO: Break down in sub-methods + + # Init the stereographic projection and compute geographical x,y + # positions of base stations + logger.info("Compute projection based on position of %d" % self.bs1.mmsi) + proj_dict = self.proj_dict + logger.debug(str(proj_dict)) + for bs in [self.bs1, self.bs2]: + bs.set_projection(proj_dict) + + # Compute average distance between base stations and compare to specified value + bs_true_distance = np.sqrt((self.bs2.geo.xc-self.bs1.geo.xc)**2. + (self.bs2.geo.yc-self.bs1.geo.yc)**2.) + bs_grid_distance = np.sqrt((self.bs2.xc_grid-self.bs1.xc_grid)**2. + (self.bs2.yc_grid-self.bs1.yc_grid)**2.) + + # Check for NaN's + invalid = np.where(bs_true_distance > 1.0e10)[0] + bs_true_distance[invalid] = np.nan + + average_true_distance = np.nanmean(bs_true_distance) + msg = "Observed base station distance: %.1f m [specified: %.1f m]" + logger.debug(msg % (average_true_distance, bs_grid_distance)) + average_distance_deviation = bs_grid_distance - average_true_distance + if average_distance_deviation > self.distance_deviation_warning_m: + msg = "High discrepancy between observed and specified base station distance: %.1f m" + logger.warning(msg % average_distance_deviation) + + # Compute the geographical azimuth between the two base stations + self.azimuth_stations_geo = self.angle_between_points(self.bs1.geo.xc, self.bs1.geo.yc, + self.bs2.geo.xc, self.bs2.geo.yc) + msg = "Geographical azimuth of base stations: %.1f deg [min: %.1f, max: %.1f]" + msg = msg % (np.nanmean(self.azimuth_stations_geo), + np.nanmin(self.azimuth_stations_geo), + np.nanmax(self.azimuth_stations_geo)) + logger.debug(msg) + + # Compute the reference azimuth of the two base stations in FloeNavi grid coordinates + self.azimuth_stations_grid = self.angle_between_points(self.bs1.xc_grid, self.bs1.yc_grid, + self.bs2.xc_grid, self.bs2.yc_grid) + logger.debug("Target azimuth of base stations: %.1f deg" % self.azimuth_stations_grid) + + # Reference axis in the FloeNavi grid is the x-axis + # self._azimuth_stations_grid = np.mod(self._azimuth_stations_grid - 90.0, 360.0) + + # Compute the heading of the x-axis + # TODO: This is currently hard-coded to the x-axis of the FloeNavi Grid, evaluate customization + self.rotation = 360. - (self.azimuth_stations_geo - self.azimuth_stations_grid) + + # Reference axes for FloeNavi grid is the x-axis + logger.info("Compute heading of FloeNavi x-axis") + self.axes_heading = self.rotation + 90.0 + self.axes_heading = self.axes_heading % 360.0 + + # Compute the origin in projection coordinates + # NOTE: This can be done for each base station + origins_xc = np.full((2, self.bs1.n_records), np.nan) + origins_yc = origins_xc.copy() + for i, bs in enumerate([self.bs1, self.bs2]): + xc, yc = self.get_origin(bs.geo.xc, bs.geo.yc, bs.range_m, bs.ref_angle, self.rotation) + invalid = np.where(np.logical_or(xc > 1.0e10, yc > 1.0e10))[0] + xc[invalid], yc[invalid] = np.nan, np.nan + origins_xc[i, :] = xc + origins_yc[i, :] = yc + + # Compute the average deviation (uncertainty) of the origin + origin_xc_delta = origins_xc[1]-origins_xc[0] + origin_yc_delta = origins_yc[1]-origins_yc[0] + self._origin_uncertainty_m = np.sqrt(origin_xc_delta**2.0 + origin_yc_delta**2.0) + invalid = np.where(self._origin_uncertainty_m > 1.0e10)[0] + self._origin_uncertainty_m[invalid] = np.nan + logger.debug("Average origin uncertainty: %.1fm" % np.nanmean(self._origin_uncertainty_m)) + + # Compute the origin in geographical coordinates + logger.info("Compute FloeNavi origin in lon/lat") + p = pyproj.Proj(**self.proj_dict) + origin_xc = np.nanmean(origins_xc, axis=0) + origin_yc = np.nanmean(origins_yc, axis=0) + self.origin_longitude, self.origin_latitude = p(origin_xc, origin_yc, inverse=True) + + @staticmethod + def angle_between_points(xc1, yc1, xc2, yc2, deg_output=True): + """ + Compute the angle in cartesian coordinates between two points + :param p1: + :param p2: + :param deg_output: + :return: + """ + + # Angle with respect to x-axis + angle = np.arctan2(yc2 - yc1, xc2 - xc1) + angle = np.mod(angle, 2.0*np.pi) + + # Convert output (optional) and return + if deg_output: + angle = np.rad2deg(angle) + return angle + + @staticmethod + def get_origin(xc, yc, range_m, ref_heading, rotation): + """ + Compute the origin from basestation data and the known rotation of the + ice reference system + :param xc: + :param yc: + :param range_m: + :param ref_heading: + :param rotation: + :return: + """ + direction_to_origin_deg = 180. - (rotation - ref_heading) + direction_to_origin_rad = np.deg2rad(direction_to_origin_deg) + vector_x = range_m * np.cos(direction_to_origin_rad) + vector_y = range_m * np.sin(direction_to_origin_rad) + return xc + vector_x, yc + vector_y + + @property + def proj_dict(self): + lon_0 = np.nanmedian(self.bs1.geo.longitude) + lat_0 = np.nanmedian(self.bs1.geo.latitude) + proj_dict = dict(proj='stere', lon_0=lon_0, lat_0=lat_0, lat_ts=lat_0, ellps='WGS84') + return proj_dict + + @property + def n_records(self): + return len(self.origin_longitude) + + +class FloeNaviBaseStationCollection(object): + """ + A container for all basestation data and their metadata in + a specific FloeNavi configuration + """ + + def __init__(self, basestations=None): + """ + Init the list of basestations + """ + if basestations is None: + self.list = [] + else: + self.list = basestations + + def add_basestation(self, basestation): + """ + Add a basestation to the collection + :param basestation: floenavi.FloeNaviBaseStation instance + :return: + """ + # TODO: input validation + self.list.append(basestation) + + def get_by_mmsi(self, mmsi, error_on_invalid=False): + """ + Return basestation data based on the mmsi + :param mmsi: + :param error_on_invalid: + :return: + """ + basestations = [basestation for basestation in self.list if basestation.mmsi == mmsi] + if len(basestations) == 0: + basestations = None + elif len(basestations) == 1: + basestations = basestations[0] + else: + if error_on_invalid: + raise ValueError("No basestation with mmsi: %d" % mmsi) + return basestations + + @property + def mmsis(self): + return list(sorted([basestation.mmsi for basestation in self.list])) + + @property + def time(self): + if self.n_basestations == 0: + return None + basestation = self.list[0] + return basestation.geo.time + + @property + def n_basestations(self): + return len(self.list) + + +class FloeNaviBaseStation(object): + """ + Container for FloeNavi basestation data including observed + geographic locations and known location on the FloeNavi grid + """ + + def __init__(self, mmsi, ref_pos, geo): + """ + Init the basestation data + :param mmsi: + :param ref_pos: + :param geo: + """ + + # Store Attributes + self.mmsi = mmsi + self.ref_pos = ref_pos + self.geo = geo + + def set_projection(self, proj_dict): + """ + Compute geographical cartesian coordinates with a given projection + (Usually polar-stereographic) + :param proj_dict: + :return: + """ + p = pyproj.Proj(**proj_dict) + self.geo.xc, self.geo.yc = p(self.geo.longitude, self.geo.latitude) + + @property + def xc_grid(self): + return self.ref_pos[0] + + @property + def yc_grid(self): + return self.ref_pos[1] + + @property + def range_m(self): + return np.sqrt(self.ref_pos[0]**2. + self.ref_pos[1]**2.) + + @property + def ref_angle(self): + return np.rad2deg(np.arctan2(self.ref_pos[1], self.ref_pos[0])) + + @property + def n_records(self): + return self.geo.n_records + + diff --git a/floenavi/ais.py b/floenavi/ais.py new file mode 100644 index 0000000000000000000000000000000000000000..6a7171c7df31cb042eedd56940e1615ed7e30d45 --- /dev/null +++ b/floenavi/ais.py @@ -0,0 +1,57 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from dateutil.parser import parse as du_parse + + +class DshipAISExtract(object): + """ + Class to ingest standard DSHIP/realtime extraction AIS station data. + """ + + def __init__(self, filepath, zero_lonlat_to_nan=True): + self.filepath = filepath + self.zero_lonlat_to_nan = zero_lonlat_to_nan + self.time = None + self._parse() + + def _parse(self): + """ + Read time, lat, lon, heading from DSHIP extract + :return: + """ + + # Read ASCII extract + self.csv_var_dtypes = dict( + names=['date', 'time', 'latitude', 'longitude'], + formats=['U10', 'U8', 'f4', 'f4']) + self.content = np.loadtxt(self.filepath, dtype=self.csv_var_dtypes) + + # Compute datetime + self.time = np.ndarray(shape=tuple([self.n_records]), dtype=object) + for i in np.arange(self.n_records): + datetime_str = self.content["date"][i].replace("/", "")+"T"+self.content["time"][i] + self.time[i] = du_parse(datetime_str) + + if self.zero_lonlat_to_nan: + lon_is_zero = np.abs(self.longitude) <= 1.e-6 + lat_is_zero = np.abs(self.latitude) <= 1.e-6 + lonlat_are_zero = np.where(np.logical_and(lon_is_zero, lat_is_zero))[0] + for target in ["longitude", "latitude"]: + self.content[target][lonlat_are_zero] = np.nan + + @property + def longitude(self): + return self.content["longitude"] + + @property + def latitude(self): + return self.content["latitude"] + + @property + def n_records(self): + return self.content.size + + @property + def time_bounds(self): + return [self.time[0], self.time[-1]] \ No newline at end of file diff --git a/floenavi/dataset.py b/floenavi/dataset.py new file mode 100644 index 0000000000000000000000000000000000000000..dcca84e6f7a6f6812668dd7ed87a72b0b144a7d3 --- /dev/null +++ b/floenavi/dataset.py @@ -0,0 +1,315 @@ +# -*- coding: utf-8 -*- + +""" +This module contains dataset classes intended for the use with the FloeNavi().get_xy_coordinates() +method. + +Convention for dataset classes: + + Required attributes + - time: array of datetime.datetime in UTC (or at least in the same timezone as the base stations) + - longitude: float array + - latitude: float array + - xc: float array or None if not yet computed + - yc: float array or None if not yet computed + - n_records: Number of records + - time_bounds: [time_coverage_start, time_coverage_end] as datetime.datetime objects + + Required methods: + - update_floenavi_coords(self, xc, yc): Setting of (x,y) by FloeNavi().get_xy_coordinates() + +""" + + +import sys +from loguru import logger +import numpy as np +from matplotlib.dates import date2num +from dateutil.parser import parse as du_parse + + +class Waypoints(object): + """ + FloeNavi dataset for GPS waypoints data with label, time, longitude, latitude. To meant be used as + + wps = Waypoints() + for label, time, lon, lat in <some-iterable>: + add(label, time, lon, lat) + floenavi = FloeNavi() + floenavi.get_xy_coordinates(wps) + + After succesful transformation, the waypoint instance will have xc, yc attributes. + + The class can also be initialized from a csv file with the class method: + + wps = Waypoints.from_csv(filename) + + given that the input file has comma-separated columns for label, time, longitude, latitude: + + Example: + + 081-Polarstern, 2019-10-06T01:43:20, 133.780012987554073, 85.109757967293262 + 082-OceanCity-PowerHub, 2019-10-06T02:04:43, 133.748755026608706, 85.110806962475181 + 083-Ridge1, 2019-10-06T02:13:21, 133.747596982866526, 85.110585009679198 + + """ + + def __init__(self): + """ + Create a Waypoint instance. Waypoint data has to be added with the method `add(label, time, lon, lat)` + """ + self._wps = [] + + @classmethod + def from_csv(cls, filepath): + """ + Read a standard waypoint csv file with 4 columns: + + waypoint label (str), datetime (str), latitude (float), longitude (float) + + :param filepath: + :return: + """ + + # --- Read the csv content --- + names = ('label', 'time', 'latitude', 'longitude') + formats = ('U96', 'U96', 'f4', 'f4') + content = np.loadtxt(filepath, skiprows=1, delimiter=',', dtype={'names': names, 'formats': formats}) + + # Construct the instance and compute waypoints + wps = cls() + for label, time_str, lat, lon in content: + time = du_parse(time_str) + wps.add(label, time, lon, lat) + + # All done, return initialized instance + return wps + + def add(self, *args): + self._wps.append(WaypointData(*args)) + + def update_floenavi_coords(self, xc, yc): + """ + Add x, y coordinates to each waypoint + :param xc: + :param yc: + :return: + """ + for i, (xc_val, yc_val) in enumerate(zip(xc, yc)): + self.wps[i].add_xy(xc_val, yc_val) + + @property + def wps(self): + return list(self._wps) + + @property + def n_waypoints(self): + return len(self.wps) + + @property + def waypoints(self): + return self.wps + + @property + def time(self): + return np.array([wp.time for wp in self.wps]) + + @property + def latitude(self): + return np.array([wp.latitude for wp in self.wps]) + + @property + def longitude(self): + return np.array([wp.longitude for wp in self.wps]) + + @property + def labels(self): + return np.array([wp.label for wp in self.wps]) + + @property + def xc(self): + return np.array([wp.xc for wp in self.wps]) + + @property + def yc(self): + return np.array([wp.yc for wp in self.wps]) + + @property + def time_bounds(self): + return [np.amin(self.time), np.amax(self.time)] + + +class GPSTrack(object): + """ + FloeNavi dataset for track data with time, longitude, latitude. To meant be used as + + track = GPSTrack(times, longitude, latitude) + floenavi = FloeNavi() + floenavi.get_xy_coordinates(track) + + After succesful transformation, the track instance will have xc, yc attributes. + + The class can also be initialized from a csv file with the class method: + + wps = GPSTrack.from_csv(filename) + + given that the input file has comma-separated columns for time, longitude, latitude: + + Example: + + 2019-10-24T05:22:18.433593, 129.20388794, 85.39130402 + 2019-10-24T05:22:18.554687, 129.20387268, 85.39130402 + 2019-10-24T05:22:18.675781, 129.20387268, 85.39130402 + + """ + + def __init__(self, time, longitude, latitude, xc=None, yc=None): + """ + Create a GPS track data set for use with FloeNavi + :param time: (datetime array) timestamp in UTC + :param longitude: (float array) longitude values + :param latitude: (float array) latitude values + :param xc: (float array) Optional: FloeNavi x coordinates (will be set by FloeNavi) + :param yc: (float array) Optional: FloeNavi y coordinates (will be set by FloeNavi) + """ + + # --- Store input --- + self.time = time + self.longitude = longitude + self.latitude = latitude + self.xc = xc + self.yc = yc + + @classmethod + def from_csv(cls, filepath, zero_latlon_to_nan=True, interpolate_gaps=True): + """ + Read a standard waypoint csv file with 4 columns: + + waypoint label (str), datetime (str), latitude (float), longitude (float) + + :param filepath: + :param zero_latlon_to_nan: + :param interpolate_gaps: + :return: + """ + + # --- Read the csv content --- + names = ('time', 'latitude', 'longitude') + formats = ('U96', 'f4', 'f4') + content = np.loadtxt(filepath, skiprows=0, delimiter=',', dtype={'names': names, 'formats': formats}) + + # Construct the instance and compute waypoints + time = np.ndarray(shape=content.shape, dtype=object) + for i, (time_str, _, _) in enumerate(content): + time[i] = du_parse(time_str) + + # Sanity check and set invalid entries to NaN + # NOTE: Currently, it is assumed that invalid entries have a value of 0.0 + # TODO: Generalize this to allow diverse nodata flags + longitude, latitude = content["latitude"], content["longitude"] + is_invalid = np.abs(latitude) < 1 + if zero_latlon_to_nan: + invalid_indices = np.where(is_invalid)[0] + longitude[invalid_indices] = np.nan + latitude[invalid_indices] = np.nan + + # Simple gap interpolation + if interpolate_gaps: + datenum = date2num(time) + valid_indices = np.where(np.logical_not(is_invalid))[0] + longitude = np.interp(datenum, datenum[valid_indices], longitude[valid_indices]) + latitude = np.interp(datenum, datenum[valid_indices], latitude[valid_indices]) + + # All done, return initialized instance + track = cls(time, longitude, latitude) + return track + + def update_floenavi_coords(self, xc, yc): + """ + Add x, y coordinates to each waypoint + :param xc: + :param yc: + :return: + """ + + # Simple Input sanity check (shape of arrays must match) + for input_var in [xc, yc]: + if input_var.shape != self.time.shape: + msg = "Dimension mismatch: GPSTrack: %s != Coordinates: %s -> terminating" + logger.error(msg % (str(self.time.shape), str(input_var.shape))) + sys.exit(1) + + # Store (x,y) coordinates + self.xc = xc + self.yc = yc + + def export_csv(self, output_filepath): + """ + Export the file as csv including xc & yc + :param output_filepath: + :return: + """ + + # Check if data set has assigned x, y coordinates + # if not -> display warning and use warning + if self.has_xy_coordinates: + xc, yc = self.xc, self.yc + else: + xc, yc = np.full(self.time.shape, np.nan), np.full(self.time.shape, np.nan) + logger.warning("Exporting GPSTrack without xc, yc: method `update_floenavi_coords` has not been used") + + # Create the track output + fmt = "%s, %.8f, %.8f, %.2f, %.2f\n" + with open(output_filepath, "w") as fh: + for i in np.arange(self.n_records): + time_str = self.time[i].isoformat() + fh.write(fmt % (time_str, self.latitude[i], self.longitude[i], xc[i], yc[i])) + + @property + def n_records(self): + return len(self.time) + + @property + def time_bounds(self): + return [np.amin(self.time), np.amax(self.time)] + + @property + def has_xy_coordinates(self): + """ + Indicates if this instnace has assigned floenavi (x,y) coordinates + :return: flag (bool) + """ + return self.xc is not None and self.yc is not None + + +class WaypointData(object): + """ + Container for individual waypoints, meant for use with dataset.Waypoints + """ + + def __init__(self, label, time, longitude, latitude, xc=None, yc=None): + """ + Create a Waypoint entry + :param label: (str) label of the waypoints + :param time: (datetime.datetime) timestamp in UTC + :param longitude: (float) longitude in degrees + :param latitude: (float) latitude in degrees + :param xc: (float) optional: FloeNavi x-coordinate (will be overwritten by `add_xy` method) + :param yc: (float) optional: FloeNavi y-coordinate (will be overwritten by `add_xy` method) + """ + self.label = label + self.time = time + self.longitude = longitude + self.latitude = latitude + self.xc = xc + self.yc = yc + + def add_xy(self, xc, yc): + """ + Add the FloeNavi grid coordinates to the waypoint + :param xc: (float) FloeNavi x-coordinate + :param yc: (float) FloeNavi x-coordinate + :return: + """ + self.xc = xc + self.yc = yc diff --git a/floenavi_cfg.yaml b/floenavi_cfg.yaml new file mode 100644 index 0000000000000000000000000000000000000000..25895bafae137d1bab02dd0c0022692747163a31 --- /dev/null +++ b/floenavi_cfg.yaml @@ -0,0 +1,319 @@ +# General idea is that the setup of the Floe Navigation system +# will change with time. e.g. by changing the AIS base stations, +# adjusting headings etc. The configuration therefore has to be +# provided in certain periods, with all config data. A new period +# has to be added, if the configuration changes + + +# List of mmsi allocated to basestations and mobile stations +# as well as the labels according to the official +# frequency allocation +# +# NOTE: This list can be expected to be fixed and will +# not change during MOSAiC +mmsis: + + basestations: + - mmsi: 211003810 + label: TEST-MOSAIC-B01 + dship_id: !!null + - mmsi: 211003811 + label: TEST-MOSAIC-B02 + dship_id: !!null + - mmsi: 211003812 + label: TEST-MOSAIC-B03 + dship_id: AIS.FLOEN.03 + - mmsi: 211003813 + label: TEST-MOSAIC-B04 + dship_id: AIS.FLOEN.04 + - mmsi: 211003814 + label: TEST-MOSAIC-B05 + dship_id: AIS.FLOEN.05 + - mmsi: 211003820 + label: TEST-MOSAIC-B06 + dship_id: AIS.FLOEN.06 + - mmsi: 211003821 + label: TEST-MOSAIC-B07 + dship_id: AIS.FLOEN.07 + - mmsi: 211003822 + label: TEST-MOSAIC-B08 + dship_id: AIS.FLOEN.08 + - mmsi: 211003823 + label: TEST-MOSAIC-B09 + dship_id: AIS.FLOEN.09 + + mobilestations: + - mmsi: 211003815 + label: TEST-MOSAIC-M01 + - mmsi: 211003816 + label: TEST-MOSAIC-M02 + - mmsi: 211003817 + label: TEST-MOSAIC-M03 + - mmsi: 211003818 + label: TEST-MOSAIC-M04 + - mmsi: 211003819 + label: TEST-MOSAIC-M05 + - mmsi: 211003824 + label: TEST-MOSAIC-M06 + - mmsi: 211003825 + label: TEST-MOSAIC-M07 + - mmsi: 211003826 + label: TEST-MOSAIC-M08 + - mmsi: 211003827 + label: TEST-MOSAIC-M09 + - mmsi: 211003828 + label: TEST-MOSAIC-M10 + - mmsi: 211003829 + label: TEST-MOSAIC-M11 + - mmsi: 211003830 + label: TEST-MOSAIC-M12 + - mmsi: 211003831 + label: TEST-MOSAIC-M13 + - mmsi: 211003832 + label: TEST-MOSAIC-M14 + + +# FloeNavi periods as a list +periods: + + # NOTE: There were two experimental setup before the start of the operational phase. These are + # ignored here, because their position with respect to the later maintained base stations is + # difficult to measure (basestation were moved without overlap, while the ship itself was moving) + # + # First experimental setup: Two station with one on starboard and one to bow + # + # Seconds experimental setup: One station at the orginal power hub and one 500 m to the starboard + # side with bearing 090. The ship turned shortly afterwards, so that this location is now + # longer at a bearing 090. + + # ---- + + # Initial setup in final floenavi configuration consisting of the origin basestation (211003814) + # at the Polarstern power hub (moved location after breakup at the bow) and the Ocean city power + # hub station (211003821). + # NOTE: The short-lived base station at ROV city is ignored + - id: 03-org-oc + valid_since: 2019-10-16T00:00:00 + valid_until: 2019-10-20T23:59:59 + basestations: + - mmsi: 211003814 + pos: [0.0, 0.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # ---- + + # Adding the base stations at remote sensing (211003822) + # and MET city (211003822). The configuration ended + # with the failure of the origin base station (211003814) + - id: 04-org-oc-rs-met + valid_since: 2019-10-21T00:00:00 + valid_until: 2019-10-27T23:59:59 + basestations: + - mmsi: 211003814 + pos: [0.0, 0.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + - mmsi: 211003822 + pos: [263.0, 43.0] + - mmsi: 211003823 + pos: [449.0, 200.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # No base station at the origin (failure of 211003814). Grid has to + # be computed from station at ocean city, remote sensing at ROV city + - id: 05-oc-rs-met + valid_since: 2019-10-28T00:00:00 + valid_until: 2019-11-04T23:59:59 + basestations: + - mmsi: 211003821 + pos: [284.0, 0.0] + - mmsi: 211003822 + pos: [263.0, 43.0] + - mmsi: 211003823 + pos: [449.0, 200.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # 211003814 has been re-deployed at the approximated old (500, 0) site + # and 211003820 has been deployed at ROV City. + # NOTE: 2011003813 has also been active at the approximated origin station + # however it is not listed here, as there is some significant + # position scatter at certain periods + - id: 06-oc-rs-met-rov-500m + valid_since: 2019-11-05T00:00:00 + valid_until: 2019-11-12T23:59:59 + basestations: + - mmsi: 211003814 + pos: [449.0, 6.0] + - mmsi: 211003820 + pos: [134.0, 138.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + - mmsi: 211003822 + pos: [263.0, 43.0] + - mmsi: 211003823 + pos: [449.0, 200.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # 211003814 has been removed from the (500, 0) site since there was no power + # NOTE: 2011003813 has also been active at the approximated origin station + # however it is not listed here, as there is some significant + # position scatter at certain periods + - id: 07-oc-rs-met-rov + valid_since: 2019-11-13T00:00:00 + valid_until: 2019-11-15T23:59:59 + basestations: + - mmsi: 211003820 + pos: [134.0, 138.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + - mmsi: 211003822 + pos: [263.0, 43.0] + - mmsi: 211003823 + pos: [449.0, 200.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Period of strong movement in the ice (mid-November shear event) + # Only the origin station (211003813) and the station at Ocean City + # (211003821) can be considered to have stayed on their position + # in the floenavi grid + - id: 08-orig-oc + valid_since: 2019-11-16T00:00:00 + valid_until: 2019-11-19T23:59:59 + basestations: + - mmsi: 211003813 + pos: [22.0, -17.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Previous setup supported by an additional battery-powered basestation + # at the unfinished watchtower tower. Basestation at the other side of the + # shear zone are not considered stable at this point and ignored in this + # period + - id: 09-orig-oc-watchtower + valid_since: 2019-11-20T00:00:00 + valid_until: 2019-11-22T23:59:59 + basestations: + - mmsi: 211003813 + pos: [22.0, -17.0] + - mmsi: 211003821 + pos: [284.0, 0.0] + - mmsi: 211003814 + pos: [181.0, -234.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Basestation at Ocean City (211003821) has been moved + - id: 10-orig-watchtower + valid_since: 2019-11-23T00:00:00 + valid_until: 2019-11-24T23:59:59 + basestations: + - mmsi: 211003813 + pos: [22.0, -17.0] + - mmsi: 211003814 + pos: [181.0, -234.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Ocean city base station (211003821) is active + # again at new location + - id: 11-orig-oc-watchtower + valid_since: 2019-11-25T00:00:00 + valid_until: 2019-11-28T23:59:59 + basestations: + - mmsi: 211003813 + pos: [22.0, -17.0] + - mmsi: 211003814 + pos: [181.0, -234.0] + - mmsi: 211003821 + pos: [293.0, -44.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Origin base station has been removed ánd the grid is held + # by Ocean City (211003821) and the watch tower station + # (211003814) + - id: 12-oc-watchtower + valid_since: 2019-11-29T00:00:00 + valid_until: 2019-11-30T23:59:59 + basestations: + - mmsi: 211003814 + pos: [181.0, -234.0] + - mmsi: 211003821 + pos: [293.0, -44.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # New origin (again at another location) added (211003820) + - id: 13-orig-oc-watchtower + valid_since: 2019-12-01T00:00:00 + valid_until: 2019-12-03T23:59:59 + basestations: + - mmsi: 211003814 + pos: [181.0, -234.0] + - mmsi: 211003820 + pos: [63.0, -20.0] + - mmsi: 211003821 + pos: [293.0, -44.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # Watch tower station replaced (211003814 -> 211003813) + - id: 14-orig-oc-watchtower + valid_since: 2019-12-04T00:00:00 + valid_until: 2019-12-11T23:59:59 + basestations: + - mmsi: 211003813 + pos: [181.0, -234.0] + - mmsi: 211003820 + pos: [63.0, -20.0] + - mmsi: 211003821 + pos: [293.0, -44.0] + icecs_transformation: + heading_axis: "positive x" + heading_offset_deg: 0.0 + dx: 0.0 + dy: 0.0 + + # ---- \ No newline at end of file diff --git a/floenavi_local_path_def.yaml b/floenavi_local_path_def.yaml new file mode 100644 index 0000000000000000000000000000000000000000..2b5e684e33c7347948fcbcbf5ef4f3bb8331f613 --- /dev/null +++ b/floenavi_local_path_def.yaml @@ -0,0 +1,7 @@ +# This config file indicates the location of the AIS base station +# data to the floenavi package. + +# General Information +data_source_id: "ais" +lookup_dir: D:\data\floenavi\basestations +filename_template: floenavi-ais_{mmsi}-{yyyymmdd}.dat \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000000000000000000000000000000000000..efd661e1d03dbb187a838500105ed8bc2f25a2b9 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +loguru +matplotlib +numpy +icedrift \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..4992f49b6d45e04cca2fac5356509142fb3fa9a9 --- /dev/null +++ b/setup.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- + +from setuptools import setup, find_packages + +import os + +# Get the readme +with open('README.rst') as f: + readme = f.read() + +# Get the licence +with open('LICENSE') as f: + license_content = f.read() + +# Get the version +mypackage_root_dir = os.path.dirname(os.path.abspath(__file__)) +with open(os.path.join(mypackage_root_dir, "floenavi", 'VERSION')) as version_file: + version = version_file.read().strip() + +# Package requirements +with open("requirements.txt") as f: + requirements = f.read().splitlines() + + +setup( + name='floenavi', + version=version, + description='python toolbox for post-processing FloeNavi data (requires icedrift)', + long_description=readme, + author='Stefan Hendricks', + author_email='stefan.hendricks@awi.de', + scripts=['bin/floenavi-set-cfgdir.py', + 'bin/floenavi_gpstrack2xy.py'], + url='', + license=license_content, + install_requires=requirements, + packages=find_packages(exclude=('tests', 'docs')), + include_package_data=True +)