diff --git a/scripts/cmems_obs/fix_obs_tracks.py b/scripts/cmems_obs/fix_obs_tracks.py new file mode 100644 index 000000000..2fd4a28e5 --- /dev/null +++ b/scripts/cmems_obs/fix_obs_tracks.py @@ -0,0 +1,79 @@ +import argparse + +import numpy as np +import pandas as pd +from sqlalchemy import create_engine +from sqlalchemy.dialects.mysql import insert + + +def mysql_upsert(table, conn, keys, data_iter): + data = [dict(zip(keys, row)) for row in data_iter] + stmt = insert(table.table).values(data) + update_stmt = stmt.on_duplicate_key_update(**{c.name: c for c in stmt.inserted}) + conn.execute(update_stmt) + + +def main(uri: str, platform_ids: str | list | None = None): + """ + Corrects station coordinates at crossing of date line so that tracks are displayed + continuously in the Navigator. Arguments should be passed on commandline. e.g: + python scripts/data_importers/get_cmems_month.py "uri" -p 12345 45678 + + :param uri: The URI string of the MariaDB Observation database + :param platform_ids: IDs of platforms to be corrected. If not provided all + platforms will be updated. + """ + + engine = create_engine(uri) + + if isinstance(platform_ids, str): + platform_ids = [platform_ids] + elif platform_ids is None: + platforms = pd.read_sql("SELECT * FROM platforms;", con=engine) + platform_ids = platforms.id.values + + for platform_id in platform_ids: + stations = pd.read_sql( + f"SELECT * FROM stations where platform_id={platform_id};", + con=engine, + index_col="id", + ).sort_values(by="time") + + corrected_lon = np.copy(stations.longitude.values) + diffs = np.diff(corrected_lon) + crossings = np.where(np.abs(diffs) > 180)[0] + + if len(crossings) > 0: + print(f"Updating track of platform {platform_id}.") + + for crossing in crossings: + if diffs[crossing] > 0: + corrected_lon[crossing + 1 :] -= 360 + else: + corrected_lon[crossing + 1 :] += 360 + stations.longitude = corrected_lon + stations.to_sql( + "stations", engine, if_exists="append", method=mysql_upsert, index=False + ) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + prog="Initialize/append CMEMS Observation data", + description=( + "Add monthly CMEMS Global Ocean In-Situ Near-Real-Time " + "Observations data to ONAV Obs database." + ), + ) + parser.add_argument("uri", type=str, help="URI of the observation database.") + parser.add_argument( + "-p", + "--platform_ids", + type=str, + help="ID of platforms to modify.", + default=None, + ) + + args = parser.parse_args() + + main(args.uri, args.platform_ids) diff --git a/scripts/cmems_obs/get_cmems_obs_day.py b/scripts/cmems_obs/get_cmems_obs_day.py new file mode 100644 index 000000000..70f4b3f4d --- /dev/null +++ b/scripts/cmems_obs/get_cmems_obs_day.py @@ -0,0 +1,88 @@ +import argparse +import os +from datetime import date, timedelta + +import copernicusmarine + +from import_cmems_obs import import_cmems_obs_mp + + +def main(uri: str, output_dir: str, date_str: str): + """ + Downloads data from CMEMS latest observation dataset matching the provided date + string to add to the observation database. + Deletes dowloaded files after data is added. Arguments should be passed on + commandline. e.g: + python scripts/data_importers/init_cmems_obs.py "uri" "/output_dir/" -d 20250101 + + :param uri: The URI string of the MariaDB Observation database + :param output_dir: The directory that you want to save the data to. + :param date_str: The date in YYYYMMDD format. + """ + + prod_id = "INSITU_GLO_PHYBGCWAV_DISCRETE_MYNRT_013_030" + meta = copernicusmarine.describe(product_id=prod_id) + version = [v.label for v in meta.products[0].datasets[0].versions][-1] + + obs_types = ["GL", "DB", "PF", "CT"] + + for obs in obs_types: + print(obs) + obs_filter = f"*/{date_str}/*_*_{obs}_*_*.nc" + resp = copernicusmarine.get( + dataset_id="cmems_obs-ins_glo_phybgcwav_mynrt_na_irr", + dataset_version=version, + output_directory=output_dir, + filter=obs_filter, + dataset_part="latest", + sync=True, + ) + file_list = [f.file_path for f in resp.files] + + if len(file_list) == 0: + continue + + try: + import_cmems_obs_mp(uri, file_list, obs) + + for file in file_list: + os.remove(file) + except Exception as e: + print(e) + with open("obs_error.log", "a") as f: + f.write(f"{obs_filter}\n") + f.write("\n") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + prog="Initialize/append CMEMS Observation data", + description=( + "Add monthly CMEMS Global Ocean In-Situ Near-Real-Time " + "Observations data to ONAV Obs database." + ), + ) + parser.add_argument("uri", type=str, help="URI of the observation database.") + parser.add_argument( + "output_dir", + type=str, + help="Output directory for observation data.", + ) + parser.add_argument( + "-d", + "--date", + type=str, + help=( + "The date to add in YYYYMMDD format. If omitted then yesterday's date " + "will be used." + ), + default=None, + ) + + args = parser.parse_args() + + if args.date: + main(args.uri, args.output_dir, args.date) + else: + date_str = (date.today() - timedelta(days=1)).strftime("%Y%m%d") + main(args.uri, args.output_dir, date_str) diff --git a/scripts/cmems_obs/get_cmems_obs_month.py b/scripts/cmems_obs/get_cmems_obs_month.py new file mode 100644 index 000000000..5e46f8c7d --- /dev/null +++ b/scripts/cmems_obs/get_cmems_obs_month.py @@ -0,0 +1,90 @@ +import argparse +import os +from datetime import date, datetime, timedelta + +import copernicusmarine + +from import_cmems_obs import import_cmems_obs_mp + + +def main(uri: str, output_dir: str, month_key: str): + """ + Downloads one month of CMEMS monthly observation data to append to the observation + database. Deletes dowloaded files after data is added. Arguments should be passed + on commandline. e.g: + python scripts/data_importers/get_cmems_month.py "uri" "/output_dir/" -m "month_key" + + :param uri: The URI string of the MariaDB Observation database + :param output_dir: The directory that you want to save the data to. + :param month_key: The month in YYYYMM format + """ + + prod_id = "INSITU_GLO_PHYBGCWAV_DISCRETE_MYNRT_013_030" + meta = copernicusmarine.describe(product_id=prod_id) + version = [v.label for v in meta.products[0].datasets[0].versions][-1] + + obs_types = ["GL", "DB", "PF", "CT"] + + for obs in obs_types: + obs_filter = f"*/{obs}/{month_key}/*.nc" + resp = copernicusmarine.get( + dataset_id="cmems_obs-ins_glo_phybgcwav_mynrt_na_irr", + dataset_version=version, + output_directory=output_dir, + filter=obs_filter, + dataset_part="monthly", + sync=True + ) + + file_list = [f.file_path for f in resp.files] + + try: + import_cmems_obs_mp(uri, file_list, obs) + + for file in file_list: + os.remove(file) + except Exception as e: + print(e) + with open("obs_error.log", "a") as f: + f.write(f"{obs_filter}\n") + f.write("\n") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + prog="Initialize/append CMEMS Observation data", + description=( + "Add monthly CMEMS Global Ocean In-Situ Near-Real-Time " + "Observations data to ONAV Obs database." + ), + ) + parser.add_argument("uri", type=str, help="URI of the observation database.") + parser.add_argument( + "output_dir", + type=str, + help="Output directory for observation data.", + ) + parser.add_argument( + "-m", + "--month", + type=str, + help=( + "The month to add in YYYYMM format. If omitted then all months from 2017 " + "onward will be downloaded." + ), + default=None, + ) + + args = parser.parse_args() + + if args.month: + main(args.uri, args.output_dir, args.month) + else: + start_month = date.today().replace(day=1) + end_month = datetime(2017, 1, 1) + while end_month >= end_month: + month_key = start_month.strftime("%Y%m") + + main(args.uri, args.output_dir, month_key) + + start_month = (start_month - timedelta(days=28)).replace(day=1) diff --git a/scripts/cmems_obs/import_cmems_obs.py b/scripts/cmems_obs/import_cmems_obs.py new file mode 100644 index 000000000..b0a2d2ea4 --- /dev/null +++ b/scripts/cmems_obs/import_cmems_obs.py @@ -0,0 +1,378 @@ +import argparse +import itertools +import multiprocessing +import os +import sys +from pathlib import PosixPath + +import gsw +import numpy as np +import pandas as pd +import xarray as xr +from sqlalchemy import create_engine, select, insert +from sqlalchemy.exc import IntegrityError +from sqlalchemy.orm import Session + +current = os.path.dirname(os.path.realpath(__file__)) +parent = os.path.dirname(os.path.dirname(current)) +sys.path.append(parent) + +from data.observational import DataType, Platform, Sample, Station + +""" +QC Flag meanings: + 0: "no_qc_performed", + 1: "good_data", + 2: "probably_good_data", + 3: "bad_data_that_are_potentially_correctable", + 4: "bad_data", + 5: "value_changed", + 6: "value_below_detection", + 7: "nominal_value", + 8: "interpolated_value", + 9: "missing_value" +Only use data with values 1, 2, 5, 7 or 8 +""" + + +def get_platform_variables(platform_type): + + variables = None + match platform_type: + case "GL": + variables = ["PRES", "PSAL", "TEMP", "CNDC"] + case "DB": + variables = ["TEMP", "ATMS"] + case "PF": + variables = ["TEMP", "PSAL"] + case "CT": + variables = ["PRES", "PSAL", "TEMP"] + + return variables + + +def insert_stations(pd_table, conn, keys, data_iter): + + data = [dict(zip(keys, d)) for d in data_iter] + stmt = insert(Station).values(data).prefix_with("IGNORE") + conn.execute(stmt) + conn.commit() + + +def insert_samples(pd_table, conn, keys, data_iter): + + data = [dict(zip(keys, d)) for d in data_iter] + stmt = insert(Sample).values(data).prefix_with("IGNORE") + conn.execute(stmt) + conn.commit() + + +def get_station_ids(row, session): + + stmt = select(Station).where( + Station.platform_id == row.platform_id, + Station.time == row.time, + Station.latitude == row.latitude, + Station.longitude == row.longitude, + ) + station_data = session.execute(stmt).first() + + row.station_id = station_data[0].id + + return row + + +def fix_drifter_depths(ds): + if "DEPH" in ds.variables and len(ds.DEPH) == 2 and -12000 in ds.DEPH.data: + inv_depth_idx = np.argwhere(ds.DEPH.data != -12000) + ds.DEPH.values[~inv_depth_idx] = ds.DEPH.values[inv_depth_idx] + return ds + + +def format_attrs(attrs): + attrs = { + "Platform Code": attrs.get("platform_code"), + "WMO Platform Code": attrs.get("wmo_platform_code"), + "Platform Name": attrs.get("platform_name"), + "Institution": attrs.get("institution"), + "PI Name": attrs.get("creator_name"), + "Data Center": attrs.get("DC_REFERENCE"), + "Platform Type": attrs.get("wmo_instrument_type"), + "WMO Instrument Type": attrs.get("wmo_instrument_type"), + } + formatted_attrs = {} + for key, value in attrs.items(): + new_value = value + if isinstance(new_value, bytes): + new_value = new_value.decode("utf-8").strip() + if isinstance(new_value, str): + uft8 = new_value.encode('utf-8', errors='ignore') + if b"\xef\xbf\xbd" in uft8: + new_value = formatted_attrs["Platform Code"] + else: + new_value = uft8.decode("utf-8").strip() + formatted_attrs[key] = new_value + formatted_attrs = {key: value for key, value in formatted_attrs.items() if value and value.strip()} + return formatted_attrs + + +def import_obs_data(file, variables, platform_type): + ds = xr.open_dataset(file).drop_duplicates("TIME").sel(TIME=slice("2020-01-01",None)) + + variables = [v for v in variables if v in ds.variables] + qc_variables = [var for var in ds.data_vars if 'QC' in var] + if len(variables) == 0 or len(ds.TIME) == 0: + return pd.DataFrame(), {}, [] + + if platform_type == "DB": + ds = fix_drifter_depths(ds) + + df = ( + ds[["TIME", "LATITUDE", "LONGITUDE", *variables, *qc_variables]] + .to_dataframe() + .reset_index() + ) + + # filter out bad data based on QC values + for qc_var in qc_variables: + if "TIME" in qc_var or "POSITION" in qc_var: + df = df[~df[qc_var].isin([0, 3, 4, 6, 9])] + else: + idx = df.loc[df[qc_var].isin([0, 3, 4, 6, 9])].index + var = qc_var.replace("_QC", "") + df.loc[idx, var] = np.nan + + df.drop(columns=qc_variables, inplace=True) + + # calculate depths if missing + if "DEPTH" not in ds.variables and "DEPH" not in ds.variables and "PRES" in ds.variables: + df = df.dropna(subset="PRES", ignore_index=True) + if "DEPTH" not in df.columns: + df["DEPTH"] = 0.0 + df = df.astype({"DEPTH": "float64"}) + df.loc[:, "DEPTH"] = gsw.conversions.z_from_p(-df["PRES"], df["LATITUDE"]) + elif "DEPH" in ds.variables: + df["DEPTH"] = df["DEPH"] + + new_cols = { + "TIME": "time", + "DEPTH": "depth", + "LATITUDE": "latitude", + "LONGITUDE": "longitude", + } + for col in df.columns: + if col in variables: + new_cols[col] = ds[col].standard_name + + df.rename(columns=new_cols, inplace=True) + + if "depth" not in df.columns and platform_type == "DB": + df["depth"] = 0.0 + + if "TRAJECTORY" in df.columns: + df.drop(columns=["TRAJECTORY"], inplace=True) + + if "DEPH" in df.columns: + df.drop(columns=["DEPH"], inplace=True) + + df = pd.melt( + df, + id_vars=["time", "depth", "latitude", "longitude"], + var_name="datatype_key", + ).sort_values(by=["time", "depth"], ignore_index=True) + df.time = df.time.dt.floor("s") + df.dropna(inplace=True) + + var_metadata = [{ + "key": ds[v].standard_name, + "name": ds[v].long_name, + "unit": ds[v].units + } for v in variables] + + attrs = format_attrs(ds.attrs) + + return df, attrs, var_metadata + + +def add_datatypes(engine, variables, var_metadata): + with Session(engine) as session: + db_datatypes = session.execute(select(DataType)).all() + db_datatype_keys = [dt[0].key for dt in db_datatypes] + + for variable in var_metadata: + if variable["key"] not in db_datatype_keys: + try: + dt = DataType( + key=variable["key"], + name=variable["name"], + unit=variable["unit"], + ) + session.add(dt) + except IntegrityError as e: + print("Error committing datatype.") + print(e) + session.rollback() + + session.commit() + + +def add_platform(engine, attrs, platform_type): + + with Session(engine) as session: + stmt = select(Platform.id).where( + Platform.unique_id == attrs["Platform Code"] + ) + platform_id = session.execute(stmt).first() + + if platform_id is not None: + platform_id = platform_id[0] + else: + platform_schema = None + + match platform_type: + case "GL": + platform_schema = Platform.Type.glider + case "DB": + platform_schema = Platform.Type.drifter + case "PF": + platform_schema = Platform.Type.argo + case "CT": + platform_schema = Platform.Type.mission + + platform = Platform(type=platform_schema, unique_id=f"{attrs["Platform Code"]}") + platform.attrs = attrs + + with Session(engine) as session: + try: + session.add(platform) + session.commit() + platform_id = platform.id + except IntegrityError as e: + print("Error committing platform.") + print(e) + session.rollback() + + return platform_id + + +def add_stations(df, engine): + station_df = df[["platform_id", "time", "latitude", "longitude"]].copy() + station_df.drop_duplicates(inplace=True) + + with Session(engine) as session: + station_df.to_sql( + name="stations", + con=session.connection(), + if_exists="append", + index=False, + method=insert_stations, + ) + + station_df = station_df.assign(station_id=None) + station_df = station_df.apply(get_station_ids, axis=1, session=session) + + df = pd.merge(df, station_df) + + return df + + +def add_samples(df, engine): + samples_df = df[["station_id", "depth", "value", "datatype_key"]] + + with Session(engine) as session: + samples_df.to_sql( + name="samples", + con=session.connection(), + if_exists="append", + index=False, + method=insert_samples, + chunksize=10000 + ) + + +def correct_platform_tracks(platform_id : str, engine): + + def mysql_upsert(table, conn, keys, data_iter): + data = [dict(zip(keys, row)) for row in data_iter] + stmt = insert(table.table).values(data) + update_stmt = stmt.on_duplicate_key_update(**{c.name: c for c in stmt.inserted}) + conn.execute(update_stmt) + + stations = pd.read_sql( + f"SELECT * FROM stations where platform_id={platform_id};", + con=engine, + index_col="id", + ).sort_values(by="time") + + corrected_lon = np.copy(stations.longitude.values) + diffs = np.diff(corrected_lon) + crossings = np.where(np.abs(diffs) > 180)[0] + + if len(crossings) > 0: + print(f"Updating track of platform {platform_id}.") + + for crossing in crossings: + if diffs[crossing] > 0: + corrected_lon[crossing + 1 :] -= 360 + else: + corrected_lon[crossing + 1 :] += 360 + stations.longitude = corrected_lon + stations.to_sql( + "stations", engine, if_exists="append", method=mysql_upsert, index=False + ) + + +def import_cmems_obs(uri : str, file_list : list | PosixPath | str, platform_type : str): + + engine = create_engine( + uri, + connect_args={"connect_timeout": 10}, + pool_recycle=3600, + ) + + if isinstance(file_list, str) | isinstance(file_list, PosixPath): + file_list = [file_list] + + variables = get_platform_variables(platform_type) + + for file in file_list: + try: + print(f"Importing file: {file}") + df, attrs, var_metadata = import_obs_data(file, variables, platform_type) + if len(df) == 0: + print(f"No valid variables found in file: {file}. Skipping.") + continue + + add_datatypes(engine, variables, var_metadata) + + platform_id = add_platform(engine, attrs, platform_type) + df["platform_id"] = platform_id + + df = add_stations(df, engine) + add_samples(df, engine) + + correct_platform_tracks(platform_id, engine) + except Exception as e: + print(f"Error importing file: {file}") + print(e) + raise + engine.dispose() + + +def import_cmems_obs_mp(uri, file_list, platform_type): + inputs = list(itertools.product([uri], file_list, [platform_type])) + with multiprocessing.Pool(processes=4) as pool: + pool.starmap(import_cmems_obs, inputs) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + prog="CMEMS Obs script", + description="Add CMEMS Observation data to ONAV Obs database.", + ) + parser.add_argument("uri", type=str, help="URI of the observation database.") + parser.add_argument("filename", type=str, help="Path to the observation NetCDF file.") + parser.add_argument("platform_type", type=str, help="Type of platform: GL, DB, PF, CT.") + args = parser.parse_args() + + import_cmems_obs([args.uri, args.filename, args.platform_type]) diff --git a/scripts/data_importers/cmems_argo.py b/scripts/data_importers/cmems_argo.py deleted file mode 100644 index 7b37f0303..000000000 --- a/scripts/data_importers/cmems_argo.py +++ /dev/null @@ -1,218 +0,0 @@ -#!/usr/bin/env python -import glob -import os -import sys - -import defopt -import pandas as pd -import gsw -import numpy as np -import xarray as xr -from sqlalchemy import create_engine, select -from sqlalchemy.exc import IntegrityError -from sqlalchemy.orm import Session - -current = os.path.dirname(os.path.realpath(__file__)) -parent = os.path.dirname(os.path.dirname(current)) -sys.path.append(parent) - -from data.observational import DataType, Platform, Sample, Station - -VARIABLES = ["TEMP", "PSAL"] - -META_FIELDS = { - "PLATFORM_NUMBER": "wmo_platform_code", - "PROJECT_NAME": "platform_name", - "PI_NAME": "creator_name", - "DATA_CENTRE": "DC_REFERENCE", - "PLATFORM_TYPE": "wmo_instrument_type", - "WMO_INST_TYPE": "wmo_instrument_type", -} - -datatype_map = {} - - -def extract_metadata(ds, metadata): - platform_number = ds.attrs.get(metadata["PLATFORM_NUMBER"]) - project_name = ds.attrs.get(metadata["PROJECT_NAME"]) - pi_name = ds.attrs.get(metadata["PI_NAME"]) - data_centre = ds[metadata["DATA_CENTRE"]].values[0] - platform_type = ds.attrs.get(metadata["PLATFORM_TYPE"]) - wmo_inst_type = ds.attrs.get(metadata["WMO_INST_TYPE"]) - - return { - "PLATFORM_NUMBER": platform_number, - "PROJECT_NAME": project_name, - "PI_NAME": pi_name, - "DATA_CENTRE": data_centre, - "PLATFORM_TYPE": platform_type, - "WMO_INST_TYPE": wmo_inst_type, - } - - -def reformat_coordinates(ds: xr.Dataset) -> xr.Dataset: - """ - Shifts coordinates so that tracks are continuous on each side of map limits - (-180,180 degrees longitude). i.e if a track crosses -180 deg such that the - first point is -178 and the next is 178 then the second coordinate will be - replaced with -182. This allows the navigator to draw the track continusouly - without bounching between points on the far sides of the map. - """ - - lons = ds.LONGITUDE.data.copy() - - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - while len(crossings) > 0: - if lons[crossings[0]] > lons[crossings[0] + 1]: - lons[crossings[0] + 1 :] = 360 + lons[crossings[0] + 1 :] - else: - lons[crossings[0] + 1 :] = -360 + lons[crossings[0] + 1 :] - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - ds.LONGITUDE.data = lons - - return ds - - -def main(uri: str, filename: str): - """Import Argo Profiles - - :param str uri: Database URI - :param str filename: Argo NetCDF Filename, or directory of files - """ - - engine = create_engine( - uri, - connect_args={"connect_timeout": 10}, - pool_recycle=3600, - ) - - with Session(engine) as session: - - if not isinstance(filename, list): - if os.path.isdir(filename): - filenames = sorted(glob.glob(os.path.join(filename, "*.nc"))) - else: - filenames = [filename] - else: - filenames = filename - - for fname in filenames: - with xr.open_dataset(fname) as ds: - print(fname) - if ds.LATITUDE.size == 1: - print("Moored instrument: skipping file.") - continue - - ds = reformat_coordinates(ds) - - times = pd.to_datetime(ds.TIME.values) - - meta_data = extract_metadata(ds, META_FIELDS) - - variables = [v for v in VARIABLES if v in ds.variables] - - platform_number = meta_data["PLATFORM_NUMBER"] - unique_id = f"argo_{platform_number}" - - platform = ( - session.query(Platform) - .filter( - Platform.unique_id == unique_id, - Platform.type == Platform.Type.argo, - ) - .first() - ) - if platform is None: - platform = Platform(type=Platform.Type.argo, unique_id=unique_id) - attrs = {} - for f in META_FIELDS: - attrs[META_FIELDS[f]] = (meta_data[f] or "").strip() - - platform.attrs = attrs - session.add(platform) - session.commit() - - for idx, time in enumerate(times): - station = Station( - time=time, - latitude=ds["LATITUDE"].values[idx], - longitude=ds["LONGITUDE"].values[idx], - platform_id=platform.id, - ) - - try: - session.add(station) - session.commit() - except IntegrityError: - print("Error committing station.") - session.rollback() - - samples = [] - for variable in variables: - if variable not in datatype_map: - statement = select(DataType).where( - DataType.key == ds[variable].standard_name - ) - dt = session.execute(statement).all() - if not dt: - dt = DataType( - key=ds[variable].standard_name, - name=ds[variable].long_name, - unit=ds[variable].units, - ) - - session.add(dt) - session.commit() - datatype_map[variable] = dt - else: - dt = dt[0][0] - else: - dt = datatype_map[variable] - - values = ds[variable].isel(TIME=idx).values - - if "DEPH" in ds.variables: - depth = ds.DEPH.isel(TIME=idx).values - elif "PRES" in ds.variables: - pres = ds["PRES"].isel(TIME=idx).values - lat = ds["LATITUDE"][idx].values - - depth = gsw.conversions.z_from_p( - -pres, - lat, - ) - - data = np.stack( - [ - depth.flatten(), - values.flatten(), - ], - axis=1, - ) - data = data[~np.isnan(data).any(axis=1)] - - samples = [ - Sample( - depth=pair[0], - datatype_key=dt.key, - value=pair[1], - station_id=station.id, - ) - for pair in data - ] - - try: - session.bulk_save_objects(samples) - except IntegrityError: - print("Error committing samples.") - session.rollback() - - session.commit() - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/cmems_ctd.py b/scripts/data_importers/cmems_ctd.py deleted file mode 100644 index 9a99a9f6f..000000000 --- a/scripts/data_importers/cmems_ctd.py +++ /dev/null @@ -1,188 +0,0 @@ -#!/usr/bin/env python - -import glob -import os -import sys - -import defopt -import gsw -import numpy as np -import pandas as pd -import xarray as xr - -from sqlalchemy import create_engine, select -from sqlalchemy.exc import IntegrityError -from sqlalchemy.orm import Session - -current = os.path.dirname(os.path.realpath(__file__)) -parent = os.path.dirname(os.path.dirname(current)) -sys.path.append(parent) - -from data.observational import DataType, Platform, Sample, Station - -def reformat_coordinates(ds: xr.Dataset) -> xr.Dataset: - """ - Shifts coordinates so that tracks are continuous on each side of map limits - (-180,180 degrees longitude). i.e if a track crosses -180 deg such that the - first point is -178 and the next is 178 then the second coordinate will be - replaced with -182. This allows the navigator to draw the track continusouly - without bounching between points on the far sides of the map. - """ - - lons = ds.LONGITUDE.data.copy() - - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - while len(crossings) > 0: - if lons[crossings[0]] > lons[crossings[0] + 1]: - lons[crossings[0] + 1 :] = 360 + lons[crossings[0] + 1 :] - else: - lons[crossings[0] + 1 :] = -360 + lons[crossings[0] + 1 :] - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - ds.LONGITUDE.data = lons - - return ds - -def main(uri: str, filename: str): - - """Import NAFC CTD - - :param str uri: Database URI - :param str filename: NetCDF file, or directory of files - """ - engine = create_engine( - uri, - connect_args={"connect_timeout": 10}, - pool_recycle=3600, - ) - - with Session(engine) as session: - datatype_map = {} - - if not isinstance(filename, list): - if os.path.isdir(filename): - filenames = sorted(glob.glob(os.path.join(filename, "*.nc"))) - else: - filenames = [filename] - else: - filenames = filename - - for fname in filenames: - print(fname) - with xr.open_dataset(fname) as ds: - - if ds.LATITUDE.size == 1: - print("Moored instrument: skipping file.") - continue - - ds = reformat_coordinates(ds) - - times = pd.to_datetime(ds.TIME.values) - - if len(datatype_map) == 0: - # Generate the DataTypes; only consider variables that have depth - for var in filter( - lambda x, dataset=ds: "DEPTH" in dataset[x].dims, - [d for d in ds.variables], - ): - variable = ds[var] - standard_name = variable.attrs.get('standard_name') - if standard_name is not None: - statement = select(DataType).where( - DataType.key == ds[var].standard_name - ) - dt = session.execute(statement).all() - if not dt : - dt = DataType( - key=ds[var].standard_name, - name=ds[var].long_name, - unit=ds[var].units, - ) - session.add(dt) - else: - dt = dt[0][0] - - datatype_map[var] = dt - session.commit() - - p = Platform( - type=Platform.Type.mission, unique_id=f"{ds.attrs["platform_code"]}" - ) - attrs = { - "Institution": ds.attrs["institution"], - } - p.attrs = attrs - - try: - session.add(p) - session.commit() - except IntegrityError: - print("Error committing platform.") - session.rollback() - stmt = select(Platform.id).where(Platform.unique_id == ds.attrs["platform_code"]) - p.id = session.execute(stmt).first()[0] - - for idx, time in enumerate(times): - # Generate the station - s = Station( - latitude=ds.LATITUDE.values[idx], - longitude=ds.LONGITUDE.values[idx], - time=time, - platform_id=p.id - ) - try: - session.add(s) - session.commit() - except IntegrityError: - print("Error committing station.") - session.rollback() - - # Generate the samples - for var, dt in datatype_map.items(): - if "DEPH" in ds.variables: - depth = ds.DEPH.isel(TIME=idx).values - elif "PRES" in ds.variables: - pres = ds["PRES"].isel(TIME=idx).values - lat = ds["LATITUDE"][idx].values - - depth = gsw.conversions.z_from_p( - -pres, - lat, - ) - - if var in ds.variables: - values = ds[var].isel(TIME=idx).values - - data = np.stack( - [ - depth.flatten(), - values.flatten(), - ], - axis=1, - ) - data = data[~np.isnan(data).any(axis=1)] - - samples = [ - Sample( - depth=pair[0], - datatype_key=dt.key, - value=pair[1], - station_id=s.id, - ) - for pair in data - ] - - try: - session.bulk_save_objects(samples) - except IntegrityError: - print("Error committing samples.") - session.rollback() - - session.commit() - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/cmems_drifter.py b/scripts/data_importers/cmems_drifter.py deleted file mode 100644 index d6f5f2d55..000000000 --- a/scripts/data_importers/cmems_drifter.py +++ /dev/null @@ -1,157 +0,0 @@ -#!/usr/bin/env python - -import glob -import os -import sys -import defopt - -import numpy as np -import xarray as xr -from sqlalchemy import create_engine, select -from sqlalchemy.exc import IntegrityError -from sqlalchemy.orm import Session - -current = os.path.dirname(os.path.realpath(__file__)) -parent = os.path.dirname(os.path.dirname(current)) -sys.path.append(parent) - -from data.observational import DataType, Platform, Sample, Station - -# Mapping of variable names to their attributes -VARIABLE_MAPPING = { - "TEMP": ("sea_water_temperature", "Sea Temperature", "degrees_C"), - "ATMS": ("air_pressure_at_sea_level", "Atmospheric Pressure", "hPa"), - # Add more variables here if necessary -} - -def reformat_coordinates(ds: xr.Dataset) -> xr.Dataset: - """ - Shifts coordinates so that tracks are continuous on each side of map limits - (-180,180 degrees longitude). i.e if a track crosses -180 deg such that the - first point is -178 and the next is 178 then the second coordinate will be - replaced with -182. This allows the navigator to draw the track continusouly - without bounching between points on the far sides of the map. - """ - - lons = ds.LONGITUDE.data.copy() - - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - while len(crossings) > 0: - if lons[crossings[0]] > lons[crossings[0] + 1]: - lons[crossings[0] + 1 :] = 360 + lons[crossings[0] + 1 :] - else: - lons[crossings[0] + 1 :] = -360 + lons[crossings[0] + 1 :] - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - ds.LONGITUDE.data = lons - - return ds - -def main(uri: str, filename: str): - """Import data from NetCDF file(s) into the database. - - :param str uri: Database URI - :param str filename: NetCDF file or directory of NetCDF files - """ - - engine = create_engine( - uri, - connect_args={"connect_timeout": 10}, - pool_recycle=3600, - ) - with Session(engine) as session: - - if not isinstance(filename, list): - if os.path.isdir(filename): - filenames = sorted(glob.glob(os.path.join(filename, "*.nc"))) - else: - filenames = [filename] - else: - filenames = filename - - for fname in filenames: - print(fname) - with xr.open_dataset(fname) as ds: - - if ds.LATITUDE.size > 1: - ds = reformat_coordinates(ds) - - df = ds.to_dataframe().reset_index().dropna(axis=1, how="all").dropna() - - # Iterate over variables defined in the mapping - for var_name, (key, name, unit) in VARIABLE_MAPPING.items(): - if var_name not in df.columns: - continue # Skip if the variable is not present in the dataset - - # Create or retrieve DataType object - statement = select(DataType).where(DataType.key == key) - data_type = session.execute(statement).first() - if data_type is None: - data_type = DataType( - key=key, - name=name, - unit=unit, - ) - session.add(data_type) - session.commit() - else: - data_type = data_type[0] - - # Create Platform object - - platform = Platform(type=Platform.Type.drifter, unique_id=f"{ds.attrs["platform_code"]}") - try: - session.add(platform) - session.commit() - except IntegrityError: - print("Error committing platform.") - session.rollback() - stmt = select(Platform.id).where(Platform.unique_id == ds.attrs["platform_code"]) - platform.id = session.execute(stmt).first()[0] - - # Iterate over rows in the DataFrame - for index, row in df.iterrows(): - time = row["TIME"] - latitude = row["LATITUDE"] - longitude = row["LONGITUDE"] - depth = row["DEPH"] - value = row[var_name] - - # Create Station object - station = Station( - time=time, - latitude=latitude, - longitude=longitude, - platform_id=platform.id, - ) - - try: - session.add(station) - session.commit() - except IntegrityError: - print("Error committing station.") - session.rollback() - - # Create Sample object - sample = Sample( - depth=depth, - value=value, - datatype_key=data_type.key, - station_id=station.id, - ) - - try: - session.add(sample) - session.commit() - except IntegrityError: - print("Error committing sample.") - session.rollback() - - session.commit() - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/cmems_glider.py b/scripts/data_importers/cmems_glider.py deleted file mode 100644 index e6fdbfe22..000000000 --- a/scripts/data_importers/cmems_glider.py +++ /dev/null @@ -1,187 +0,0 @@ -#!/usr/bin/env python - -import glob -import os -import sys - -import defopt -import numpy as np -import xarray as xr -from sqlalchemy import create_engine, select -from sqlalchemy.exc import IntegrityError -from sqlalchemy.orm import Session - -current = os.path.dirname(os.path.realpath(__file__)) -parent = os.path.dirname(os.path.dirname(current)) -sys.path.append(parent) - -from data.observational import DataType, Platform, Sample, Station - -VARIABLES = ["PRES", "PSAL", "TEMP", "CNDC"] - - -def reformat_coordinates(ds: xr.Dataset) -> xr.Dataset: - """ - Shifts coordinates so that tracks are continuous on each side of map limits - (-180,180 degrees longitude). i.e if a track crosses -180 deg such that the - first point is -178 and the next is 178 then the second coordinate will be - replaced with -182. This allows the navigator to draw the track continusouly - without bounching between points on the far sides of the map. - """ - - lons = ds.LONGITUDE.data.copy() - - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - while len(crossings) > 0: - if lons[crossings[0]] > lons[crossings[0] + 1]: - lons[crossings[0] + 1 :] = 360 + lons[crossings[0] + 1 :] - else: - lons[crossings[0] + 1 :] = -360 + lons[crossings[0] + 1 :] - lon_diff = np.diff(lons) - crossings = np.where(np.abs(lon_diff) > 180)[0] - - ds.LONGITUDE.data = lons - - return ds - -def main(uri: str, filename: str): - - """Import Glider NetCDF - - :param str uri: Database URI - :param str filename: Glider Filename, or directory of NetCDF files - """ - engine = create_engine( - uri, - connect_args={"connect_timeout": 10}, - pool_recycle=3600, - ) - - with Session(engine) as session: - - if not isinstance(filename, list): - if os.path.isdir(filename): - filenames = sorted(glob.glob(os.path.join(filename, "*.nc"))) - else: - filenames = [filename] - else: - filenames = filename - - datatype_map = {} - for fname in filenames: - print(fname) - with xr.open_dataset(fname).drop_duplicates("TIME") as ds: - ds = reformat_coordinates(ds) - variables = [v for v in VARIABLES if v in ds.variables] - df = ( - ds[["TIME", "LATITUDE", "LONGITUDE", *variables]] - .to_dataframe() - .reset_index() - .dropna(axis=1, how='all') - .dropna() - ) - - # remove missing variables from variables list - variables = [v for v in VARIABLES if v in df.columns] - - for variable in variables: - if variable not in datatype_map: - statement = select(DataType).where( - DataType.key == ds[variable].standard_name - ) - dt = session.execute(statement).all() - if not dt: - dt = DataType( - key=ds[variable].standard_name, - name=ds[variable].long_name, - unit=ds[variable].units, - ) - session.add(dt) - else: - dt = dt[0][0] - - datatype_map[variable] = dt - - session.commit() - - p = Platform( - type=Platform.Type.glider, unique_id=f"{ds.attrs["platform_code"]}" - ) - attrs = { - "Glider Platform": ds.attrs["platform_code"], - "WMO": ds.attrs["wmo_platform_code"], - "Institution": ds.attrs["institution"], - } - p.attrs = attrs - - try: - session.add(p) - session.commit() - except IntegrityError: - print("Error committing platform.") - session.rollback() - stmt = select(Platform.id).where(Platform.unique_id == ds.attrs["platform_code"]) - p.id = session.execute(stmt).first()[0] - - n_chunks = np.ceil(len(df)/1e4) - - if n_chunks < 1: - continue - - for chunk in np.array_split(df, n_chunks): - chunk["STATION_ID"] = 0 - - stations = [ - Station( - platform_id=p.id, - time=row.TIME, - latitude=row.LATITUDE, - longitude=row.LONGITUDE, - ) - for idx, row in chunk[["TIME","LATITUDE", "LONGITUDE"]].drop_duplicates().iterrows() - ] - - # Using return_defaults=True here so that the stations will get - # updated with id's. It's slower, but it means that we can just - # put all the station ids into a pandas series to use when - # constructing the samples. - try: - session.bulk_save_objects(stations, return_defaults=True) - session.commit() - except IntegrityError: - print("Error committing station.") - session.rollback() - - stmt = select(Station).where(Station.platform_id==p.id) - station_data = session.execute(stmt).all() - - for station in station_data: - chunk.loc[chunk["TIME"]==station[0].time,"STATION_ID"] = station[0].id - - samples = [ - [ - Sample( - station_id=row.STATION_ID, - depth=row.DEPTH, - value=row[variable], - datatype_key=datatype_map[variable].key, - ) - for variable in variables - ] - for idx, row in chunk.iterrows() - ] - try: - session.bulk_save_objects( - [item for sublist in samples for item in sublist] - ) - except IntegrityError: - print("Error committing samples.") - session.rollback() - - session.commit() - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/concepts_drifter.py b/scripts/data_importers/concepts_drifter.py deleted file mode 100755 index f5c3f5690..000000000 --- a/scripts/data_importers/concepts_drifter.py +++ /dev/null @@ -1,124 +0,0 @@ -#!/usr/bin/env python - -import glob -import os - -import defopt -import numpy as np -import pandas as pd -import xarray as xr - -import data.observational -from data.observational import DataType, Platform, Sample, Station - -# These files were created by our scripts so they don't have nice attributes -# like standard_name, etc. -# Ideally this script won't stick around as anything other than an example for -# bulk loading drifter data. Going forward a new importer should be written to -# read from the incoming csv files and just add the new samples. -DATATYPE_MAPPING = { - "sst": ("sea_water_temperature", "Water Temperature", "degree_Celsius"), - "vbat": ("battery_voltage", "Battery Voltage", "Volts"), - "bp": ("air_pressure", "Barometric Pressure", "mbar"), - "bpt": ("tendency_of_air_pressure", "Barometric Pressure Tendency", "mbar"), -} - - -def main(uri: str, filename: str): - """Import CONCEPTS drifter NetCDF - - :param str uri: Database URI - :param str filename: Drifter Filename, or directory of NetCDF files - """ - data.observational.init_db(uri, echo=False) - data.observational.create_tables() - - if os.path.isdir(filename): - filenames = sorted(glob.glob(os.path.join(filename, "*.nc"))) - else: - filenames = [filename] - - for fname in filenames: - print(fname) - with xr.open_dataset(fname) as ds: - df = ds.to_dataframe().drop(["wmo", "deployment", "imei"], axis=1) - columns = list(filter(lambda c: c in DATATYPE_MAPPING, df.columns)) - - dt_map = {} - for c in columns: - # First check our local cache for the DataType object, if - # that comes up empty, check the db, and failing that, - # create a new one. - if c not in dt_map: - dt = DataType.query.get(DATATYPE_MAPPING[c][0]) - if dt is None: - dt = DataType( - key=DATATYPE_MAPPING[c][0], - name=DATATYPE_MAPPING[c][1], - unit=DATATYPE_MAPPING[c][2], - ) - - data.observational.db.session.add(dt) - - dt_map[c] = dt - - # Commit to make sure all the variables are in the db so we don't - # get any foreign key errors - data.observational.db.session.commit() - - p = Platform(type=Platform.Type.drifter) - attrs = dict(ds.attrs) - attrs["wmo"] = ds.wmo.values[0] - attrs["deployment"] = ds.deployment.values[0] - attrs["imei"] = ds.imei.values[0] - p.attrs = attrs - data.observational.db.session.add(p) - data.observational.db.session.commit() - - samples = [] - for index, row in df.iterrows(): - time = index[0] - lat = row["latitude"] - lon = row["longitude"] - - station = Station( - time=time, latitude=lat, longitude=lon, platform_id=p.id - ) - data.observational.db.session.bulk_save_objects( - [station], return_defaults=True - ) - - for c in columns: - value = row[c] - if isinstance(value, pd.Timestamp): - value = value.value / 10**9 - - if np.isfinite(value): - samples.append( - Sample( - depth=0, - datatype_key=DATATYPE_MAPPING[c][0], - value=value, - station_id=station.id, - ) - ) - - # Commit every 1000 samples, that's a decent balance between - # locking the db for too long and performance - if len(samples) > 1000: - data.observational.db.session.bulk_save_objects(samples) - data.observational.db.session.commit() - samples = [] - - # If there are any samples that haven't been committed yet, do so - # now. - if samples: - data.observational.db.session.bulk_save_objects(samples) - data.observational.db.session.commit() - samples = [] - - data.observational.db.session.commit() - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/init_cmems_obs.py b/scripts/data_importers/init_cmems_obs.py deleted file mode 100644 index 1cadfb13a..000000000 --- a/scripts/data_importers/init_cmems_obs.py +++ /dev/null @@ -1,52 +0,0 @@ -import defopt -import os - -import copernicusmarine - -import cmems_argo -import cmems_drifter -import cmems_glider -import cmems_ctd - - -def main(uri: str, output_dir: str): - """ - Downloads historical CMEMS observation data for the observation types - given in obs_types to initialize the observation database. Deletes - dowloaded files after data is added. - Arguments should be passed on commandline. e.g. - python scripts/data_importers/init_cmems_obs.py "mysql://usr:pwd@address/database" "/data/nrt.cmems-du.eu/" - - :param uri: The URI string of the MariaDB Observation database - :param output_dir: The directory that you want to save the data to. - """ - - obs_types = ["GL", "DB", "PF", "CT"] - - for obs in obs_types: - obs_filter = f"{obs}/*.nc" - file_list = copernicusmarine.get( - dataset_id="cmems_obs-ins_glo_phybgcwav_mynrt_na_irr", - output_directory=output_dir, - filter=obs_filter, - force_download=True, - overwrite_output_data=True, - dataset_part="history", - ) - - # Call the appropriate function for the observation type - if "GL" in obs_filter: - cmems_glider.main(uri, file_list) - elif "DB" in obs_filter: - cmems_drifter.main(uri, file_list) - elif "PF" in obs_filter: - cmems_argo.main(uri, file_list) - elif "CT" in obs_filter: - cmems_ctd.main(uri, file_list) - - for file in file_list: - os.remove(file) - - -if __name__ == "__main__": - defopt.run(main) diff --git a/scripts/data_importers/sync_cmems_obs.py b/scripts/data_importers/sync_cmems_obs.py deleted file mode 100644 index a424eae17..000000000 --- a/scripts/data_importers/sync_cmems_obs.py +++ /dev/null @@ -1,50 +0,0 @@ -import defopt -from datetime import datetime - -import copernicusmarine - -import cmems_argo -import cmems_drifter -import cmems_glider -import cmems_ctd - - -def main(uri: str, output_dir: str): - """ - Downloads today's data for the observation types given in obs_types then adds that - data to the observation database. - Arguments should be passed on commandline. e.g. - python scripts/data_importers/sync_cmems_obs.py "mysql://usr:pwd@address/database" "/data/nrt.cmems-du.eu/" - - :param uri: The URI string of the MariaDB Observation database - :param output_dir: The directory that you want to save the data to. - """ - - obs_types = ["GL", "DB", "PF", "CT"] - - timestamp = datetime.today().strftime("%Y%m%d") - - for obs in obs_types: - obs_filter = f"{timestamp}/*_*_{obs}*.nc" - file_list = copernicusmarine.get( - dataset_id="cmems_obs-ins_glo_phybgcwav_mynrt_na_irr", - output_directory=output_dir, - filter=obs_filter, - force_download=True, - overwrite_output_data=True, - dataset_part="latest", - ) - - # Call the appropriate function for the observation type - if "GL" in obs_filter: - cmems_glider.main(uri, file_list) - elif "DB" in obs_filter: - cmems_drifter.main(uri, file_list) - elif "PF" in obs_filter: - cmems_argo.main(uri, file_list) - elif "CT" in obs_filter: - cmems_ctd.main(uri, file_list) - - -if __name__ == "__main__": - defopt.run(main)