Skip to content
Snippets Groups Projects
Commit 4c28b03d authored by Rudolf, Michael's avatar Rudolf, Michael
Browse files

Adapted Script for EGMS Data

parent d0ed6e4a
Branches
No related tags found
No related merge requests found
...@@ -638,39 +638,46 @@ def get_envlen(eee: str) -> int: ...@@ -638,39 +638,46 @@ def get_envlen(eee: str) -> int:
def gen_queries_psi_gpkg( def gen_queries_psi_gpkg(
file_path: os.PathLike, direction: str = "vertikal" file_path: os.PathLike, table: str = "vertikal"
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list]: ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list]:
""" """
Opens a gpkg and generates queries to extract the data limited to psi files Opens a gpkg and generates queries to extract the data limited to psi files
:param file_path: The database as a gpkg file. :param file_path: The database as a gpkg file.
:type file_path: os.PathLike :type file_path: os.PathLike
:param direction: The direction to use, defaults to "vertikal" :param table: The direction to use, defaults to "vertikal"
:type direction: str, optional :type table: str, optional
:return: Some data and preformatted queries to extract the data from the database. :return: Some data and preformatted queries to extract the data from the database.
:rtype: Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list] :rtype: Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list]
""" """
logging.info(f"Generating queries to extract {direction} from {file_path}") logging.info(f"Generating queries to extract {table} from {file_path}")
# Get number of rows and names of columns # Get number of rows and names of columns
con = sqlite3.connect(file_path) con = sqlite3.connect(file_path)
cur = con.cursor() cur = con.cursor()
info = read_info(cur, direction) tables = get_table_names(file_path)
if table not in tables:
raise KeyError(
f"Given table does not exist in file! Specify one of {tables}"
)
info = read_info(cur, table)
# Get Coordinates # Get Coordinates
logging.info("Loading coordinates")
try: try:
xx = np.array(select(cur, "X", direction)) xx = np.array(select(cur, "X", table))
except sqlite3.OperationalError as _: except sqlite3.OperationalError as _:
xx = np.array(select(cur, "easting", direction)) xx = np.array(select(cur, "easting", table))
try: try:
yy = np.array(select(cur, "Y", direction)) yy = np.array(select(cur, "Y", table))
except sqlite3.OperationalError as _: except sqlite3.OperationalError as _:
yy = np.array(select(cur, "northing", direction)) yy = np.array(select(cur, "northing", table))
try: try:
zz = np.array(select(cur, "Z", direction)) zz = np.array(select(cur, "Z", table))
except sqlite3.OperationalError as _: except sqlite3.OperationalError as _:
zz = np.array(select(cur, "height", direction)) zz = np.array(select(cur, "height", table))
logging.info("Generating timeseries queries")
if info["has_time"]: if info["has_time"]:
time, queries = gen_timeseries_queries(file_path, direction, info) time, queries = gen_timeseries_queries(file_path, table, info)
con.close() con.close()
......
...@@ -4,10 +4,12 @@ Plots all PSI from the gpkg file into geotiffs and saves them in a folder. ...@@ -4,10 +4,12 @@ Plots all PSI from the gpkg file into geotiffs and saves them in a folder.
import os import os
from datetime import datetime from datetime import datetime
from pathlib import Path
import numpy as np import numpy as np
from tqdm import tqdm from tqdm import tqdm
import u4py.analysis.processing as u4proc
import u4py.io.sql as u4sql import u4py.io.sql as u4sql
import u4py.io.tiff as u4tiff import u4py.io.tiff as u4tiff
import u4py.utils.config as u4config import u4py.utils.config as u4config
...@@ -18,39 +20,77 @@ u4config.start_logger() ...@@ -18,39 +20,77 @@ u4config.start_logger()
def main(): def main():
project = u4proj.get_project( project = u4proj.get_project(
required=["psi_path", "output_path"], required=["base_path"],
# interactive=False, proj_path=Path(
"~/Documents/umwelt4/Convert_EGMS_to_Tiff.u4project"
).expanduser(),
) )
direction = "vertikal" folder_list = [
"Data_EGMS_2015_2021",
"Data_EGMS_2018_2022",
"Data_EGMS_2019_2023",
]
file_list = []
for fp in folder_list:
folder_path = os.path.join(project["paths"]["base_path"], fp)
file_list.extend(
[
os.path.join(folder_path, gpkg_fp)
for gpkg_fp in os.listdir(folder_path)
if gpkg_fp.endswith(".gpkg")
]
)
for fp in tqdm(file_list, desc="Converting gpkg files", leave=False):
egms_gpkg_to_tiff(fp, project["paths"]["base_path"])
def egms_gpkg_to_tiff(gpkg_file_path: os.PathLike, output_path: os.PathLike):
"""
Converts the timeseries data stored in the egms gpkg file to a set of
geotiffs for faster display. Extracts data from all directions in the
given file.
:param gpkg_file_path: The path to the gpkg file.
:type gpkg_file_path: os.PathLike
:param output_path: The path where the folder with the tiffs shall be stored.
:type output_path: os.PathLike
"""
directions = u4sql.get_table_names(gpkg_file_path)
for direction in tqdm(
directions, desc="Extracting directions", leave=False
):
# Setting up paths # Setting up paths
output_dir = os.path.join( output_dir = os.path.join(output_path, f"PSI_Tiffs_{direction}")
project["paths"]["output_path"], f"PSI_Tiffs_{direction}"
)
os.makedirs(output_dir, exist_ok=True) os.makedirs(output_dir, exist_ok=True)
# Reading common data, elevation data (zs) not used atm # Reading common data, elevation data (zs) not used atm
xs, ys, _, time, queries, info = u4sql.gen_queries_psi_gpkg( xs, ys, _, time, queries, info = u4sql.gen_queries_psi_gpkg(
os.path.join(project["paths"]["psi_path"], "hessen_l3_clipped.gpkg"), gpkg_file_path,
direction, direction,
) )
xq = np.array((xs - np.min(xs)) / 50, dtype=int) spacing = np.min(
yq = np.array((ys - np.min(ys)) / 50, dtype=int) [
np.min(np.diff(np.sort(np.unique(xs)))),
np.min(np.diff(np.sort(np.unique(ys)))),
]
)
xq = np.array((xs - np.min(xs)) / spacing, dtype=int)
yq = np.array((ys - np.min(ys)) / spacing, dtype=int)
extent = ( # Adjustments necessary to match up coordinates extent = ( # Adjustments necessary to match up coordinates
np.min(xs) + 25, np.min(xs) + spacing / 2,
np.max(xs) + 75, np.max(xs) + spacing * 1.5,
np.min(ys) - 25, np.min(ys) - spacing / 2,
np.max(ys) + 25, np.max(ys) + spacing / 2,
) )
# Build arguments list (for easier parallelization, WIP) # Build arguments list (for easier parallelization, WIP)
arguments = [(xq, yq, time, extent, output_dir, info, q) for q in queries] arguments = [
for arg in tqdm( (xq, yq, time, extent, output_dir, info, q) for q in queries
arguments, ]
desc="Generating GeoTiffs", u4proc.batch_mapping(
leave=False, arguments, multi_geotiff, desc="Generating GeoTiffs"
): )
multi_geotiff(arg)
def multi_geotiff(arguments: tuple): def multi_geotiff(arguments: tuple):
...@@ -80,8 +120,8 @@ def multi_geotiff(arguments: tuple): ...@@ -80,8 +120,8 @@ def multi_geotiff(arguments: tuple):
vels, vels,
extent, extent,
os.path.join(output_dir, fname), os.path.join(output_dir, fname),
crs=info["proj_crs"], crs="EPSG:3035",
time=time_str, compress="lzw",
) )
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment