diff --git a/environment-dev-geod-lev.yml b/environment-dev-geod-lev.yml new file mode 100644 index 00000000..0395790e --- /dev/null +++ b/environment-dev-geod-lev.yml @@ -0,0 +1,36 @@ +name: fire-dev-geod-lev +channels: + - conda-forge +dependencies: + - black + - click + - fiona + - gama + - graphviz + - matplotlib + - numpy + - openpyxl + - oracledb + - pandas + - pyarrow + - pip + - pre-commit + - pylint + - pyproj + - pytest-cov + - pytest-freezegun + - pytest-mock + - pytest + - python= + - rich + - scipy + - shapely + - sphinx-click + - sphinx + - sphinx_rtd_theme + - sphinxcontrib-programoutput + - sqlalchemy + - xmltodict + - astropy + - jplephem + - cartopy \ No newline at end of file diff --git a/src/fire/api/geodetic_levelling/__init__.py b/src/fire/api/geodetic_levelling/__init__.py new file mode 100644 index 00000000..bd3b0960 --- /dev/null +++ b/src/fire/api/geodetic_levelling/__init__.py @@ -0,0 +1,13 @@ +""" +geodetic_levelling +""" + +import fire.api.geodetic_levelling.contour_plot +import fire.api.geodetic_levelling.geodetic_correction_levelling_obs +import fire.api.geodetic_levelling.histogram +import fire.api.geodetic_levelling.metric_to_gpu_transformation +import fire.api.geodetic_levelling.tidal_transformation +import fire.api.geodetic_levelling.time_propagation +import fire.api.geodetic_levelling.geophysical_parameters + +__version__ = "1.0.0" diff --git a/src/fire/api/geodetic_levelling/contour_plot.py b/src/fire/api/geodetic_levelling/contour_plot.py new file mode 100644 index 00000000..5027d8f9 --- /dev/null +++ b/src/fire/api/geodetic_levelling/contour_plot.py @@ -0,0 +1,171 @@ +"""This module contains functions for the generation of contour plots showing +3rd precision levelling recalculated vs. DVR90 +""" + +from pathlib import Path + +from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER +import cartopy +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker +import matplotlib.tri as tri +import numpy as np +import pandas as pd + + +def extract_heights_from_fire_project( + fire_project: str, + excel_inputfolder: Path, +) -> pd.DataFrame: + """Extract heights from a FIRE project. + + Extracts adjusted heights and database heights from a FIRE project and returns the + extracted heights and corresponding height differences in a DataFrame. + + Note that only the heights of DVR90 defining points are extracted from the FIRE project. + + Args: + fire_project: str, name of FIRE project with adjusted heights and database heights, + e.g. "3prs" + excel_inputfolder: Path, folder containing FIRE project/excel file + + Returns: + pd.DataFrame, a DataFrame with extracted heights and height differences + + Raises: + ? Hvis input mapper eller filer ikke findes + + Input file: + FIRE project/excel file with adjusted heights and database heights, e.g. "3prs.xlsx" + """ + excel_inputfile = excel_inputfolder / f"{fire_project}.xlsx" + + punktoversigt_df = pd.read_excel(excel_inputfile, sheet_name="Punktoversigt") + kontrolberegning_df = pd.read_excel(excel_inputfile, sheet_name="Kontrolberegning") + + # DataFrame with DVR90 defining points in sheet "Punktoversigt" + # (i.e. all G.I. or G.M. points with a database height calculated 2000-02-11 13:30:00) + punktoversigt_df = punktoversigt_df[ + (punktoversigt_df["Hvornår"] == "2000-02-11 13:30:00") + & ( + punktoversigt_df["Punkt"].str.startswith("G.I.") + | punktoversigt_df["Punkt"].str.startswith("G.M.") + ) + ] + + # DataFrame with point id, database height and latitude/longitude of DVR90 defining points + punktoversigt_df = punktoversigt_df[["Punkt", "Kote", "Nord", "Øst"]] + + # DataFrame with point id and adjusted heights of all points in sheet "Kontrolberegning" + kontrolberegning_df = kontrolberegning_df[["Punkt", "Ny kote"]] + + kontrolberegning_df = kontrolberegning_df.set_index("Punkt") + + # DataFrame with point id, database height, adjusted height and latitude/longitude of + # DVR90 defining points + points_df = punktoversigt_df.join(kontrolberegning_df, on="Punkt", how="inner") + + # Add column with height differences (adjusted height - database height) to DataFrame + points_df["Diff"] = points_df["Ny kote"] - points_df["Kote"] + + return points_df + + +def contour_plot_recalc_vs_dvr90( + fire_project: str, + excel_inputfolder: Path, + outputfolder: Path, +) -> None: + """Generate contour plot showing 3rd precision levelling recalculated vs. DVR90. + + Generates contour plot showing 3rd precision levelling recalculated vs. DVR90 + (recalculated height - DVR90 height). + + Args: + fire_project: str, name of fire project with recalculated/adjusted heights and + DVR90/database heights, e.g. "3prs" + excel_inputfolder: Path, folder containing FIRE project/excel file + outputfolder: Path, folder for output contour plot + + Returns: + None + + Raises: + ? + + Input file: + FIRE project/excel file with adjusted heights and database heights, e.g. "3prs.xlsx" + + Output file: + Png file with contour plot + """ + # Make sure that the output folder exists + outputfolder.mkdir(parents=True, exist_ok=True) + + # DataFrame with height differences of DVR90 defining points + # (adjusted heights - database heights) + points_df = extract_heights_from_fire_project(fire_project, excel_inputfolder) + + # Generation of contour plot + fig = plt.figure(figsize=(7, 4), dpi=500) + crs = ccrs.PlateCarree() + ax = plt.axes(projection=crs) + + ax.set_title("3rd precision levelling recalculated vs. DVR90") + ax.add_feature(cartopy.feature.OCEAN, zorder=0) + ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor="black") + ax.add_feature(cartopy.feature.COASTLINE) + + # Levels for contour plot + min_level = np.floor(points_df["Diff"].min() * 1000) + max_level = np.ceil(points_df["Diff"].max() * 1000) + levels = np.arange(min_level, max_level + 1) + + triangulation = tri.Triangulation(points_df["Øst"], points_df["Nord"]) + triangles = triangulation.triangles + + # Mask off unwanted triangles across the Great Belt + east_tri = points_df.Øst.values[triangles] - np.roll( + points_df.Øst.values[triangles], 1, axis=1 + ) + north_tri = points_df.Nord.values[triangles] - np.roll( + points_df.Nord.values[triangles], 1, axis=1 + ) + maxi = np.max(np.sqrt(east_tri**2 + north_tri**2), axis=1) + + triangulation.set_mask(maxi > 0.80) + + ax.tricontour( + triangulation, + points_df["Diff"] * 1000, + levels, + linewidths=0.5, + colors="k", + transform=crs, + ) + + contour_regions = ax.tricontourf( + triangulation, + points_df["Diff"] * 1000, + levels, + transform=crs, + ) + + ax.plot(points_df["Øst"], points_df["Nord"], "ko", ms=2, transform=crs) + + fig.colorbar(contour_regions, ax=ax, label=r"$\Delta$H [mm] (recalculated - DVR90)") + + gl = ax.gridlines(crs=crs, draw_labels=True) + gl.top_labels = False + gl.right_labels = False + gl.xlocator = mticker.FixedLocator([8, 9, 10, 11, 12, 13, 14]) + gl.ylocator = mticker.FixedLocator([55, 56, 57, 58, 59]) + gl.xformatter = LONGITUDE_FORMATTER + gl.yformatter = LATITUDE_FORMATTER + + ax.set_xlim([7.9, 13.1]) + ax.set_ylim([54.3, 58.1]) + + plt.savefig(outputfolder / f"{fire_project}.png") + plt.close() diff --git a/src/fire/api/geodetic_levelling/geodetic_correction_levelling_obs.py b/src/fire/api/geodetic_levelling/geodetic_correction_levelling_obs.py new file mode 100644 index 00000000..e28c2a5b --- /dev/null +++ b/src/fire/api/geodetic_levelling/geodetic_correction_levelling_obs.py @@ -0,0 +1,283 @@ +"""This module contains functions for geodetic correction of height differences/levelling observations.""" + +from pathlib import Path + +import pandas as pd + +from fire.api.geodetic_levelling.tidal_transformation import ( + apply_tidal_corrections_to_height_diff, +) + +from fire.api.geodetic_levelling.time_propagation import ( + propagate_height_diff_from_epoch_to_epoch, +) + +from fire.api.geodetic_levelling.metric_to_gpu_transformation import ( + convert_metric_height_diff_to_geopotential_height_diff, +) + + +def apply_geodetic_corrections_to_height_diffs( + observations_df: pd.DataFrame, + points_df: pd.DataFrame, + grid_inputfolder: Path = None, + height_diff_unit: str = "metric", + tidal_system: str = None, + epoch_target: pd.Timestamp = None, + deformationmodel: str = None, + gravitymodel: str = None, +) -> tuple[pd.DataFrame, pd.DataFrame]: + """Apply various geodetic corrections to the height differences of a FIRE project. + + Applies various geodetic corrections to the height differences of a FIRE project. + + The metric height differences of a FIRE project are tidally corrected if and only if + the function is called with an argument for parameter tidal_system. + + The metric height differences of a FIRE project are propagated to a target epoch if and only if + the function is called with arguments for all three parameters epoch_target, deformationmodel + and grid_inputfolder. + + The metric height differences of a FIRE project are converted to geopotential units if and only + if the function is called with argument "gpu" for parameter height_diff_unit and with arguments + for both parameter gravitymodel and grid_inputfolder. + + Args: + fire_project: str, name of fire project with metric height differences to be corrected/converted, + e.g. "asmei_temp" + excel_inputfolder: Path, folder containing fire project/excel file + outputfolder: Path, folder for output fire project/excel file with corrected/converted height differences + grid_inputfolder: Path = None, optional parameter, folder for input grid, i.e. deformation model + and/or gravity model + height_diff_unit: str = "metric", optional parameter, determines whether or not metric + input height differences are converted to geopotential units, "metric" for no conversion, + "gpu" for conversion to gpu, default value is "metric" + tidal_system: str = None, optional parameter, system for tidal corrections of metric height + differences, "non", "mean" or "zero" for non-tidal, mean tide or zero tide + epoch_target: pd.Timestamp = None, optional parameter, target epoch for the propagation + of metric height differences (format: yyyy-mm-dd hh:mm:ss) + deformationmodel: str = None, optional parameter, deformation model used for the propagation + of metric height differences, must be in GeoTIFF or GTX file format, e.g. "NKG2016_lev.tif" + gravitymodel: str = None, optional parameter, gravity model used for the conversion of metric + height differences to gpu, must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + None + + Raises: + ? Hvis input mapper eller filer ikke findes + + Input file: + FIRE project/excel-file with height differences to be corrected/converted, e.g. "asmei_temp.xlsx" + + Output file: + Fire project/excel file with corrected/converted height differences + + TO DO: Samle "underparametre" i en eller flere dicts, fx tidal_parameters: dict = {}, + """ + + # TO DO: Flyt if-sætningerne, der kontrollerer hvilke korrektioner der foretages, + # foran for-loopet over observations_df.index og loop i stedet op til 3 gange over + # observations_df.index? + + # print(observations_df) + # print(type(observations_df)) + # print(points_df) + # print(type(points_df)) + # print(grid_inputfolder) + # print(type(grid_inputfolder)) + # print(height_diff_unit) + # print(type(height_diff_unit)) + # print(tidal_system) + # print(type(tidal_system)) + # print(epoch_target) + # print(type(epoch_target)) + # print(deformationmodel) + # print(type(deformationmodel)) + # print(gravitymodel) + # print(type(gravitymodel)) + + # DataFrame for applied geodetic corrections + corrections_df = observations_df[["Journal", "Fra", "Til"]].copy() + + for index in observations_df.index: + height_diff = observations_df.at[index, "ΔH"] + point_from = observations_df.at[index, "Fra"] + point_to = observations_df.at[index, "Til"] + epoch_obs = observations_df.at[index, "Hvornår"] + + # print(height_diff) + + # Geographic coordinates of point_from and point_to + # If no information on geographic coordinates is available the observation is skipped + try: + # fmt: off + point_from_lat = points_df.loc[points_df["Punkt"] == point_from, "Nord"].values[0] + point_from_long = points_df.loc[points_df["Punkt"] == point_from, "Øst"].values[0] + point_to_lat = points_df.loc[points_df["Punkt"] == point_to, "Nord"].values[0] + point_to_long = points_df.loc[points_df["Punkt"] == point_to, "Øst"].values[0] + # fmt: on + + except IndexError: + observations_df.at[index, "Sluk"] = "x" + + print( + "index:", + index, + "\n", + "point_from:", + point_from, + "\n", + "point_to:", + point_to, + "\n", + "The geographic coordinates of point_from and/or point_to are missing.\n", + "The observation is skipped.", + ) + + continue + + # The metric height differences of a FIRE project are tidally corrected if the + # function apply_geodetic_corrections_to_height_diffs is called with an argument for + # parameter tidal_system + if tidal_system is not None: + (height_diff, tidal_corr) = apply_tidal_corrections_to_height_diff( + height_diff, + point_from_lat, + point_from_long, + point_to_lat, + point_to_long, + epoch_obs, + tidal_system, + grid_inputfolder=grid_inputfolder, + gravitymodel=gravitymodel, + ) + # print(height_diff) + # observations_df.at[ + # index, f"ΔH tidal correction (tidal system: {tidal_system}) [m]" + # ] = tidal_corr + + corrections_df.at[ + index, f"ΔH tidal correction (tidal system: {tidal_system}) [m]" + ] = tidal_corr + + # The metric height differences of a FIRE project are propagated to a target epoch if + # the function apply_geodetic_corrections_to_height_diffs is called with arguments for + # all three parameters epoch_target, deformationmodel and grid_inputfolder + if ( + (epoch_target is not None) + and (deformationmodel is not None) + and (grid_inputfolder is not None) + ): + (height_diff, epoch_corr) = propagate_height_diff_from_epoch_to_epoch( + height_diff, + point_from_lat, + point_from_long, + point_to_lat, + point_to_long, + epoch_obs, + epoch_target, + grid_inputfolder, + deformationmodel, + ) + # print(height_diff) + # observations_df.at[ + # index, f"ΔH epoch correction (target epoch: {epoch_target}) [m]" + # ] = epoch_corr + + corrections_df.at[ + index, f"ΔH epoch correction (target epoch: {epoch_target}) [m]" + ] = epoch_corr + + # The metric height differences of a FIRE project are converted to geopotential units if + # the function apply_geodetic_corrections_to_height_diffs is called with argument "gpu" + # for parameter height_diff_unit and with arguments for both parameter gravitymodel + # and grid_inputfolder + if ( + height_diff_unit == "gpu" + and (gravitymodel is not None) + and (grid_inputfolder is not None) + ): + (height_diff, m2gpu_factor) = ( + convert_metric_height_diff_to_geopotential_height_diff( + height_diff, + point_from_lat, + point_from_long, + point_to_lat, + point_to_long, + tidal_system, + grid_inputfolder, + gravitymodel, + ) + ) + + # observations_df.at[ + # index, + # f"ΔH m2gpu multiplication factor (tidal system: {tidal_system}) [m/s^2]", + # ] = m2gpu_factor + + corrections_df.at[ + index, + f"ΔH m2gpu multiplication factor (tidal system: {tidal_system}) [m/s^2]", + ] = m2gpu_factor + + elif height_diff_unit == "metric": + pass + + else: + exit( + "Function apply_geodetic_corrections_to_height_diffs: Wrong arguments for\n\ + parameter height_diff_unit and/or gravitymodel and/or grid_inputfolder." + ) + # print(height_diff) + # Update of observations_df with corrected height difference + observations_df.at[index, "ΔH"] = height_diff + + return (observations_df, corrections_df) + + # DataFrame with parameters of output fire project + parameters_df = pd.read_excel(excel_inputfile, sheet_name="Parametre") + + parameters_new_df = pd.DataFrame( + { + "Navn": [ + "Tidal system", + "Target epoch", + "Unit of height differences", + "Deformationmodel", + "Gravitymodel", + ], + "Værdi": [ + tidal_system, + epoch_target, + height_diff_unit, + deformationmodel, + gravitymodel, + ], + }, + ) + + parameters_df = pd.concat([parameters_df, parameters_new_df], ignore_index=True) + + # Generation of output fire project/excel file with corrected/converted height differences + with pd.ExcelWriter( + outputfolder / f"{fire_project}.xlsx" + ) as writer: # pylint: disable=abstract-class-instantiated + pd.read_excel(excel_inputfile, sheet_name="Projektforside").to_excel( + writer, "Projektforside", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Sagsgang").to_excel( + writer, "Sagsgang", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Nyetablerede punkter").to_excel( + writer, "Nyetablerede punkter", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Notater").to_excel( + writer, "Notater", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Filoversigt").to_excel( + writer, "Filoversigt", index=False + ) + observations_df.to_excel(writer, "Observationer", index=False) + points_df.to_excel(writer, "Punktoversigt", index=False) + parameters_df.to_excel(writer, "Parametre", index=False) diff --git a/src/fire/api/geodetic_levelling/geophysical_parameters.py b/src/fire/api/geodetic_levelling/geophysical_parameters.py new file mode 100644 index 00000000..3c1dcd52 --- /dev/null +++ b/src/fire/api/geodetic_levelling/geophysical_parameters.py @@ -0,0 +1,53 @@ +"""This module contains various geophysical parameters and constants required +for tidal correction and transformation of gravity and height etc. +""" + +from math import pi + +# Inclination of the ecliptic/the lunar orbit in radians +epsilon = (23 + (27 / 60)) * (1 / 360) * 2 * pi + +# Geocentric distance to the Moon in units of m +moon_dist = 3.84399 * 1e8 + +# Mass of the Moon in units of kg +moon_mass = 7.346 * 1e22 + +# Mean radius of the Earth in units of m +radius_earth = 6371000 + +# Love numbers +h = 0.62 +k = 0.30 + +# Tilt factor +gamma = 1 + k - h + +# Gravimetric factor +delta = 1 + h - (3 / 2) * k + +# Defining constants for Geodetic Reference System 1980 (GRS80) +# Major semi-axis of the reference ellipsoid in units of m +a_GRS80 = 6378137 + +# Gravitational mass constant of the Earth in units of m^3/s^2 +GM_GRS80 = 3986005 * 1e8 + +# Dynamic form factor +J2_GRS80 = 108263 * 1e-8 + +# Angular velocity of the Earth’s rotation in units of rad/s +omega_GRS80 = 7292115 * 1e-11 + +# Derived constants for Geodetic Reference System 1980 (GRS80) +# Minor semi-axis of the reference ellipsoid in units of m +b_GRS80 = 6356752.3141 + +# Flattening of the reference ellipsoid +f_GRS80 = 0.00335281068118 + +# m = (normal_gravity^2*a^2*b)/(G*M) +m_GRS80 = 0.00344978600308 + +# Normal gravity at equator in units of m/s^2 +normal_gravity_equator_GRS80 = 9.7803267715 diff --git a/src/fire/api/geodetic_levelling/histogram.py b/src/fire/api/geodetic_levelling/histogram.py new file mode 100644 index 00000000..c4e1809f --- /dev/null +++ b/src/fire/api/geodetic_levelling/histogram.py @@ -0,0 +1,102 @@ +"""This module contains functions for generation of histograms visualizing +the temporal distribution of levelling observations in a FIRE project. +""" + +from pathlib import Path +import datetime + +from astropy.time import Time +from matplotlib.ticker import MaxNLocator +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +def datetime_to_decimalyear(datetime: datetime.datetime | pd.Timestamp) -> np.float64: + """Convert a datetime to decimal year. + + Converts a datetime to decimal year. + + Args: + datetime: datetime.datetime | pd.Timestamp, input datetime or Timestamp to be converted + + Returns: + np.float64, input datetime converted to decimal year + + Raises: + ? + """ + return Time(datetime, format="datetime").decimalyear + + +def generate_histogram_temporal_distr_levelling_obs( + fire_project: str, + excel_inputfolder: Path, + outputfolder: Path, + bins: list[int | float] = None, +) -> None: + """Generate histogram visualizing the temporal distribution of levelling observations + in a FIRE project. + + Generates histogram visualizing the temporal distribution of levelling observations + in a FIRE project. + + Args: + fire_project: str, name of FIRE project with levelling observations, e.g. "3prs" + excel_inputfolder: Path, folder containing FIRE project/excel file + outputfolder: Path, folder for output histogram + bins: list[int | float] = None, optional parameter, list defining bins for histogram + + Returns: + None + + Raises: + ? + + Input file: + FIRE project/excel file with levelling observations, e.g. "3prs.xlsx" + + Output file: + Png file with histogram + """ + # Make sure that the output folder exists + outputfolder.mkdir(parents=True, exist_ok=True) + + excel_inputfile = excel_inputfolder / f"{fire_project}.xlsx" + + observationer_df = pd.read_excel(excel_inputfile, sheet_name="Observationer") + + # Observation epochs (date and time) in a Series + observation_epochs = observationer_df["Hvornår"].dropna() + + # Conversion of observation epochs to decimal year + observation_epochs = observation_epochs.apply(datetime_to_decimalyear) + + fig, ax = plt.subplots(1, 1, sharey=True, tight_layout=True) + str = f"mean = {observation_epochs.mean():4.2f}" + ax.hist(observation_epochs, bins=bins, color="red") + plt.axvline( + observation_epochs.mean(), color="k", linestyle="dashed", linewidth=1, label=str + ) + fig.suptitle( + f"Observation epochs {fire_project}", + fontsize=16, + ) + ax.set( + xlabel="Observation epochs [year]", + ylabel="#Observations", + title="Observation epochs", + ) + ax.yaxis.set_major_locator(MaxNLocator(integer=True)) + ax.legend( + loc="upper right", + fontsize=8, + fancybox=False, + framealpha=1, + edgecolor="black", + ) + fig.savefig( + outputfolder / f"Histogram_observation_epochs_{fire_project}.png", + dpi=400, + ) + plt.close(fig) diff --git a/src/fire/api/geodetic_levelling/metric_to_gpu_transformation.py b/src/fire/api/geodetic_levelling/metric_to_gpu_transformation.py new file mode 100644 index 00000000..0db40295 --- /dev/null +++ b/src/fire/api/geodetic_levelling/metric_to_gpu_transformation.py @@ -0,0 +1,690 @@ +"""This module contains functions for transformation of metric heights or height differences +to geopotential units or vice versa. +""" + +from math import sin, pi, isnan +from pathlib import Path + +import pandas as pd +import pyproj + +from fire.api.geodetic_levelling.tidal_transformation import ( + transform_gravity_from_tidal_system_to_tidal_system, +) + +import fire.api.geodetic_levelling.geophysical_parameters as geo_p + + +def interpolate_gravity( + latitude: float, + longitude: float, + grid_inputfolder: Path, + gravitymodel: str, +) -> float: + """Interpolate in gravity model. + + Interpolates bilinearly in a grid-based gravity model and returns the result as a float. + + Args: + latitude: float, latitude for which gravity is interpolated, in units of degrees + longitude: float, longitude for which gravity is interpolated, in units of degrees + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, grid-based model providing gravity in units of mGal (1 mGal = 10^-5 m/s^2), + must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + float, interpolated gravity in units of m/s^2 + + Raises: + ? + + TO DO: Lav seperat funktion create_pyproj_transformer, således at Transformer-objektet + ikke oprettes hver gang der interpoleres i grid-modellen? transformer: pyproj.Transformer + """ + pyproj.datadir.append_data_dir(grid_inputfolder) + + # Transformer object for interpolation in gravity model + transformer = pyproj.Transformer.from_pipeline( + f"+proj=vgridshift +grids={gravitymodel}" + ) + + # Interpolated gravity in units of m/s^2 + gravity = transformer.transform(longitude, latitude, 0)[2] * 1e-5 * -1 + + return gravity + + +def convert_metric_height_diff_to_geopotential_height_diff( + height_diff: float, + point_from_lat: float, + point_from_long: float, + point_to_lat: float, + point_to_long: float, + tidal_system: str | None, + grid_inputfolder: Path, + gravitymodel: str, +) -> tuple[float, float]: + """Convert a metric height difference to a geopotential height difference. + + Converts a metric height difference to a geopotential height difference (in units of gpu) + and returns the converted height difference and the m2gpu multiplication factor in a tuple. + + The gravity model used for the conversion is assumed to be in zero tide system as this is + the conventional tide system for gravity. + + If the input height difference is in the zero tide system, the gravity interpolated from the + gravity model is not tidally transformed. + + If the input height difference is in non-tidal or mean tide system, the gravity interpolated + from the gravity model is transformed from the zero tide system to the tidal system of the + input height difference. + + If the input height difference is not corrected for tidal effects, the gravity interpolated + from the gravity model is transformed from the zero tide system to the mean tide system. + + Args: + height_diff: float, metric height difference to be converted to gpu + point_from_lat: float, latitude of from point in units of degrees + point_from_long: float, longitude of from point in units of degrees + point_to_lat: float, latitiude of to point in units of degrees + point_to_long: float, longitude of to point in units of degrees + tidal_system: str|None, tidal system of input height difference, i.e. "non", "mean" or "zero" + for non-tidal, mean tide or zero tide or None if the input height difference is not corrected + for tidal effects + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, gravity model used for the conversion of a height difference to gpu, + must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + tuple[float, float], a tuple containing the converted height difference + in units of gpu (1 gpu = 10 m^2/s^2) and the m2gpu multiplication factor in units of m/s^2 + + Raises: + ? + """ + # Point_from and point_to gravity in units of m/s^2 + point_from_gravity = interpolate_gravity( + point_from_lat, + point_from_long, + grid_inputfolder, + gravitymodel, + ) + + point_to_gravity = interpolate_gravity( + point_to_lat, + point_to_long, + grid_inputfolder, + gravitymodel, + ) + + # Interpolated gravity is tidally transformed if tidal system of input height difference + # is different than zero tide + if tidal_system == "zero": + pass + + elif tidal_system == "non": + point_from_gravity = transform_gravity_from_tidal_system_to_tidal_system( + point_from_gravity, point_from_lat, "zero_to_non" + ) + + point_to_gravity = transform_gravity_from_tidal_system_to_tidal_system( + point_to_gravity, point_to_lat, "zero_to_non" + ) + + elif tidal_system == "mean" or tidal_system == None: + point_from_gravity = transform_gravity_from_tidal_system_to_tidal_system( + point_from_gravity, point_from_lat, "zero_to_mean" + ) + + point_to_gravity = transform_gravity_from_tidal_system_to_tidal_system( + point_to_gravity, point_to_lat, "zero_to_mean" + ) + + # Mean gravity in units of m/s^2 + mean_gravity = (point_from_gravity + point_to_gravity) / 2 + + # Conversion of height_diff to geopotential units (1 gpu = 10 m^2/s^2) + m2gpu_factor = mean_gravity * 0.1 + height_diff = height_diff * m2gpu_factor + + return (height_diff, m2gpu_factor) + + +def calculate_normal_gravity( + latitude: float, +) -> float: + """Calculate normal gravity at the GRS80 ellipsoid. + + Calculates normal gravity at the GRS80 ellipsoid. + + Reference: + Johannes Ihde et al., Conventions for the Definition and Realization of a + European Vertical Reference System (EVRS) - EVRS Conventions 2007, p. 10, eq. (A-1). + EUREF, 2019 + + H. Moritz, GEODETIC REFERENCE SYSTEM 1980 + + TO DO: Tidal system of calculated normal gravity? + + Args: + latitude: float, latitude for which normal gravity is calculated, in units of degrees + + Returns: + float, calculated normal gravity in units of m/s^2 + + Raises: + ? + """ + # Conversion of latitude to radians + latitude = (latitude / 360) * 2 * pi + + # Coefficients of series expansion for calculation of normal gravity + a = 0.0052790414 + b = 0.0000232718 + c = 0.0000001262 + d = 0.0000000007 + + normal_gravity = geo_p.normal_gravity_equator_GRS80 * ( + 1 + + a * sin(latitude) ** 2 + + b * sin(latitude) ** 4 + + c * sin(latitude) ** 6 + + d * sin(latitude) ** 8 + ) + + return normal_gravity + + +def calculate_average_normal_gravity( + latitude: float, + normal_height: float, +) -> float: + """Calculate average normal gravity. + + Calculates the average normal gravity along the normal plumb line between + the GRS80 ellipsoid and the telluroid. + + Reference: + Johannes Ihde et al., Conventions for the Definition and Realization of a + European Vertical Reference System (EVRS) - EVRS Conventions 2007, p. 10, eq. (A-2). + EUREF, 2019 + + H. Moritz, GEODETIC REFERENCE SYSTEM 1980 + + Args: + latitude: float, latitude for which average normal gravity is calculated, in units of degrees + normal_height: float, approximate normal height, in units of meters + + Returns: + float, calculated average normal gravity in units of m/s^2 + + Raises: + ? + + TO DO: Handling of tidal system?, Tidal system of calculated average normal gravity? + """ + # Calculation of normal gravity at the ellipsoid + normal_gravity = calculate_normal_gravity(latitude) + + # Conversion of latitude to radians + latitude = (latitude / 360) * 2 * pi + + # Calculation of average normal gravity + r = 1 + geo_p.f_GRS80 + geo_p.m_GRS80 - 2 * geo_p.f_GRS80 * sin(latitude) ** 2 + s = normal_height / geo_p.a_GRS80 + + average_normal_gravity = normal_gravity * (1 - r * s + s**2) + + return average_normal_gravity + + +def convert_geopotential_height_to_normal_height( + height: float, + latitude: float, + conversion: str, + approx_normal_height: float = 0, + iterate: bool = True, +) -> tuple[float, float]: + """Convert a geopotential height to normal height or vice versa. + + Converts a geopotential height to GRS80 normal height or vice versa. + + References: + Johannes Ihde et al., Conventions for the Definition and Realization of a + European Vertical Reference System (EVRS) - EVRS Conventions 2007. EUREF, 2019 + + H. Moritz, GEODETIC REFERENCE SYSTEM 1980 + + If a geopotential height is to be converted to normal height, an a priori + normal height is required. If no argument is passed for the parameter approx_normal_height, + a default value of zero will be used as a priori normal height. By default, the + normal height is calculated iteratively until the difference between the current and previous + iteration step is less than 0.01 mm. If the normal height is calculated in only one step using + the default a priori normal height, an error of several millimeters can occur. Therefore, + it is recommended to calculate the normal height iteratively if no a priori normal height is + passed. + + TO DO: Handling of tidal corrections/systems? + TO DO: Change the condition for iteration to a looser one, based on the formulae for normal height + + Args: + height: float, input/source height to be converted, in units of gpu (1 gpu = 10 m^2/s^2) if a + geopotential height or in units of m if a normal height + latitude: float, latitude of input/source height, in units of degrees + conversion: str, specification of source and target height, "geopot_to_normal" or + "normal_to_geopot" + approx_normal_height: float = 0, optional parameter, approx normal height in units of m, + only relevant if a geopotential height is to be converted to normal height, default value is 0 + iterate: bool = True, optional parameter, determines whether or not the output/target + normal height is calculated iteratively, default value is True + + Returns: + tuple[float, float], a tuple containing the converted height (in units of gpu + (1 gpu = 10 m^2/s^2) if a geopotential height or in units of m if a normal height) + and the average normal gravity (in units of 10 m/s^2) + + Raises: + ? + """ + # Conversion of a geopotential height to normal height + if conversion == "geopot_to_normal": + # Calculation of average normal gravity in units of 10 m/s^2 + average_normal_gravity = ( + calculate_average_normal_gravity(latitude, approx_normal_height) * 0.1 + ) + + # Normal height in units of meters + height_converted = height / average_normal_gravity + + # Iterative calculation of normal height + # The iteration is started if height_converted is not nan + if iterate == True and isnan(height_converted) == False: + while not ( + -0.00001 <= (height_converted - approx_normal_height) <= 0.00001 + ): + approx_normal_height = height_converted + + # Calculation of average normal gravity in units of 10 m/s^2 + average_normal_gravity = ( + calculate_average_normal_gravity(latitude, approx_normal_height) + * 0.1 + ) + + # Normal height in units of meters + height_converted = height / average_normal_gravity + + # Conversion of a normal height to geopotential height + elif conversion == "normal_to_geopot": + # Calculation of average normal gravity in units of 10 m/s^2 + average_normal_gravity = ( + calculate_average_normal_gravity(latitude, height) * 0.1 + ) + + # Geopotential height in units of gpu (1 gpu = 10 m^2/s^2) + height_converted = height * average_normal_gravity + + else: + exit( + "Function convert_geopotential_height_to_normal_height: Wrong argument for\n\ + parameter conversion." + ) + + return (height_converted, average_normal_gravity) + + +def convert_geopotential_height_to_helmert_height( + height: float, + latitude: float, + longitude: float, + grid_inputfolder: Path, + gravitymodel: str, + conversion: str, + tidal_system: str = None, + approx_helmert_height: float = 0, + iterate: bool = True, +) -> tuple[float, float]: + """Convert a geopotential height to Helmert height or vice versa. + + Converts a geopotential height to Helmert height or vice versa. + + Reference: + Klaus Schmidt, The Danish height system DVR90, pp. app 14-15. + National Survey and Cadastre, 2000 + + If a geopotential height is to be converted to Helmert height, an a priori + Helmert height is required. If no argument is passed for the parameter approx_helmert_height, + a default value of zero will be used as a priori Helmert height. By default, the + Helmert height is calculated iteratively until the difference between the current and previous + iteration step is less than 0.01 mm. If the Helmert height is calculated in only one step using + the default a priori Helmert height, an error of several millimeters can occur. Therefore, + it is recommended to calculate the Helmert height iteratively if no a priori Helmert height is + passed. + + The gravity model used for the conversion of a geopotential height to Helmert height + (or vice versa) is assumed to be in zero tide system as this is the conventional tide system + for gravity. + + If the input height is given in the zero tide system, the gravity value interpolated from the + gravity model is not tidally transformed. + + If the input height is given in the non-tidal or mean tide system, the gravity value interpolated + from the gravity model is transformed from the zero tide system to the tidal system of the + input height. + + If the input height is not corrected for tidal effects, the gravity value interpolated + from the gravity model is transformed from the zero tide system to the mean tide system. + + TO DO: Change the condition for iteration to a looser one, based on the formulae for Helmert height + + Args: + height: float, input/source height to be converted, in units of gpu (1 gpu = 10 m^2/s^2) if a + geopotential height or in units of m if a Helmert height + latitude: float, latitude of input/source height, in units of degrees + longitude: float, longitude of input/source height, in units of degrees + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, gravity model used for the height conversion, must be in GeoTIFF + or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + conversion: str, specification of source and target height, "geopot_to_helmert" or + "helmert_to_geopot" + tidal_system: str = None, optional parameter, tidal system of input height, i.e. "non", "mean" + or "zero" for non-tidal, mean tide or zero tide. If no argument is passed it is assumed that + the input height is not corrected for tidal effects + approx_helmert_height: float = 0, optional parameter, approx Helmert height in units of m, + only relevant if a geopotential height is to be converted to Helmert height, default value is 0 + iterate: bool = True, optional parameter, determines whether or not the output/target + Helmert height is calculated iteratively, default value is True + + Returns: + tuple[float, float], a tuple containing the converted height (in units of gpu + (1 gpu = 10 m^2/s^2) if a geopotential height or in units of m if a Helmert height) + and the conversion factor (Helmert height to geopotential height) in units of 10 m/s^2 + + Raises: + ? + """ + # Interpolated gravity in units of m/s^2 + gravity = interpolate_gravity( + latitude, + longitude, + grid_inputfolder, + gravitymodel, + ) + + # Interpolated gravity is tidally transformed if tidal system of input height + # is different than zero tide + if tidal_system == "zero": + pass + + elif tidal_system == "non": + gravity = transform_gravity_from_tidal_system_to_tidal_system( + gravity, latitude, "zero_to_non" + ) + + elif tidal_system == "mean" or tidal_system == None: + gravity = transform_gravity_from_tidal_system_to_tidal_system( + gravity, latitude, "zero_to_mean" + ) + + # Conversion of a geopotential height to Helmert height + if conversion == "geopot_to_helmert": + # Conversion factor (metric Helmert height to geopotential height) in units of 10 m/s^2 + conversion_factor = (gravity * 0.1) + (0.07045 * 1e-6 * approx_helmert_height) + + # Helmert height in units of meters + height_converted = height / conversion_factor + + # Iterative calculation of Helmert height + # The iteration is started if height_converted is not nan + if iterate == True and isnan(height_converted) == False: + while not ( + -0.00001 <= (height_converted - approx_helmert_height) <= 0.00001 + ): + approx_helmert_height = height_converted + + # Conversion factor (metric Helmert height to geopotential height) in units of 10 m/s^2 + conversion_factor = (gravity * 0.1) + ( + 0.07045 * 1e-6 * approx_helmert_height + ) + + # Helmert height in units of meters + height_converted = height / conversion_factor + + # Conversion of a Helmert height to geopotential height + elif conversion == "helmert_to_geopot": + # Conversion factor (metric Helmert height to geopotential height) in units of 10 m/s^2 + conversion_factor = (gravity * 0.1) + (0.07045 * 1e-6 * height) + + # Geopotential height in units of gpu (1 gpu = 10 m^2/s^2) + height_converted = height * conversion_factor + + else: + exit( + "Function convert_geopotential_height_to_helmert_height: Wrong argument for\n\ + parameter conversion." + ) + + return (height_converted, conversion_factor) + + +def convert_geopotential_heights_to_metric_heights( + points_df: pd.DataFrame, + # fire_project: str, + # excel_inputfolder: Path, + # outputfolder: Path, + conversion: str, + grid_inputfolder: Path = None, + gravitymodel: str = None, + tidal_system: str = None, + iterate: bool = True, +) -> None: + """Convert geopotential heights to metric heights or vice versa. + + Converts geopotential heights of a FIRE project to metric heights or vice versa. + + If geopotential heights are to be converted to metric heights (Helmert heights or + normal heights) the input heights are taken from column "Ny kote" in the sheet + "Kontrolberegning" in the input excel-file and the converted heights are written to + column "Ny kote" in the sheet "Kontrolberegning" in the output excel-file. + + The conversion of geopotential heights to metric heights (Helmert heights or + normal heights) requires a priori metric heights, which are taken from column "Kote" + in the sheet "Kontrolberegning" in the input excel-file. + + If Helmert heights or normal heights are to be converted to geopotential heights + the input heights are taken from column "Kote" in the sheet "Kontrolberegning" in the + input excel-file and the converted heights are written to column "Ny kote" in the + sheet "Kontrolberegning" in the output excel-file. + + Args: + fire_project: str, name of FIRE project with heights to be converted, must be in accordance + with the name of the input excel-file, e.g. "asmei_temp" + excel_inputfolder: Path, folder with input FIRE project/excel-file with heights to be converted + outputfolder: Path, folder for output FIRE project/excel-file with converted heights + conversion: str, specification of source and target height, "geopot_to_helmert", + "helmert_to_geopot", "geopot_to_normal" or "normal_to_geopot" + grid_inputfolder: Path = None, optional parameter, folder for input grid, i.e. gravity model, + only relevant if geopotential heights are to be converted to Helmert heights or vice versa + gravitymodel: str = None, optional parameter, gravity model used for the conversion of heights, + must be in GeoTIFF or GTX file format, only relevant if geopotential heights are to be + converted to Helmert heights or vice versa + tidal_system: str = None, optional parameter, tidal system of input heights, i.e. "non", "mean" + or "zero" for non-tidal, mean tide or zero tide. If no argument is passed it is assumed that + the input heights are not corrected for tidal effects, only relevant if geopotential heights + are to be converted to Helmert heights or vice versa + iterate: bool = True, optional parameter, determines whether or not output/target + metric heights are calculated iteratively, default value is True + + Returns: + None + + Raises: + ? Hvis grid_inputfolder ikke findes, hvis grid-fil ikke findes, hvis input excel-fil ikke findes + + Input file: + FIRE project/excel-file with heights to be converted, e.g. "asmei_temp.xlsx" + + Output file: + Excel-file with converted heights. This file contains the converted heights in column "Ny kote" + as well as the conversion factors or average normal gravity values used for height conversion. + Except for that the file is identical to the input excel-file. + + TO DO: Håndtering manglende a priori værdi? + TO DO: Håndtering manglende inputhøjde? + """ + # # Make sure that the output folder exists + # outputfolder.mkdir(parents=True, exist_ok=True) + + # excel_inputfile = excel_inputfolder / f"{fire_project}.xlsx" + + # # DataFrame with heights etc. from input fire project + # points_df = pd.read_excel(excel_inputfile, sheet_name="Kontrolberegning") + + if conversion == "geopot_to_normal": + for index in points_df.index: + h_adjusted = points_df.at[index, "Ny kote"] + h_db = points_df.at[index, "Kote"] + latitude = points_df.at[index, "Nord"] + + (height_converted, average_normal_gravity) = ( + convert_geopotential_height_to_normal_height( + h_adjusted, + latitude, + "geopot_to_normal", + approx_normal_height=h_db, + iterate=iterate, + ) + ) + points_df.at[index, "Ny kote"] = height_converted + # points_df.at[index, "Average normal gravity [10 m/s^2]"] = ( + # average_normal_gravity + # ) + + elif conversion == "normal_to_geopot": + for index in points_df.index: + h_db = points_df.at[index, "Kote"] + latitude = points_df.at[index, "Nord"] + + (height_converted, average_normal_gravity) = ( + convert_geopotential_height_to_normal_height( + h_db, + latitude, + "normal_to_geopot", + ) + ) + points_df.at[index, "Kote"] = height_converted + # points_df.at[index, "Average normal gravity [10 m/s^2]"] = ( + # average_normal_gravity + # ) + + elif ( + conversion == "geopot_to_helmert" + and (grid_inputfolder is not None) + and (gravitymodel is not None) + ): + for index in points_df.index: + h_adjusted = points_df.at[index, "Ny kote"] + h_db = points_df.at[index, "Kote"] + latitude = points_df.at[index, "Nord"] + longitude = points_df.at[index, "Øst"] + + (height_converted, conversion_factor) = ( + convert_geopotential_height_to_helmert_height( + h_adjusted, + latitude, + longitude, + grid_inputfolder, + gravitymodel, + "geopot_to_helmert", + tidal_system=tidal_system, + approx_helmert_height=h_db, + iterate=iterate, + ) + ) + points_df.at[index, "Ny kote"] = height_converted + # points_df.at[index, "Conversion factor [10 m/s^2]"] = conversion_factor + + elif ( + conversion == "helmert_to_geopot" + and (grid_inputfolder is not None) + and (gravitymodel is not None) + ): + for index in points_df.index: + h_db = points_df.at[index, "Kote"] + latitude = points_df.at[index, "Nord"] + longitude = points_df.at[index, "Øst"] + + (height_converted, conversion_factor) = ( + convert_geopotential_height_to_helmert_height( + h_db, + latitude, + longitude, + grid_inputfolder, + gravitymodel, + "helmert_to_geopot", + tidal_system=tidal_system, + ) + ) + points_df.at[index, "Kote"] = height_converted + # points_df.at[index, "Conversion factor [10 m/s^2]"] = conversion_factor + + else: + exit( + "Function convert_geopotential_heights_to_metric_heights: Wrong arguments for\n\ + parameter conversion and/or grid_inputfolder and/or gravitymodel." + ) + + return points_df + + # DataFrame with parameters of output fire project + parameters_df = pd.read_excel(excel_inputfile, sheet_name="Parametre") + + if conversion == "geopot_to_normal" or conversion == "normal_to_geopot": + parameters_new_df = pd.DataFrame( + { + "Navn": [ + "Conversion of heights", + ], + "Værdi": [conversion], + }, + ) + + if conversion == "geopot_to_helmert" or conversion == "helmert_to_geopot": + parameters_new_df = pd.DataFrame( + { + "Navn": [ + "Conversion of heights", + "Gravitymodel for conversion of heights", + ], + "Værdi": [conversion, gravitymodel], + }, + ) + + parameters_df = pd.concat([parameters_df, parameters_new_df], ignore_index=True) + + # Generation of output fire project/excel file with converted heights + with pd.ExcelWriter( + outputfolder / f"{fire_project}.xlsx" + ) as writer: # pylint: disable=abstract-class-instantiated + pd.read_excel(excel_inputfile, sheet_name="Projektforside").to_excel( + writer, "Projektforside", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Sagsgang").to_excel( + writer, "Sagsgang", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Nyetablerede punkter").to_excel( + writer, "Nyetablerede punkter", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Notater").to_excel( + writer, "Notater", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Filoversigt").to_excel( + writer, "Filoversigt", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Observationer").to_excel( + writer, "Observationer", index=False + ) + pd.read_excel(excel_inputfile, sheet_name="Punktoversigt").to_excel( + writer, "Punktoversigt", index=False + ) + points_df.to_excel(writer, "Kontrolberegning", index=False) + parameters_df.to_excel(writer, "Parametre", index=False) diff --git a/src/fire/api/geodetic_levelling/tidal_transformation.py b/src/fire/api/geodetic_levelling/tidal_transformation.py new file mode 100644 index 00000000..6a995273 --- /dev/null +++ b/src/fire/api/geodetic_levelling/tidal_transformation.py @@ -0,0 +1,541 @@ +"""This module contains functions for tidal correction and transformation of gravity and height etc.""" + +from math import cos, sin, pi +from pathlib import Path + +from astropy import constants as const +from astropy.coordinates import AltAz, solar_system_ephemeris, EarthLocation, get_body +from astropy.time import Time +import astropy.units as u +import pandas as pd +import pyproj + +import fire.api.geodetic_levelling.geophysical_parameters as geo_p + + +def apply_tidal_corrections_to_height_diff( + height_diff: float, + point_from_lat: float, + point_from_long: float, + point_to_lat: float, + point_to_long: float, + epoch_obs: pd.Timestamp, + tidal_system: str, + grid_inputfolder: Path = None, + gravitymodel: str = None, +) -> tuple[float, float]: + """Apply tidal corrections to a metric height difference. + + Applies tidal corrections to a metric height difference and returns the corrected + height difference and the correction itself in a tuple. + + The application of tidal corrections to a metric height difference implies that all periodic + tidal effects are removed, whereas the permanent tidal effects are removed or retained + (in whole or in part) depending on the specified tidal system. + + Reference: + Klaus Schmidt, The Danish height system DVR90, pp. app 15-16. + National Survey and Cadastre, 2000 + + Args: + height_diff: float, metric height difference to be tidally corrected + point_from_lat: float, latitude of from point in units of degrees + point_from_long: float, longitude of from point in units of degrees + point_to_lat: float, latitiude of to point in units of degrees + point_to_long: float, longitude of to point in units of degrees + epoch_obs: pd.Timestamp, epoch/time of observation (format: yyyy-mm-dd hh:mm:ss) + tidal_system: str, tidal system of output height difference, "non", "mean" or "zero" + for non-tidal, mean tide or zero tide + grid_inputfolder: Path = None, optional parameter, folder for input grid, i.e. gravity model, + only relevant if tidal system is mean tide or zero tide + gravitymodel: str = None, optional parameter, grid-based model providing gravity in units of + mGal (1 mGal = 10^-5 m/s^2), must be in GeoTIFF or GTX file format, only relevant if + tidal system is mean tide or zero tide + + Returns: + tuple[float, float], a tuple containing the corrected height difference and + the correction itself in units of meters + + Raises: + ? + + TO DO: Højere grad af konsistens fsva. vinkelmål deg/rad? + TO DO: Split op i mindre underfunktioner + TO DO: Implementer option vedr. brug af approksimative formler til transformation til mean og + zero tide + TO DO: Take daylight saving time into account + TO DO: Handling of small inconsistence regarding tilt factor/diminution coefficient (0.7 vs. 0.68) + """ + # Calculation of levelling section length and azimuth + # Azimuth is calculated clockwise from north; interval of azimuth: [-180 deg; 180 deg] + geod = pyproj.Geod(ellps="GRS80") + azimuth_forward, azimuth_back, section_length = geod.inv( + point_from_long, point_from_lat, point_to_long, point_to_lat + ) + + # Conversion of levelling section azimuth to radians + azimuth_forward = azimuth_forward * 2 * (pi / 360) + + # Observation epoch in UTC + # Daylight saving time is not taken into account + epoch_obs = epoch_obs + pd.Timedelta(hours=-1) + t = Time(epoch_obs, scale="utc") + + # Mean geographic coordinates + mean_lat = (point_from_lat + point_to_lat) / 2 + mean_long = (point_from_long + point_to_long) / 2 + + # Conversion of mean geographic coordinates to cartesian ITRS coordinates + loc = EarthLocation(lat=mean_lat * u.deg, lon=mean_long * u.deg, height=0 * u.m) + + # Apparent positions (ra, dec, dist) of the Moon and Sun in GCRS + with solar_system_ephemeris.set("jpl"): + moon = get_body("moon", t, loc) + sun = get_body("sun", t, loc) + + # Altitude and azimuth of the Moon and Sun without refraction effects + # Interval of altitude: [-pi/2; pi/2] + # Azimuth is calculated clockwise from north; interval of azimuth: [0; 2*pi] + altazframe = AltAz(obstime=t, location=loc, pressure=0) + moon_altaz = moon.transform_to(altazframe) + sun_altaz = sun.transform_to(altazframe) + altitude_moon = moon_altaz.alt.radian + azimuth_moon = moon_altaz.az.radian + altitude_sun = sun_altaz.alt.radian + azimuth_sun = sun_altaz.az.radian + + # Conversion of altitude to zenith distance + zenith_dist_moon = (pi / 2) - altitude_moon + zenith_dist_sun = (pi / 2) - altitude_sun + + # Calculation of tidal corrections due to the Moon and Sun in units of 1e-8 m + tidal_corr_moon = ( + section_length + * 8.5 + * sin(2 * zenith_dist_moon) + * cos(azimuth_moon - azimuth_forward) + ) + tidal_corr_sun = ( + section_length + * 3.9 + * sin(2 * zenith_dist_sun) + * cos(azimuth_sun - azimuth_forward) + ) + # Tidal correction of height_diff (0.7 is diminution coefficient corresponding to a yielding + # Earth). This correction removes all tidals effects caused by the Moon and Sun, i.e. both + # periodic and permanent tidal effects and both direct and indirect effects (caused by + # a yielding/deforming Earth). Consequently, the corrected height difference is in + # non-tidal system + tidal_corr = (tidal_corr_moon + tidal_corr_sun) * 0.7 * 1e-8 + height_diff_corrected = height_diff + tidal_corr + + if tidal_system == "non": + pass + + # If the specified tidal system is mean tide, the corrected height difference is transformed + # from non-tidal to mean tide + elif ( + tidal_system == "mean" + and (gravitymodel is not None) + and (grid_inputfolder is not None) + ): + height_diff_corrected = transform_height_diff_from_tidal_system_to_tidal_system( + height_diff_corrected, + point_from_lat, + point_from_long, + point_to_lat, + point_to_long, + "non_to_mean", + grid_inputfolder, + gravitymodel, + ) + + tidal_corr = height_diff_corrected - height_diff + + # If the specified tidal system is zero tide, the corrected height difference is transformed + # from non-tidal to zero tide + elif ( + tidal_system == "zero" + and (gravitymodel is not None) + and (grid_inputfolder is not None) + ): + height_diff_corrected = transform_height_diff_from_tidal_system_to_tidal_system( + height_diff_corrected, + point_from_lat, + point_from_long, + point_to_lat, + point_to_long, + "non_to_zero", + grid_inputfolder, + gravitymodel, + ) + + tidal_corr = height_diff_corrected - height_diff + + else: + exit( + "Function apply_tidal_corrections_to_height_diff: Wrong arguments for\n\ + parameter tidal_system and/or gravitymodel and/or grid_inputfolder." + ) + + return (height_diff_corrected, tidal_corr) + + +def calculate_perm_tidal_gravitation( + latitude: float, + celestial_body: str, +) -> float: + """Calculate the permanent tidal gravitation assuming a rigid Earth. + + Calculates the permanent tidal gravitation caused by the Moon or the Sun assuming a rigid Earth + and returns the result as a float. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + p. 120, eq. (8). Nordic Geodetic Commission, 1988 + + Args: + latitude: float, latitude for which the permanent tidal gravitation is calculated, + in units of degrees + celestial_body: str, celestial body for which the permanent tidal gravitation is calculated, + "moon" or "sun" + + Returns: + float, the permanent tidal gravitation in units of m/s^2 + + Raises: + ? + """ + # The permanent tidal gravitation is given by -dW/dr, i.e. + # minus the partial derivative of the permanent tidal potential wrt Earth radius + perm_tidal_gravitation = -(2 / geo_p.radius_earth) * calculate_perm_tidal_potential( + latitude, celestial_body + ) + + return perm_tidal_gravitation + + +def calculate_perm_tidal_potential( + latitude: float, + celestial_body: str, +) -> float: + """Calculate the permanent tidal potential assuming a rigid Earth. + + Calculates the permanent tidal potential caused by the Moon or the Sun assuming a rigid Earth + and returns the result as a float. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + p. 119, eq. (5). Nordic Geodetic Commission, 1988 + + Args: + latitude: float, latitude for which the permanent tidal potential is calculated, + in units of degrees + celestial_body: str, celestial body for which the permanent tidal potential is calculated, + "moon" or "sun" + + Returns: + float, the permanent tidal potential in units of m^2/s^2 + + Raises: + ? + """ + # Conversion of latitude to radians + latitude = (latitude / 360) * 2 * pi + + # Inclination dependent term + inclination_term = (3 / 2) * (sin(geo_p.epsilon)) ** 2 - 1 + + # Latitude dependent term + latitude_term = 3 * (sin(latitude)) ** 2 - 1 + + if celestial_body == "moon": + + perm_tidal_potential = ( + ( + (const.G.value * geo_p.moon_mass * geo_p.radius_earth**2) + / (4 * geo_p.moon_dist**3) + ) + * inclination_term + * latitude_term + ) + + if celestial_body == "sun": + + perm_tidal_potential = ( + ( + (const.G.value * const.M_sun.value * geo_p.radius_earth**2) + / (4 * const.au.value**3) + ) + * inclination_term + * latitude_term + ) + + return perm_tidal_potential + + +def calculate_perm_tidal_deformation_geoid( + latitude: float, + longitude: float, + celestial_body: str, + grid_inputfolder: Path, + gravitymodel: str, +) -> float: + """Calculate the permanent tidal deformation of the geoid assuming a rigid Earth. + + Calculates the permanent tidal deformation of the geoid caused by the Moon or the Sun + assuming a rigid Earth and returns the result as a float. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + p. 120, eq. (6). Nordic Geodetic Commission, 1988 + + Args: + latitude: float, latitude for which the permanent tidal deformation of the geoid is calculated, + in units of degrees + longitude: float, longitude for which the permanent tidal deformation of the geoid is calculated, + in units of degrees + celestial_body: str, celestial body for which the permanent tidal deformation of the geoid + is calculated, "moon" or "sun" + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, grid-based model providing gravity in units of mGal (1 mGal = 10^-5 m/s^2), + must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + float, the permanent tidal deformation of the geoid in units of m + + Raises: + ? + + NB:Handling of tidal systems? Bør der tages højde for at input gravity er zero tide? + Giver det mening at bruge tyngden ved jordoverfladen? + Burde det være tyngden på geoiden? + """ + gravity = interpolate_gravity(latitude, longitude, grid_inputfolder, gravitymodel) + + # Permanent tidal potential + perm_tidal_potential = calculate_perm_tidal_potential(latitude, celestial_body) + + # Permanent tidal deformation of the geoid + perm_tidal_deformation_geoid = perm_tidal_potential / gravity + + return perm_tidal_deformation_geoid + + +def transform_gravity_from_tidal_system_to_tidal_system( + gravity: float, + latitude: float, + transformation: str, +) -> float: + """Transform gravity from one tidal system to another tidal system. + + Transforms gravity from one tidal system to another tidal system and returns the + result as a float. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + pp. 124-128, eq. (19), (20), (22). Nordic Geodetic Commission, 1988 + + Args: + gravity: float, gravity to be transformed from one tidal system to another, in units of m/s^2 + latitude: float, latitude at which the input gravity is measured, in units of degrees + transformation: str, specification of source and target tidal system, e.g. "non_to_mean" + + Returns: + float, the transformed gravity in units of m/s^2 + + Raises: + ? + """ + # The permanent tidal gravitation assuming a rigid Earth in units of m/s^2 + perm_tidal_gravitation = (calculate_perm_tidal_gravitation(latitude, "moon")) + ( + calculate_perm_tidal_gravitation(latitude, "sun") + ) + + # The permanent tidal gravitation for a deformable Earth in units of m/s^2 + perm_tidal_gravitation_deform = geo_p.delta * perm_tidal_gravitation + + if transformation == "non_to_mean": + gravity_transformed = gravity + perm_tidal_gravitation_deform + + elif transformation == "non_to_zero": + gravity_transformed = ( + gravity + perm_tidal_gravitation_deform - perm_tidal_gravitation + ) + + elif transformation == "mean_to_non": + gravity_transformed = gravity - perm_tidal_gravitation_deform + + elif transformation == "mean_to_zero": + gravity_transformed = gravity - perm_tidal_gravitation + + elif transformation == "zero_to_non": + gravity_transformed = ( + gravity - perm_tidal_gravitation_deform + perm_tidal_gravitation + ) + + elif transformation == "zero_to_mean": + gravity_transformed = gravity + perm_tidal_gravitation + + return gravity_transformed + + +def transform_height_from_tidal_system_to_tidal_system( + height: float, + latitude: float, + longitude: float, + transformation: str, + grid_inputfolder: Path, + gravitymodel: str, +) -> float: + """Transform a geophysical height from one tidal system to another tidal system. + + Transforms a geophysical height above the geoid (e.g. a levelled height) from one tidal system + to another tidal system and returns the result as a float. + + The height to be transformed is assumed to have been tidally corrected + (i.e. referred to a specific tidal system) before being transformed to another tidal system + using this function. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + p. 128-129, eq. (24), (25). Nordic Geodetic Commission, 1988 + + Args: + height: float, height to be transformed from one tidal system to another, in units of m + latitude: float, latitude of input height, in units of degrees + longitude: float, longitude of input height, in units of degrees + transformation: str, specification of source and target tidal system, e.g. "non_to_mean" + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, grid-based model providing gravity in units of mGal (1 mGal = 10^-5 m/s^2), + must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + float, the transformed height in units of m + + Raises: + ? + + NB: Hvilken rolle spiller det om højden er Helmert-højde, normalhøjde? Nok ingen? + """ + # The permanent tidal deformation of the geoid assuming a rigid Earth in units of m + perm_tidal_deformation_geoid_moon = calculate_perm_tidal_deformation_geoid( + latitude, + longitude, + "moon", + grid_inputfolder, + gravitymodel, + ) + + perm_tidal_deformation_geoid_sun = calculate_perm_tidal_deformation_geoid( + latitude, + longitude, + "sun", + grid_inputfolder, + gravitymodel, + ) + + perm_tidal_deformation_geoid = ( + perm_tidal_deformation_geoid_moon + perm_tidal_deformation_geoid_sun + ) + + # Nedentående skal verificeres + # The permanent tidal deformation of the geoid for a deformable Earth in units of m + perm_tidal_deformation_geoid_deform = geo_p.delta * perm_tidal_deformation_geoid + + if transformation == "non_to_mean": + height_transformed = height + perm_tidal_deformation_geoid_deform + + elif transformation == "non_to_zero": + height_transformed = ( + height + perm_tidal_deformation_geoid_deform - perm_tidal_deformation_geoid + ) + + elif transformation == "mean_to_non": + height_transformed = height - perm_tidal_deformation_geoid_deform + + elif transformation == "mean_to_zero": + height_transformed = height - perm_tidal_deformation_geoid + + elif transformation == "zero_to_non": + height_transformed = ( + height - perm_tidal_deformation_geoid_deform + perm_tidal_deformation_geoid + ) + + elif transformation == "zero_to_mean": + height_transformed = height + perm_tidal_deformation_geoid + + return height_transformed + + +def transform_height_diff_from_tidal_system_to_tidal_system( + height_diff: float, + point_from_lat: float, + point_from_long: float, + point_to_lat: float, + point_to_long: float, + transformation: str, + grid_inputfolder: Path, + gravitymodel: str, +) -> float: + """Transform a geophysical height difference from one tidal system to another tidal system. + + Transforms a geophysical height difference above the geoid (e.g. a levelled height) from one + tidal system to another tidal system and returns the result as a float. + + The height difference to be transformed is assumed to have been tidally corrected + (i.e. referred to a specific tidal system) before being transformed to another tidal system + using this function. + + Reference: + Martin Ekman, The impact of geodynamic phenomena on systems for height and gravity, + pp. 128-129, eq. (24), (25). Nordic Geodetic Commission, 1988 + + Args: + height_diff: float, geophysical height difference to be transformed from one tidal system + to another, in units of m + point_from_lat: float, latitude of from point in units of degrees + point_from_long: float, longitude of from point in units of degrees + point_to_lat: float, latitiude of to point in units of degrees + point_to_long: float, longitude of to point in units of degrees + transformation: str, specification of source and target tidal system, e.g. "non_to_mean" + grid_inputfolder: Path, folder for input grid, i.e. gravity model + gravitymodel: str, grid-based model providing gravity in units of mGal (1 mGal = 10^-5 m/s^2), + must be in GeoTIFF or GTX file format, e.g. "dk-g-direkte-fra-gri-thokn.tif" + + Returns: + float, the transformed height difference in units of m + + Raises: + ? + + NB: Hvor meget betyder højden af højdeforskellens fra-/til-punkter? Nok ikke ret meget + eftersom de numeriske/approksimative formler ikke tager højde herfor + """ + # Transformation of height difference + height_diff_transformed = ( + height_diff + + transform_height_from_tidal_system_to_tidal_system( + 0, + point_to_lat, + point_to_long, + transformation, + grid_inputfolder, + gravitymodel, + ) + - transform_height_from_tidal_system_to_tidal_system( + 0, + point_from_lat, + point_from_long, + transformation, + grid_inputfolder, + gravitymodel, + ) + ) + + return height_diff_transformed + + +from fire.api.geodetic_levelling.metric_to_gpu_transformation import ( + interpolate_gravity, +) diff --git a/src/fire/api/geodetic_levelling/time_propagation.py b/src/fire/api/geodetic_levelling/time_propagation.py new file mode 100644 index 00000000..61fb2afc --- /dev/null +++ b/src/fire/api/geodetic_levelling/time_propagation.py @@ -0,0 +1,68 @@ +"""This module contains functions for time propagation of height differences.""" + +from pathlib import Path + +import pandas as pd +import pyproj + + +def propagate_height_diff_from_epoch_to_epoch( + height_diff: float, + point_from_lat: float, + point_from_long: float, + point_to_lat: float, + point_to_long: float, + epoch_source: pd.Timestamp, + epoch_target: pd.Timestamp, + grid_inputfolder: Path, + deformationmodel: str, +) -> tuple[float, float]: + """Propagate a metric height difference from one epoch to another. + + Propagates a metric height difference from one epoch to another and returns the propagated + height difference and the correction itself in a tuple. + + Args: + height_diff: float, metric height difference to be propagated + point_from_lat: float, latitude of from point in units of degrees + point_from_long: float, longitude of from point in units of degrees + point_to_lat: float, latitiude of to point in units of degrees + point_to_long: float, longitude of to point in units of degrees + epoch_source: pd.Timestamp, source epoch, e.g. epoch of observation (format: yyyy-mm-dd hh:mm:ss) + epoch_target: pd.Timestamp, target epoch (format: yyyy-mm-dd hh:mm:ss) + grid_inputfolder: Path, folder for input grid, i.e. deformation model + deformationmodel: str, deformation model used for the propagation of height differences, + must be in GeoTIFF or GTX file format, e.g. "NKG2016_lev.tif" + + Returns: + tuple[float, float], a tuple containing the propagated height difference and + the correction itself in units of meters + + Raises: + ? + """ + pyproj.datadir.append_data_dir(grid_inputfolder) + + # Up velocities of point_from and point_to in units of mm/yr + # Flyt til overfunktionen height_diff_corr? + transformer = pyproj.Transformer.from_pipeline( + f"+proj=vgridshift +grids={deformationmodel}" + ) + + point_from_velocity_up = ( + transformer.transform(point_from_long, point_from_lat, 0)[2] * -1 + ) + + point_to_velocity_up = transformer.transform(point_to_long, point_to_lat, 0)[2] * -1 + + # Difference in up velocities in units of mm/yr + velocity_up_diff = point_to_velocity_up - point_from_velocity_up + + # Difference in epochs in units of years + epoch_diff = (epoch_target - epoch_source).days / 365.2425 + + # Propagation of height_diff to epoch_target + epoch_corr = velocity_up_diff * epoch_diff * 0.001 + height_diff = height_diff + epoch_corr + + return (height_diff, epoch_corr) diff --git a/src/fire/api/niv/regnemotor.py b/src/fire/api/niv/regnemotor.py index f37c3bf8..f366ef0d 100644 --- a/src/fire/api/niv/regnemotor.py +++ b/src/fire/api/niv/regnemotor.py @@ -15,6 +15,13 @@ import pandas as pd +from fire.api.geodetic_levelling.geodetic_correction_levelling_obs import ( + apply_geodetic_corrections_to_height_diffs, +) +from fire.api.geodetic_levelling.metric_to_gpu_transformation import ( + convert_geopotential_heights_to_metric_heights +) + # Smarte type hints PunktNavn = str # NivSubnet er en forsimplet udgave af et rigtigt net, som bare indeholder navnene på punkter som indgår @@ -240,6 +247,7 @@ def til_dataframe(self) -> pd.DataFrame: return df_out + # KREBSLW: kan være den her ikke skal være cached_property når vi nu laver konvertering af koterne @cached_property def fastholdte(self) -> dict[PunktNavn, float]: """Find fastholdte punkter og koter til en beregning""" @@ -385,6 +393,11 @@ def filer(self) -> list: """En liste af filnavne som motoren producerer""" pass + @property + @abstractmethod + def parametre(self) -> dict: + """En liste af relevante parametre som motoren bruger""" + pass class GamaRegn(RegneMotor): """ @@ -417,6 +430,11 @@ def filer(self, nye_filnavne): """Sæt nye filnavne""" self.xml_in, self.xml_out, self.html_out = nye_filnavne + @property + def parametre(self) -> dict: + """Beret om diverse parametre brugt i gama-local""" + return dict() + def skriv_gama_inputfil(self): """ Skriv gama-inputfil i XML-format @@ -546,9 +564,252 @@ def udjævn(self): @property def filer(self) -> list: """En liste af filer som DumRegn producerer""" - return None + return [] + + @property + def parametre(self) -> dict: + """En dict af parametre brugt i DumRegn""" + return dict() + + +class GeodætiskRegn(GamaRegn): + """ + En geodætisk regnemotor + + Her ville være et godt sted at dokumentere hvad GeodætiskRegn kan, i stil med det som + er lavet for den overordnede RegneMotor. + + """ + def __init__( + self, + tidal_system: str, + epoch_target: pd.Timestamp, + height_diff_unit: str, + output_height: str, + deformationmodel: str, + gravitymodel: str, + grid_inputfolder: Path, + **kwargs, + ): + # intitialiser parametre + self.tidal_system = tidal_system + self.epoch_target = epoch_target + self.height_diff_unit = height_diff_unit + self.output_height = output_height + self.deformationmodel = deformationmodel + self.gravitymodel = gravitymodel + self.grid_inputfolder = grid_inputfolder + + # do some stuff + + # start GamaRegn med de resterende parametre + super().__init__(**kwargs) + + + @property + def parametre(self): + + return dict() + + + def korriger(self): + """ + Eksempel hvor apply_geodetic_corrections_to_height_diffs antages at virke med inputs: + point_from_lat + point_from_lon + point_to_lat + point_to_lon + epoch + height_diff + ... + etc + + """ + self.observationer + self.gamle_koter + for i, obs in enumerate(self.observationer): + + fra = obs.fra + til = obs.til + epoch = obs.dato + height_diff = obs.deltaH + + (observationer, korrektioner) = apply_geodetic_corrections_to_height_diffs( + fra, + til, + epoch, + height_diff, + self.tidal_system, + Path( + self.grid_inputfolder + ), # Nødvendigt med Path? Hvis der manuelt angives en sti? + self.height_diff_unit, + self.tidal_system, + self.epoch_target, + self.deformationmodel, + self.gravitymodel, + ) + + # KREBSLW: derefter skal RegneMotorens observations korrigeres, så det virker med GamaRegn.udjævn metoden + # Kan måske gøres som nedenfor, men har ikke testet det? + self.observationer[i].deltaH += korrektioner + + return korrektioner + + def korriger(self): + """ + Eksempel hvor apply_geodetic_corrections_to_height_diffs antages at virke med inputs: + + observationer: list[InternNivObservation] + punkter: list[InternKote] + + og antages at spytte en liste af korrektioner ud: list[float] + """ + + if ( + self.tidal_system is not None + or self.epoch_target is not None + or self.height_diff_unit == "gpu" + ): + print("Højdeforskelle påføres geodætiske korrektioner inden udjævning") + + (observationer, korrektioner) = apply_geodetic_corrections_to_height_diffs( + # KREBSLW: Observationer og arbejdssæt erstattes af InternNivObservation og InternKote + # observationer, + # arbejdssæt, + self.observationer, + self.gamle_koter, + Path( + self.grid_inputfolder + ), # Nødvendigt med Path? Hvis der manuelt angives en sti? + self.height_diff_unit, + self.tidal_system, + self.epoch_target, + self.deformationmodel, + self.gravitymodel, + ) + + # KREBSLW: Herefter skal RegneMotorens observations korrigeres, så det virker med GamaRegn.udjævn metoden + # Kan måske gøres som nedenfor, men har ikke testet det? + for i,(o,k) in enumerate(zip(self.observationer, korrektioner)): + self.observationer[i].deltaH = o.deltaH + k + + # KREBSLW: man kunne godt gemme korrektionerne .... + self.korrektioner = korrektioner + + def gem_korrektioner_som_geojson(self): + # ... så man kan skrive dem ud som en geojson + with open("obs_corrected.geojson", "w") as fil: + fil.write( + f""" + [ + Min geojson + Line: [ + x: 55, y :10, korrektion: {self.korrektioner}, bla:bla + ] + ] + """ + ) + + + def konverter(self): + """ + Samme princip som ovenstående eksempler med `korriger` + """ + # Hvis udjævningen skal foretages i gpu konverteres Helmert-højder til geopotentielle højder + # Pt. foretages konverteringen til geopotentielle højder både ifm. kontrolberegning og ifm. + # den endelige beregning - det er min plan i stedet at gemme de konverterede højder i et faneblad + # mhp. genbrug ifm. den endelige beregning + if self.height_diff_unit == "gpu": + print( + "Højder konverteres fra Helmert-højder til geopotentielle højder inden udjævning" + ) + + # KREBSLW: antag at convert_geopotential_heights_to_metric_heights arbejder med list[InterKote] + liste_af_konverterede_højder = convert_geopotential_heights_to_metric_heights( + # arbejdssæt, + self.gamle_koter, + "helmert_to_geopot", + Path( + self.grid_inputfolder + ), # Nødvendigt med Path? Hvis der manuelt angives en sti? + self.gravitymodel, + self.tidal_system, + ) + + # KREBSLW: I samme stil som i `korriger` skal regnemotoren interne repræsentation af koterne opdateres + for i, (punkt, konverteret_højde) in enumerate(zip(self.gamle_koter, liste_af_konverterede_højder)): + self.gamle_koter[i].H = konverteret_højde + + + def konverter_tilbage(self): + # Geopotentielle højder fra udjævningen konverteres til Helmert- eller normalhøjder + if self.height_diff_unit == "gpu": + print( + "Højder konverteres fra geopotentielle højder til Helmert-højder (eller normalhøjder) efter udjævning" + ) + + liste_af_konverterede_højder = convert_geopotential_heights_to_metric_heights( + # beregning, + self.nye_koter, # antages at indeholde de udjævnede koter + f"geopot_to_{self.output_height}", + Path( + self.grid_inputfolder + ), # Nødvendigt med Path? Hvis der manuelt angives en sti? + self.gravitymodel, + self.tidal_system, + ) + + # KREBSLW: I samme stil som i ovenstående skal regnemotoren interne repræsentation af de nye koter opdateres + for i, (punkt, konverteret_højde) in enumerate(zip(self.nye_koter, liste_af_konverterede_højder)): + self.nye_koter[i].H = konverteret_højde + + def udjævn(self): + """ + Den her samler det hele og er dén funktion som bliver kaldt i fire niv regn + """ + self.korriger() + self.konverter() + + # Brug udjævningsfunktionalitet fra GamaRegn + GamaRegn.udjævn(self) + # GamaRegn.udjævn sætter værdien af self.nye_koter, som så bruges i konverter_tilbage + + self.konverter_tilbage() + + + +class DVR90Regn(GeodætiskRegn): + """ Den her klasse gør i princippet ikke andet end at fastsætte nogle parametre """ + def __init__( + self, + epoch_target: pd.Timestamp = '1990', + height_diff_unit: str = 'helmert', + output_height: str = 'helmert', + deformationmodel: str = 'knudsen', + gravitymodel: str = 'forsberg', + **kwargs + ): + + fastlåste_params = dict ( + tidal_system ='non', + grid_inputfolder = 'C:\min\mappe', + ) + + super().__init__( + epoch_target = epoch_target, + height_diff_unit=height_diff_unit, + output_height = output_height, + deformationmodel = deformationmodel, + gravitymodel = gravitymodel, + **fastlåste_params + # ... + **kwargs + ) + +# KREBSLW: Har du tænkt over hvordan spredningerne skal håndteres når der korrigeres og konverteres? def _spredning( observationstype: str, afstand_i_m: float, diff --git "a/src/fire/api/niv/udtr\303\246k_observationer.py" "b/src/fire/api/niv/udtr\303\246k_observationer.py" index 2aebdf6c..1360c9e5 100644 --- "a/src/fire/api/niv/udtr\303\246k_observationer.py" +++ "b/src/fire/api/niv/udtr\303\246k_observationer.py" @@ -113,7 +113,8 @@ def klargør_geometrifiler(geometrifiler: Iterable[str]) -> List[Geometry]: # Gentagne punkter fjernes på forhånd da det er en tydelig fejl delgeometrier = [ shapely.remove_repeated_points( - shapely.geometry.shape(delgeometri.get("geometry")), 0) + shapely.geometry.shape(delgeometri.get("geometry")), 0 + ) for delgeometri in geometri_data ] @@ -130,24 +131,6 @@ def klargør_geometrifiler(geometrifiler: Iterable[str]) -> List[Geometry]: return klargjorte_geometrier -def punkter_til_geojson(data: pd.DataFrame) -> dict: - """Konvertér punkter til geojson-tekststreng.""" - return { - "type": "FeatureCollection", - "features": [ - { - "type": "Feature", - "properties": {k: v for k, v in row.items()}, - "geometry": { - "type": "Point", - "coordinates": row[["Øst", "Nord"]].tolist(), - }, - } - for _, row in data.iterrows() - ], - } - - def søgefunktioner_med_valgte_metoder( forberedt_søgefunktion: Callable, metoder: List[NivMetode] ) -> Iterator[Callable]: @@ -165,7 +148,9 @@ def polygoner(punkter: Iterable[Punkt], buffer: int) -> Iterator[Geometry]: koordinatsæt = (punkt.geometri.koordinater for punkt in punkter) # Opbyg geometri for punkt-koordinater til søgning. - shapely_punkter = (shapely.geometry.Point(*koordinater) for koordinater in koordinatsæt) + shapely_punkter = ( + shapely.geometry.Point(*koordinater) for koordinater in koordinatsæt + ) # Lav den endelige søge-geometri ved at bruge den angivne buffer som # radius i en forsimplet cirkel (polygon) omkring koordinaterne. diff --git a/src/fire/cli/info.py b/src/fire/cli/info.py index 57df4930..e9ddac80 100644 --- a/src/fire/cli/info.py +++ b/src/fire/cli/info.py @@ -14,6 +14,7 @@ import fire.cli from fire.ident import klargør_ident_til_søgning +from fire.io.geojson import skriv_sagsrapport_geojson from fire.api.model import ( EventType, EVENTTYPER, @@ -296,6 +297,8 @@ def punktsamlingsrapport(punktsamlinger: list[PunktSamling], id: str = None): fire.cli.print(header, bold=True) fire.cli.print(subheader) + punktsamlinger = [ps for ps in punktsamlinger if ps.registreringtil is None] + # Sortér Punktsamlinger efter Jessennummer, dernæst efter Punktsamlingsnavn punktsamlinger.sort(key=lambda x: (x.jessenpunkt.jessennummer, x.navn)) @@ -341,6 +344,8 @@ def tidsserietype(tstype): return "Højde" for ts in tidsserier: + if ts.registreringtil is not None: + continue navn_ombrudt = textwrap.wrap(str(ts.navn), kolonnebredder[0]) for navn_del in navn_ombrudt[:-1]: fire.cli.print(navn_del) @@ -911,16 +916,18 @@ def _grupper_sagsevents(sagsevent_liste: list[Sagsevent]) -> dict[set]: @click.option( "-r", "--rapport", - is_flag=True, - default=False, - help="Udskriv rapport over fremsøgte sager.", + type=str, + default=None, + is_flag=False, + flag_value="", + help="Udskriv rapport over fremsøgte sager. Gives et navn på rapporten skrives også en geojson-fil som oversigt.", ) def sag( sagsid: str, fra: datetime.datetime, til: datetime.datetime, aktive: bool, - rapport: bool, + rapport: str, **kwargs, ): """ @@ -951,6 +958,7 @@ def sag( eller indgår i nogen af ovenstående kategorier. **NB!** Optællingen kan ved store sager eller mange fremsøgte sager godt tage lidt tid. + Angiver man et filnavn ``--rapport sagsrapport.geojson``, skrives en rapporten som geojson med det angivne filnavn. **EKSEMPEL** @@ -968,7 +976,7 @@ def sag( Vis alle sager indeholdende søgeteksten "2024_VEDL" og lav samlet optælling:: - fire info sag 2024_VEDL --rapport + fire info sag 2024_VEDL --rapport "2024_VEDLIGEHOLD_rapport.geojson" """ sager = [] @@ -1027,11 +1035,22 @@ def sag( sagseventid = sagsevent.id[0:8] fire.cli.print(f"[{tid}|{sagseventid}] {eventtype}: {beskrivelse}") - if rapport: - fire.cli.print(f"\n Sagsoptælling :\n") - stats = _optæl_punkter_i_sagsevents(sag.sagsevents) - for k, v in stats.items(): - fire.cli.print(f" Antal {k}: {len(v)}") + # hvis --rapport er sat, så tæller vi sagen op + if rapport is None: + return + + fire.cli.print(f"\n Sagsoptælling :\n") + stats = _optæl_punkter_i_sagsevents(sag.sagsevents) + for k, v in stats.items(): + fire.cli.print(f" Antal {k}: {len(v)}") + + # hvis --rapport "filnavn" er sat så skriver vi også en geojson fil ud + if rapport == "": + return + + punktider = list(set().union(*stats.values())) + punkter = fire.cli.firedb.hent_punkter_fra_uuid_liste(punktider) + skriv_sagsrapport_geojson(rapport, punkter, stats) return @@ -1046,7 +1065,8 @@ def sag( beskrivelse = sag.beskrivelse[0:70].strip().replace("\n", " ").replace("\r", "") fire.cli.print(f"{sag.id[0:8]}: {sag.behandler:20}{beskrivelse}...") - if not rapport: + # hvis --rapport er sat, så tæller vi sagerne op + if rapport is None: return fire.cli.print(f"\n--- Optælling af alle fremsøgte sager ---") @@ -1056,6 +1076,14 @@ def sag( for k, v in stats.items(): print(f"Antal {k}: {len(v)}") + # hvis --rapport "filnavn" er sat så skriver vi også en geojson fil ud + if rapport == "": + return + + punktider = list(set().union(*stats.values())) + punkter = fire.cli.firedb.hent_punkter_fra_uuid_liste(punktider) + skriv_sagsrapport_geojson(rapport, punkter, stats) + @info.command() @fire.cli.default_options() diff --git a/src/fire/cli/niv/__init__.py b/src/fire/cli/niv/__init__.py index 0a23b578..a5f17744 100644 --- a/src/fire/cli/niv/__init__.py +++ b/src/fire/cli/niv/__init__.py @@ -1,6 +1,4 @@ import datetime -import json -import math import os import os.path from pathlib import Path @@ -12,7 +10,6 @@ import click import pandas as pd from sqlalchemy.orm.exc import NoResultFound -from sqlalchemy.exc import DatabaseError import packaging.version from fire.api.model import ( @@ -23,7 +20,6 @@ HøjdeTidsserie, ) from fire.io.regneark import arkdef -from fire.ident import kan_være_gi_nummer import fire.cli from fire.cli import firedb, grøn @@ -382,178 +378,6 @@ def find_sagsid(sagsgang: pd.DataFrame) -> str: return "" -# ------------------------------------------------------------------------------ -def _geojson_filnavn(projektnavn: str, infiks: str, variant: str): - """Generer filnavn på geojson-fil""" - return f"{projektnavn}{infiks}-{variant}.geojson" - - -def punkt_feature(punkter: pd.DataFrame) -> Dict[str, str]: - """Omsæt punktinformationer til JSON-egnet dict""" - - def _none_eller_nan(værdi: float) -> bool: - """ - Check om værdi er None eller NaN. - - Vi checker både for None og NaN, da Pandas og Numpy kan være lidt - drilske på dette område, og har udvist skiftende adfærd gennem tiden. - """ - return værdi is None or math.isnan(værdi) - - for i in range(punkter.shape[0]): - # Nye punkter har hverken ny eller gammel kote. - # Vi rammer ind i denne situation ved læsning af observationer til nye punkter, - # der endnu ikke er regnet en kote for. - if _none_eller_nan(punkter.at[i, "Kote"]) and _none_eller_nan( - punkter.at[i, "Ny kote"] - ): - fastholdt = False - delta = None - kote = None - sigma = None - - # Fastholdte punkter har ingen ny kote, så vi viser den gamle - elif _none_eller_nan(punkter.at[i, "Ny kote"]) and not _none_eller_nan( - punkter.at[i, "Kote"] - ): - fastholdt = True - delta = 0.0 - kote = float(punkter.at[i, "Kote"]) - sigma = float(punkter.at[i, "σ"]) - - # Gamle punkter med nye koter er "standardtilfældet" - else: - fastholdt = False - delta = float(punkter.at[i, "Δ-kote [mm]"]) - kote = float(punkter.at[i, "Ny kote"]) - sigma = float(punkter.at[i, "Ny σ"]) - - # Ignorerede ændringer (under 1 um) - if _none_eller_nan(delta): - delta = None - - # Forbered punktnumre til attributtabellen. Hvis muligt finder vi information - # i databasen og bruger punktets landsnummer som ID, ellers bruges strengen - # der kommer fra Dataframe'n. - try: - punkt = firedb.hent_punkt(punkter.at[i, "Punkt"]) - landsnr = punkt.landsnummer - gi_nummer = punkt.ident if kan_være_gi_nummer(punkt.ident) else None - except (NoResultFound, DatabaseError): - landsnr = punkter.at[i, "Punkt"] - gi_nummer = None - - feature = { - "type": "Feature", - "properties": { - "id": landsnr, - "GI": gi_nummer, - "H": kote, - "sH": sigma, - "Δ": delta, - "fastholdt": fastholdt, - }, - "geometry": { - "type": "Point", - "coordinates": [punkter.at[i, "Øst"], punkter.at[i, "Nord"]], - }, - } - yield feature - - -def punkter_geojson( - punkter: pd.DataFrame, -) -> str: - """Returner punkter/koordinater som geojson-streng""" - # with open(f"{projektnavn}{infiks}-punkter.geojson", "wt") as punktfil: - til_json = { - "type": "FeatureCollection", - "Features": list(punkt_feature(punkter)), - } - return json.dumps(til_json, indent=4) - - -def skriv_punkter_geojson(projektnavn: str, punkter: pd.DataFrame, infiks: str = ""): - """Skriv geojson-fil med punktdata til disk""" - geojson = punkter_geojson(punkter) - filnavn = _geojson_filnavn(projektnavn, infiks, "punkter") - with open(filnavn, "wt") as punktfil: - punktfil.write(geojson) - - -# ------------------------------------------------------------------------------ -def obs_feature( - punkter: pd.DataFrame, observationer: pd.DataFrame, antal_målinger: Dict[Tuple, int] -) -> Dict[str, str]: - """Omsæt observationsinformationer til JSON-egnet dict""" - for i in range(observationer.shape[0]): - fra = observationer.at[i, "Fra"] - til = observationer.at[i, "Til"] - feature = { - "type": "Feature", - "properties": { - "Fra": fra, - "Til": til, - "Målinger": antal_målinger[tuple(sorted([fra, til]))], - "Afstand": observationer.at[i, "L"], - "ΔH": observationer.at[i, "ΔH"], - "Observationstidspunkt": str(observationer.at[i, "Hvornår"]), - # konvertering, da json.dump ikke uderstøtter int64 - "Opstillinger": int(observationer.at[i, "Opst"]), - "Journal": observationer.at[i, "Journal"], - "Type": observationer.at[i, "Type"], - "Slukket": observationer.at[i, "Sluk"], - }, - "geometry": { - "type": "LineString", - "coordinates": [ - [punkter.at[fra, "Øst"], punkter.at[fra, "Nord"]], - [punkter.at[til, "Øst"], punkter.at[til, "Nord"]], - ], - }, - } - yield feature - - -def observationer_geojson( - punkter: pd.DataFrame, - observationer: pd.DataFrame, -) -> None: - """Skriv observationer til geojson-fil""" - - fra = observationer["Fra"] - til = observationer["Til"] - - # Optæl antal frem-og/eller-tilbagemålinger pr. strækning: Vi starter - # med en dict med et nul for hver strækning - par = [tuple(p) for p in zip(fra, til)] - antal_målinger = dict((tuple(sorted(p)), 0) for p in par) - # ...og så tæller vi det relevante element op for hver observation - for p in par: - # Indeksering med tuple(sorted(p)) da set(p) ikke kan hashes - antal_målinger[tuple(sorted(p))] += 1 - - til_json = { - "type": "FeatureCollection", - "Features": list(obs_feature(punkter, observationer, antal_målinger)), - } - - return json.dumps(til_json, indent=4) - - -def skriv_observationer_geojson( - projektnavn: str, - punkter: pd.DataFrame, - observationer: pd.DataFrame, - infiks: str = "", -) -> None: - """Skriv geojson-fil med observationsdata til disk""" - filnavn = _geojson_filnavn(projektnavn, infiks, "observationer") - geojson = observationer_geojson(punkter, observationer) - with open(filnavn, "wt") as obsfil: - obsfil.write(geojson) - - def bekræft(spørgsmål: str, gentag=True) -> bool: """ Anmod bruger om at tage stilling til spørgsmålet. diff --git "a/src/fire/cli/niv/_l\303\246s_observationer.py" "b/src/fire/cli/niv/_l\303\246s_observationer.py" index 09113295..3fbd1ddf 100644 --- "a/src/fire/cli/niv/_l\303\246s_observationer.py" +++ "b/src/fire/cli/niv/_l\303\246s_observationer.py" @@ -14,14 +14,17 @@ nyt_ark, arkdef, ) +from fire.io.geojson import ( + skriv_punkter_geojson, + skriv_observationer_geojson, +) + import fire.io.dataframe as frame import fire.cli from fire.cli.niv import ( find_faneblad, niv, - skriv_punkter_geojson, - skriv_observationer_geojson, skriv_ark, er_projekt_okay, KOTESYSTEMER, diff --git a/src/fire/cli/niv/_opret_sag.py b/src/fire/cli/niv/_opret_sag.py index 4507875b..be90e992 100644 --- a/src/fire/cli/niv/_opret_sag.py +++ b/src/fire/cli/niv/_opret_sag.py @@ -162,7 +162,16 @@ def opret_sag(projektnavn: str, beskrivelse: str, sagsbehandler: str, **kwargs) filoversigt = nyt_ark(arkdef.FILOVERSIGT) param = pd.DataFrame( columns=["Navn", "Værdi"], - data=[("Version", fire.__version__), ("Database", fire.cli.firedb.db)], + data=[ + ("Version", fire.__version__), + ("Database", fire.cli.firedb.db), + ("Tidesystem", None), + ("Epoke", None), + ("Enhed højdedifferencer", None), + ("Output højde", None), + ("Deformationsmodel", None), + ("Tyngdemodel", None), + ], ) resultater = { diff --git a/src/fire/cli/niv/_regn.py b/src/fire/cli/niv/_regn.py index c1efb5d5..5c49eab7 100644 --- a/src/fire/cli/niv/_regn.py +++ b/src/fire/cli/niv/_regn.py @@ -2,13 +2,15 @@ import webbrowser import click -from pandas import DataFrame +from pandas import DataFrame, Timestamp, to_datetime from fire.api.model import ( HøjdeTidsserie, Koordinat, ) from fire.api.niv.regnemotor import ( + GeodætiskRegn, + DVR90Regn, RegneMotor, GamaRegn, DumRegn, @@ -17,17 +19,20 @@ ) from fire.io.regneark import arkdef +from fire.io.geojson import ( + skriv_punkter_geojson, + skriv_observationer_geojson, +) import fire.cli from fire.cli.ts.plot_ts import ( plot_tidsserier, ) + from fire.cli.niv import ( find_faneblad, niv, - skriv_punkter_geojson, - skriv_observationer_geojson, skriv_ark, er_projekt_okay, hent_relevante_tidsserier, @@ -36,9 +41,13 @@ from fire.cli.niv._netoversigt import byg_netgeometri_og_singulære +DEFAULT_STI_GRIDS = Path(__file__).parents[2] / Path("data/") + motorvælger = { "gama": GamaRegn, "dum": DumRegn, + "geod": GeodætiskRegn, + "dvr90": DVR90Regn, } @@ -62,9 +71,30 @@ default=False, help="Angiv om beregnede koter skal plottes som forlængelse af en tidsserie", ) -def regn(projektnavn: str, plot: bool, MotorKlasse: type[RegneMotor], **kwargs) -> None: +# KREBSLW: option taget fra Evers' PR. har slettet de andre +@click.option( + "-r", + "--regneparameter", + "regneparametre", + multiple=True, + help="Regnemotorspecifikke parametre. Sættes på formen 'parameter=værdi'. Flere parametre kan sættes i samme kommando." +) +def regn( + projektnavn: str, + MotorKlasse: type[RegneMotor], + plot: bool, + regneparametre: tuple[str], + **kwargs, +) -> None: """Beregn nye koter. + Forudsætninger vedr. geodætisk korrektion/konvertering: + Input højder ("Højde" i fanebladet "Punktoversigt") er altid Helmert-højder, + output højder ("Ny højde" i fanebladene "Kontrolberegning" samt "Endelig beregning") skal enten + være Helmert-højder eller normalhøjder. Især sidstnævnte forudsætning skyldes, at det umiddelbart + ville blive for "kringlet" at håndetere, hvis output højder også skulle kunne være geopotentielle + højder. + Forudsat nivellementsobservationer allerede er indlæst i sagsregnearket kan der beregnes nye koter på baggrund af disse observationer. Beregning af koter med dette program er en totrinsprocedure. Først udføres en @@ -175,10 +205,87 @@ def regn(projektnavn: str, plot: bool, MotorKlasse: type[RegneMotor], **kwargs) infiks = "" beregningstype = "endelig" - # Håndter fastholdte punkter og slukkede observationer. + # Indlæsning af observationer, punktoversigt samt parametre observationer = find_faneblad(projektnavn, "Observationer", arkdef.OBSERVATIONER) punktoversigt = find_faneblad(projektnavn, "Punktoversigt", arkdef.PUNKTOVERSIGT) arbejdssæt = find_faneblad(projektnavn, aktuelt_faneblad, arkdef.PUNKTOVERSIGT) + parametre = find_faneblad(projektnavn, "Parametre", arkdef.PARAM) + + + # KREBSLW: Taget fra Ever's PR. + # Tjek om der er parametre til regnemotoren. Parametre der ikke kan læses + # springes over. + motorkwargs = {} + for regneparameter in regneparametre: + try: + parameter, værdi = regneparameter.split("=") + except ValueError: + fire.cli.print( + ( + f"ADVARSEL: regneparameteren '{regneparameter} kan ikke tolkes. " + "Skal være på formen 'parameter=værdi'." + ), + bold=True, + bg="yellow", + ) + + # konverter til float hvis der er givet en talværdi + try: + værdi = float(værdi) + except ValueError: + pass + + motorkwargs[parameter] = værdi + + # KREBSLW : Hmm tror det her kommer til at blive besværligt også for brugeren. + # Hvis man gerne vil lave mange kontrol-beregninger i træk, så skal man huske at slette/omdøbe Fanebladet "Kontrolberegning", da den ellers vil brokke sig + + # Hvis der skal foretages en kontrolberegning gemmes parametre vedr. geodætisk korrektion/konvertering + # if kontrol: + # [ + # parametre.loc[parametre["Navn"] == "Tidesystem", "Værdi"], + # parametre.loc[parametre["Navn"] == "Epoke", "Værdi"], + # parametre.loc[parametre["Navn"] == "Enhed højdedifferencer", "Værdi"], + # parametre.loc[parametre["Navn"] == "Output højde", "Værdi"], + # parametre.loc[parametre["Navn"] == "Deformationsmodel", "Værdi"], + # parametre.loc[parametre["Navn"] == "Tyngdemodel", "Værdi"], + # ] = [ + # tidal_system, + # epoch_target, + # height_diff_unit, + # output_height, + # deformationmodel, + # gravitymodel, + # ] + + # # I modsat fald kontrolleres at parametrene er konsistente med kontrolberegningen + # # Pt. fejler kontrollen fejlagtigt, hvis ingen parametre angives i kontrolberegningen/den endelige beregning + # else: + # # fmt: off + # if ( + # [ + # parametre.loc[parametre["Navn"] == "Tidesystem", "Værdi"].item(), + # to_datetime(parametre.loc[parametre["Navn"] == "Epoke", "Værdi"]).item(), + # parametre.loc[parametre["Navn"] == "Enhed højdedifferencer", "Værdi"].item(), + # parametre.loc[parametre["Navn"] == "Output højde", "Værdi"].item(), + # parametre.loc[parametre["Navn"] == "Deformationsmodel", "Værdi"].item(), + # parametre.loc[parametre["Navn"] == "Tyngdemodel", "Værdi"].item(), + # ] + # ) != ( + # [ + # tidal_system, + # epoch_target, + # height_diff_unit, + # output_height, + # deformationmodel, + # gravitymodel, + # ] + # ): + # # fmt: on + # fire.cli.print( + # "Fejl: Beregningsparametrene for den endelige beregning skal være konsistente med kontrolberegningen" + # ) + # raise SystemExit(1) # Til den endelige beregning skal vi bruge de oprindelige observationsdatoer if not kontrol: @@ -188,9 +295,24 @@ def regn(projektnavn: str, plot: bool, MotorKlasse: type[RegneMotor], **kwargs) observationer_uden_slukkede = observationer[observationer["Sluk"] != "x"] # Start regnemotoren! - motor = MotorKlasse.fra_dataframe( - observationer_uden_slukkede, arbejdssæt, projektnavn=projektnavn - ) + # KREBSLW: kopieret fra Evers' eksempel! + try: + motor = MotorKlasse.fra_dataframe( + observationer_uden_slukkede, arbejdssæt, projektnavn=projektnavn, **motorkwargs + ) + except TypeError as error: + # Fejlbeskeden vi kan få er på formen: + # + # TypeError: RegneMotor.__init__() got an unexpected keyword argument 'param' + # + # ... og derfor kan vi slippe afsted med at splitte stringen på ' + parameter_navn = str(error).split("'")[1] + fire.cli.print( + f"FEJL: regneparameteren '{parameter_navn}' er ukendt.", + bold=True, + bg="red", + ) + raise SystemExit # Tilføj "-kontrol" eller "-endelig" til alle filnavne motor.filer = [ @@ -237,8 +359,34 @@ def regn(projektnavn: str, plot: bool, MotorKlasse: type[RegneMotor], **kwargs) # Opdater arbejdssæt med udjævningsresultat beregning = opdater_arbejdssæt(arbejdssæt, nye_punkter_df) beregning = beregning.reset_index() + + # Endvidere genindføres oprindelige Helmert-højder + beregning["Kote"] = punktoversigt["Kote"] + + # Beregningsresultater gemmes resultater[næste_faneblad] = beregning + # Hvis der er tale om en kontrolberegning gemmes endvidere parametre vedr. geodætisk + # korrektion/konvertering + if kontrol: + resultater["Parametre"] = parametre + # Hvis endvidere observationerne er påført geodætiske korrektioner gemmes de korrigerede + # observationer samt korrektionerne + + # KREBSLW : Måske motoren skal gemme det i nogle geojson filer (se ovre i GeodætiskRegn)? + # if ( + # tidal_system is not None + # or epoch_target is not None + # or height_diff_unit == "gpu" + # ): + # korrigerede_observationer = observationer[["ΔH"]].copy() + # korrigerede_observationer = korrigerede_observationer.join(korrektioner) + # resultater["Korrigerede observationer"] = korrigerede_observationer + + # korrigerede_observationer = korrektioner.insert( + # 3, "ΔH", korrigerede_observationer + # ) + # Plot tidsserier forlænget med de nyberegnede koter. if plot == True: kotesystem = fire.cli.firedb.hent_srid(beregning["System"][0]) @@ -314,6 +462,11 @@ def regn(projektnavn: str, plot: bool, MotorKlasse: type[RegneMotor], **kwargs) observationer, infiks=infiks, ) + + # KREBSLW: Man kan gøre så Parametre altid gemmes, hvis motoren har nogle som den vil fortælle om + # Noget i stil med nedenstående + resultater["Parametre"] = DataFrame.from_dict(motor.parametre) + skriv_ark(projektnavn, resultater) if fire.cli.firedb.config.getboolean("general", "niv_open_files"): # åbn html output hvis motoren producerer et diff --git "a/src/fire/cli/niv/_udtr\303\246k_observationer.py" "b/src/fire/cli/niv/_udtr\303\246k_observationer.py" index 53b369b8..d7202028 100644 --- "a/src/fire/cli/niv/_udtr\303\246k_observationer.py" +++ "b/src/fire/cli/niv/_udtr\303\246k_observationer.py" @@ -47,10 +47,13 @@ from fire.cli.niv import ( niv as niv_command_group, er_projekt_okay, - skriv_observationer_geojson, - skriv_punkter_geojson, KOTESYSTEMER, ) +from fire.io.geojson import ( + skriv_punkter_geojson, + skriv_observationer_geojson, +) + from fire.typologi import ( adskil_filnavne, adskil_identer, @@ -272,12 +275,10 @@ def udtræk_observationer( fire.cli.print("Søg med punkterne som opstillingspunkt", fg="yellow") objekter = punkter # Forbereder søgefunktion med argumenter fastsat. - DVR90 = db.hent_srid(SRID.DVR90) funktion = db.hent_observationer_fra_opstillingspunkt fastholdte_argumenter = dict( tid_fra=fra, tid_til=til, - srid=DVR90, kun_aktive=True, sigtepunkter=punkter, ) diff --git a/src/fire/data/DKup24geo_DTU2024_PK.tif b/src/fire/data/DKup24geo_DTU2024_PK.tif new file mode 100644 index 00000000..4b499ffe Binary files /dev/null and b/src/fire/data/DKup24geo_DTU2024_PK.tif differ diff --git a/src/fire/data/NKG2016_lev.tif b/src/fire/data/NKG2016_lev.tif new file mode 100644 index 00000000..3f8cc33b Binary files /dev/null and b/src/fire/data/NKG2016_lev.tif differ diff --git a/src/fire/data/README.md b/src/fire/data/README.md new file mode 100644 index 00000000..6754e26f --- /dev/null +++ b/src/fire/data/README.md @@ -0,0 +1,46 @@ +This folder contains the following grids: + + +Grid/file name: absg_DTU2016_PK.tif +Model name: +Version: absg +Released: 2016 +Author/institution: Per Knudsen, DTU Space +Description: land uplift relative to geoid in mm/year +Grid extent/spacing: 54.099998 57.999999 7.700000 13.100000 0.100000001 0.200000003 +Epoch: 1900-2000 (in principle) +Reference: F:\GRF\Data\GEO\Uplift-modeller\DK-model\2016_model_DTU\Upliftmodel2016.pdf + + +Grid/file name: DKup24geo_DTU2024_PK.tif +Model name: DKup24 uplift model for Denmark +Version: DKup24geo +Released: 2024 +Author/institution: Per Knudsen, Karsten E. Engsager and S. Abbas Khan, DTU Space +Description: land uplift relative to geoid in mm/year +Grid extent/spacing: 54.10000 58.000000 7.700000 13.100000 0.100000 0.200000 +Epoch: 1900-2000 (in principle) +Reference: https://doi.org/10.11583/DTU.27277848.v1 + + +Grid/file name: NKG2016_lev.tif +Model name: NKG2016LU +Version: NKG2016LU_lev +Released: 2016 +Author/institution: Olav Vestøl et al., NKG +Description: land uplift relative to geoid in mm/year +Grid extent/spacing: 49.00000000 75.00000000 0.00000000 50.00000000 0.083333333300 0.166666666700 +Epoch: +Reference: https://link.springer.com/article/10.1007/s00190-019-01280-8 + + +Grid/file name: dk-g-direkte-fra-gri-thokn.tif +Model name: +Version: +Released: 2023 +Author/institution: René Forsberg, DTU Space +Description: gravity at surface level in units of mGal +Grid extent/spacing: 54.500000 57.800000 7.800000 15.500000 0.005000 0.010000 +Epoch: +Tidal system: zero tide +Reference: F:\GRF\Data\GEO\TYNGDE\FRA_DTUSPACE\GRID_2023\tyngdekraft-grid.pdf diff --git a/src/fire/data/absg_DTU2016_PK.tif b/src/fire/data/absg_DTU2016_PK.tif new file mode 100644 index 00000000..8b7a416a Binary files /dev/null and b/src/fire/data/absg_DTU2016_PK.tif differ diff --git a/src/fire/data/dk-g-direkte-fra-gri-thokn.tif b/src/fire/data/dk-g-direkte-fra-gri-thokn.tif new file mode 100644 index 00000000..94ad281a Binary files /dev/null and b/src/fire/data/dk-g-direkte-fra-gri-thokn.tif differ diff --git a/src/fire/io/geojson.py b/src/fire/io/geojson.py new file mode 100644 index 00000000..60283c30 --- /dev/null +++ b/src/fire/io/geojson.py @@ -0,0 +1,234 @@ +""" +Modul til håndtering af læsning og skrivning geojson-filer +""" + +import json +import math + +import pandas as pd +from sqlalchemy.orm.exc import NoResultFound +from sqlalchemy.exc import DatabaseError + +from fire.api.model import ( + Punkt, +) +from fire.ident import kan_være_gi_nummer +from fire.cli import firedb + + +def _geojson_filnavn(projektnavn: str, infiks: str, variant: str): + """Generer filnavn på geojson-fil""" + return f"{projektnavn}{infiks}-{variant}.geojson" + + +def punkt_feature(punkter: pd.DataFrame) -> dict[str, str]: + """Omsæt punktinformationer til JSON-egnet dict""" + + def _none_eller_nan(værdi: float) -> bool: + """ + Check om værdi er None eller NaN. + + Vi checker både for None og NaN, da Pandas og Numpy kan være lidt + drilske på dette område, og har udvist skiftende adfærd gennem tiden. + """ + return værdi is None or math.isnan(værdi) + + for i in range(punkter.shape[0]): + # Nye punkter har hverken ny eller gammel kote. + # Vi rammer ind i denne situation ved læsning af observationer til nye punkter, + # der endnu ikke er regnet en kote for. + if _none_eller_nan(punkter.at[i, "Kote"]) and _none_eller_nan( + punkter.at[i, "Ny kote"] + ): + fastholdt = False + delta = None + kote = None + sigma = None + + # Fastholdte punkter har ingen ny kote, så vi viser den gamle + elif _none_eller_nan(punkter.at[i, "Ny kote"]) and not _none_eller_nan( + punkter.at[i, "Kote"] + ): + fastholdt = True + delta = 0.0 + kote = float(punkter.at[i, "Kote"]) + sigma = float(punkter.at[i, "σ"]) + + # Gamle punkter med nye koter er "standardtilfældet" + else: + fastholdt = False + delta = float(punkter.at[i, "Δ-kote [mm]"]) + kote = float(punkter.at[i, "Ny kote"]) + sigma = float(punkter.at[i, "Ny σ"]) + + # Ignorerede ændringer (under 1 um) + if _none_eller_nan(delta): + delta = None + + # Forbered punktnumre til attributtabellen. Hvis muligt finder vi information + # i databasen og bruger punktets landsnummer som ID, ellers bruges strengen + # der kommer fra Dataframe'n. + try: + punkt = firedb.hent_punkt(punkter.at[i, "Punkt"]) + landsnr = punkt.landsnummer + gi_nummer = punkt.ident if kan_være_gi_nummer(punkt.ident) else None + except (NoResultFound, DatabaseError): + landsnr = punkter.at[i, "Punkt"] + gi_nummer = None + + feature = { + "type": "Feature", + "properties": { + "id": landsnr, + "GI": gi_nummer, + "H": kote, + "sH": sigma, + "Δ": delta, + "fastholdt": fastholdt, + }, + "geometry": { + "type": "Point", + "coordinates": [punkter.at[i, "Øst"], punkter.at[i, "Nord"]], + }, + } + yield feature + + +def punkter_geojson( + punkter: pd.DataFrame, +) -> str: + """Returner punkter/koordinater som geojson-streng""" + til_json = { + "type": "FeatureCollection", + "Features": list(punkt_feature(punkter)), + } + return json.dumps(til_json, indent=4) + + +def skriv_punkter_geojson(projektnavn: str, punkter: pd.DataFrame, infiks: str = ""): + """Skriv geojson-fil med punktdata til disk""" + geojson = punkter_geojson(punkter) + filnavn = _geojson_filnavn(projektnavn, infiks, "punkter") + with open(filnavn, "wt") as punktfil: + punktfil.write(geojson) + + +def obs_feature( + punkter: pd.DataFrame, observationer: pd.DataFrame, antal_målinger: dict[tuple, int] +) -> dict[str, str]: + """Omsæt observationsinformationer til JSON-egnet dict""" + for i in range(observationer.shape[0]): + fra = observationer.at[i, "Fra"] + til = observationer.at[i, "Til"] + feature = { + "type": "Feature", + "properties": { + "Fra": fra, + "Til": til, + "Målinger": antal_målinger[tuple(sorted([fra, til]))], + "Afstand": observationer.at[i, "L"], + "ΔH": observationer.at[i, "ΔH"], + "Observationstidspunkt": str(observationer.at[i, "Hvornår"]), + # konvertering, da json.dump ikke uderstøtter int64 + "Opstillinger": int(observationer.at[i, "Opst"]), + "Journal": observationer.at[i, "Journal"], + "Type": observationer.at[i, "Type"], + "Slukket": observationer.at[i, "Sluk"], + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [punkter.at[fra, "Øst"], punkter.at[fra, "Nord"]], + [punkter.at[til, "Øst"], punkter.at[til, "Nord"]], + ], + }, + } + yield feature + + +def observationer_geojson( + punkter: pd.DataFrame, + observationer: pd.DataFrame, +) -> None: + """Skriv observationer til geojson-fil""" + + fra = observationer["Fra"] + til = observationer["Til"] + + # Optæl antal frem-og/eller-tilbagemålinger pr. strækning: Vi starter + # med en dict med et nul for hver strækning + par = [tuple(p) for p in zip(fra, til)] + antal_målinger = dict((tuple(sorted(p)), 0) for p in par) + # ...og så tæller vi det relevante element op for hver observation + for p in par: + # Indeksering med tuple(sorted(p)) da set(p) ikke kan hashes + antal_målinger[tuple(sorted(p))] += 1 + + til_json = { + "type": "FeatureCollection", + "Features": list(obs_feature(punkter, observationer, antal_målinger)), + } + + return json.dumps(til_json, indent=4) + + +def skriv_observationer_geojson( + projektnavn: str, + punkter: pd.DataFrame, + observationer: pd.DataFrame, + infiks: str = "", +) -> None: + """Skriv geojson-fil med observationsdata til disk""" + filnavn = _geojson_filnavn(projektnavn, infiks, "observationer") + geojson = observationer_geojson(punkter, observationer) + with open(filnavn, "wt") as obsfil: + obsfil.write(geojson) + + +def skriv_sagsrapport_geojson(filnavn: str, punkter: list[Punkt], attributter: dict): + """ + Skriv geojson-fil med sagsstatistik for punkter til disk + + Er en hjælpefunktion til ``fire info sag`` + + Punkternes attributter forventes indeholdt i en dict med attributternes navne som + nøgle og lister af punkt-uuider. Eksempelvis: + + attributter = { + tabtgået = [UUID_A, UUID_B, UUID_C], + oprettet = [UUID_D, UUID_E, UUID_F], + observeret = [UUID_A, UUID_B], + min_attribut = [UUID_X, UUID_Y], + } + + """ + + def _punkt_feature(punkter: list[Punkt], attributter: dict): + for pkt in punkter: + properties = dict(id=str(pkt.landsnummer)) + properties |= { + attribut: (pkt.id in punkter_med_attribut) + for attribut, punkter_med_attribut in attributter.items() + } + + feature = { + "type": "Feature", + "properties": properties, + "geometry": { + "type": "Point", + "coordinates": [ + pkt.geometri.koordinater[0], + pkt.geometri.koordinater[1], + ], + }, + } + yield feature + + til_json = { + "type": "FeatureCollection", + "Features": list(_punkt_feature(punkter, attributter)), + } + geojson = json.dumps(til_json, indent=4) + + with open(filnavn, "wt") as punktfil: + punktfil.write(geojson) diff --git "a/test/api/niv/test_udtr\303\246k_observationer.py" "b/test/api/niv/test_udtr\303\246k_observationer.py" index 3bd9302d..260e558d 100644 --- "a/test/api/niv/test_udtr\303\246k_observationer.py" +++ "b/test/api/niv/test_udtr\303\246k_observationer.py" @@ -31,7 +31,6 @@ observationer_inden_for_spredning, filtrer_præcisionsnivellement, timestamp, - punkter_til_geojson, ) @@ -223,37 +222,3 @@ def test_timestamp_string(): result = timestamp() expected = "2021-11-01T212100" assert result == expected, f"Forventede, at {result!r} var {expected!r}." - - -def test_punkter_til_geojson(ark_punktoversigt, række_punktoversigt): - # Arrange - øst = 8.0 - nord = 55.0 - - række = { - **række_punktoversigt, - "Øst": øst, - "Nord": nord, - } - # `data`-parameteren i pd.DataFrame - # skal være en iterable med rækker. - data = [række.values()] - columns = række.keys() - df_række = pd.DataFrame(data=data, columns=columns) - ark = pd.concat([ark_punktoversigt, df_række], ignore_index=True) - - # Act - geojson = punkter_til_geojson(ark) - print(geojson) - - # Assert - properties = geojson["features"][0]["properties"] - geometry = geojson["features"][0]["geometry"] - check = ( - (properties["Øst"], øst), - (properties["Nord"], nord), - (geometry["coordinates"][0], øst), - (geometry["coordinates"][1], nord), - ) - for result, expected in check: - assert result == expected, f"Forventede, at {result!r} var {expected!r}." diff --git a/testprojekt/geod-lev-test.xlsx b/testprojekt/geod-lev-test.xlsx new file mode 100644 index 00000000..d2d03c1b Binary files /dev/null and b/testprojekt/geod-lev-test.xlsx differ