Kish Basin halite data

https://hyss.ie/

import glob
import os
from textwrap import wrap
from zipfile import BadZipFile, ZipFile

import cartopy.crs as ccrs
import contextily as cx
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
from geocube.api.core import make_geocube
from shapely.geometry import Polygon

from h2ss import data as rd

plt.rcParams["xtick.major.size"] = 0
plt.rcParams["ytick.major.size"] = 0

# base data download directory
DATA_DIR = os.path.join("data", "kish-basin")

FILE_NAME = "Kish-Basin-dat-files.zip"

URL = f"https://hyss.ie/wp-content/uploads/2023/07/{FILE_NAME}"

DATA_FILE = os.path.join(DATA_DIR, FILE_NAME)

# basemap cache directory
cx.set_cache_dir(os.path.join("data", "basemaps"))

rd.download_data(url=URL, data_dir=DATA_DIR, file_name=FILE_NAME)

Data 'Kish-Basin-dat-files.zip' already exists in 'data/kish-basin'.
Data downloaded on: 2023-09-15 16:03:58.362445+00:00
Download URL: https://hyss.ie/wp-content/uploads/2023/07/Kish-Basin-dat-files.zip
SHA256 hash: 73b3c35059540638bb4f2f85a8417a7bdff051c7f9564880fa71e8c54fcb9f35

ZipFile(DATA_FILE).namelist()

['Flyde Halite Base Depth XYZ Meters.dat',
 'Flyde Halite Thickness - Zone Of Interest - XYZ Meters.dat',
 'Flyde Halite Thickness XYZ Meters.dat',
 'Flyde Halite Top Depth XYZ Meters.dat',
 'Flyde Halite Top TWT XYZ Milliseconds.dat',
 'Kish GIS Map Extent - Square.csv',
 'Mythop Halite Base Depth XYZ Meters.dat',
 'Mythop Halite Thickness XYZ Meters.dat',
 'Mythop Halite Top Depth XYZ Meters.dat',
 'Mythop Halite Top TWT XYZ Milliseconds.dat',
 'Presall Halite Base Depth XYZ Meters-corrected.dat',
 'Presall Halite Thickness - Zone Of Interest - XYZ Meters-corrected.dat',
 'Presall Halite Thickness XYZ Meters-corrected.dat',
 'Presall Halite Top Depth XYZ Meters-corrected.dat',
 'Presall Halite Top TWT XYZ Milliseconds.dat',
 'Rossall Halite Base Depth XYZ Meters.dat',
 'Rossall Halite Thickness - Zone Of Interest - XYZ Meters.dat',
 'Rossall Halite Thickness XYZ Meters.dat',
 'Rossall Halite Top Depth XYZ Meters.dat',
 'Rossall Halite Top TWT XYZ Milliseconds.dat']

# extract archive
try:
    z = ZipFile(DATA_FILE)
    z.extractall(DATA_DIR)
except BadZipFile:
    print("There were issues with the file", DATA_FILE)

Map extent

with open(
    os.path.join(DATA_DIR, "Kish GIS Map Extent - Square.csv"),
    encoding="utf-8",
) as f:
    print(f.read())

Polygon exported from Kingdom on Thu Nov  3 20:00:33 2022
Coordinate System: ED50 / UTM zone 29N   EPSG: 23029
Polygon Name, Index, X, Y
AREA OF INTEREST,1,742000.000000,5937000.000000
AREA OF INTEREST,2,687000.000000,5937000.000000
AREA OF INTEREST,3,687000.000000,5888000.000000
AREA OF INTEREST,4,742000.000000,5888000.000000
AREA OF INTEREST,5,742000.000000,5937000.000000

CRS = 23029

# create extent polygon
extent = pd.read_csv(
    os.path.join(DATA_DIR, "Kish GIS Map Extent - Square.csv"), skiprows=2
)
extent = gpd.GeoSeries(
    Polygon(
        [
            (extent[" X"][0], extent[" Y"][0]),
            (extent[" X"][1], extent[" Y"][1]),
            (extent[" X"][2], extent[" Y"][2]),
            (extent[" X"][3], extent[" Y"][3]),
        ]
    ),
    crs=CRS,
)

extent.bounds

	minx	miny	maxx	maxy
0	687000.0	5888000.0	742000.0	5937000.0

extent.crs

<Projected CRS: EPSG:23029>
Name: ED50 / UTM zone 29N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Europe - between 12°W and 6°W - Faroe Islands - onshore; Spain - mainland onshore; Ireland offshore.
- bounds: (-12.0, 36.13, -6.0, 62.41)
Coordinate Operation:
- name: UTM zone 29N
- method: Transverse Mercator
Datum: European Datum 1950
- Ellipsoid: International 1924
- Prime Meridian: Greenwich

plt.figure(figsize=(7, 7))
ax = plt.axes(projection=ccrs.epsg(3857))
plt.xlim(-1.2e6, -0.55e6)
plt.ylim(6.65e6, 7.475e6)
extent.to_crs(3857).boundary.plot(ax=ax, color="darkslategrey")
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, zoom=7)
plt.title("Kish GIS Map Extent")
plt.tight_layout()
plt.show()

XYZ data

def read_dat_file(dat_path: str, dat_crs):
    """
    Read XYZ data layers into an Xarray dataset
    """

    gdf = {}
    for dat_file in glob.glob(os.path.join(dat_path, "*.dat")):
        # read each layer as individual dataframes into a dictionary
        gdf[os.path.split(dat_file)[1][:-4]] = pd.read_fwf(
            dat_file, header=None, names=["X", "Y", "Z"]
        )

        # assign layer name to a column
        gdf[os.path.split(dat_file)[1][:-4]]["data"] = os.path.split(dat_file)[
            1
        ][:-4]
    dat_file = glob.glob(os.path.join(dat_path, "*.dat"))[-1]

    # find data resolution
    gdf_xr = (
        gdf[os.path.split(dat_file)[1][:-4]].set_index(["X", "Y"]).to_xarray()
    )
    resx = gdf_xr["X"][1] - gdf_xr["X"][0]
    resy = gdf_xr["Y"][1] - gdf_xr["Y"][0]

    # combine dataframes
    gdf = pd.concat(gdf.values())

    # convert dataframe to geodataframe
    gdf["wkt"] = (
        "POINT (" + gdf["X"].astype(str) + " " + gdf["Y"].astype(str) + ")"
    )
    gdf = gpd.GeoDataFrame(
        gdf, geometry=gpd.GeoSeries.from_wkt(gdf["wkt"]), crs=dat_crs
    )
    gdf.drop(columns=["wkt", "X", "Y"], inplace=True)

    # convert to Xarray dataset
    gdf = make_geocube(
        vector_data=gdf,
        resolution=(-abs(resy), abs(resx)),
        align=(abs(resy / 2), abs(resx / 2)),
        group_by="data",
    )

    return gdf

ds = read_dat_file(DATA_DIR, dat_crs=CRS)

ds

<xarray.Dataset> Size: 8MB
Dimensions:      (data: 19, y: 237, x: 218)
Coordinates:
  * y            (y) float64 2kB 5.936e+06 5.936e+06 ... 5.889e+06 5.889e+06
  * x            (x) float64 2kB 6.966e+05 6.968e+05 ... 7.398e+05 7.4e+05
  * data         (data) <U66 5kB 'Flyde Halite Base Depth XYZ Meters' ... 'Ro...
    spatial_ref  int64 8B 0
Data variables:
    Z            (data, y, x) float64 8MB nan nan nan nan ... nan nan nan nan

xarray.Dataset

• Dimensions:
  • data: 19
  • y: 237
  • x: 218
• Coordinates: (4)
  • y
    (y)
    float64
    5.936e+06 5.936e+06 ... 5.889e+06

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([5936400., 5936200., 5936000., ..., 5889600., 5889400., 5889200.])

  • x
    (x)
    float64
    6.966e+05 6.968e+05 ... 7.4e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([696600., 696800., 697000., ..., 739600., 739800., 740000.])

  • data
    (data)
    <U66
    'Flyde Halite Base Depth XYZ Met...

    array(['Flyde Halite Base Depth XYZ Meters',
           'Flyde Halite Thickness - Zone Of Interest - XYZ Meters',
           'Flyde Halite Thickness XYZ Meters',
           'Flyde Halite Top Depth XYZ Meters',
           'Flyde Halite Top TWT XYZ Milliseconds',
           'Mythop Halite Base Depth XYZ Meters',
           'Mythop Halite Thickness XYZ Meters',
           'Mythop Halite Top Depth XYZ Meters',
           'Mythop Halite Top TWT XYZ Milliseconds',
           'Presall Halite Base Depth XYZ Meters-corrected',
           'Presall Halite Thickness - Zone Of Interest - XYZ Meters-corrected',
           'Presall Halite Thickness XYZ Meters-corrected',
           'Presall Halite Top Depth XYZ Meters-corrected',
           'Presall Halite Top TWT XYZ Milliseconds',
           'Rossall Halite Base Depth XYZ Meters',
           'Rossall Halite Thickness - Zone Of Interest - XYZ Meters',
           'Rossall Halite Thickness XYZ Meters',
           'Rossall Halite Top Depth XYZ Meters',
           'Rossall Halite Top TWT XYZ Milliseconds'], dtype='<U66')

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["ED50 / UTM zone 29N",GEOGCS["ED50",DATUM["European_Datum_1950",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],AUTHORITY["EPSG","6230"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4230"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","23029"]]

    semi_major_axis :

    6378388.0

    semi_minor_axis :

    6356911.9461279465

    inverse_flattening :

    297.0

    reference_ellipsoid_name :

    International 1924

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    ED50

    horizontal_datum_name :

    European Datum 1950

    projected_crs_name :

    ED50 / UTM zone 29N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -9.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["ED50 / UTM zone 29N",GEOGCS["ED50",DATUM["European_Datum_1950",SPHEROID["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],AUTHORITY["EPSG","6230"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4230"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","23029"]]

    GeoTransform :

    696500.0 200.0 0.0 5936500.0 0.0 -200.0

    array(0)

• Data variables: (1)
  • Z
    (data, y, x)
    float64
    nan nan nan nan ... nan nan nan nan

    name :

    Z

    long_name :

    Z

    _FillValue :

    nan

    array([[[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
    ...
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]]])

• Indexes: (3)
  • y
    PandasIndex

    PandasIndex(Index([5936400.0, 5936200.0, 5936000.0, 5935800.0, 5935600.0, 5935400.0,
           5935200.0, 5935000.0, 5934800.0, 5934600.0,
           ...
           5891000.0, 5890800.0, 5890600.0, 5890400.0, 5890200.0, 5890000.0,
           5889800.0, 5889600.0, 5889400.0, 5889200.0],
          dtype='float64', name='y', length=237))

  • x
    PandasIndex

    PandasIndex(Index([696600.0, 696800.0, 697000.0, 697200.0, 697400.0, 697600.0, 697800.0,
           698000.0, 698200.0, 698400.0,
           ...
           738200.0, 738400.0, 738600.0, 738800.0, 739000.0, 739200.0, 739400.0,
           739600.0, 739800.0, 740000.0],
          dtype='float64', name='x', length=218))

  • data
    PandasIndex

    PandasIndex(Index(['Flyde Halite Base Depth XYZ Meters',
           'Flyde Halite Thickness - Zone Of Interest - XYZ Meters',
           'Flyde Halite Thickness XYZ Meters',
           'Flyde Halite Top Depth XYZ Meters',
           'Flyde Halite Top TWT XYZ Milliseconds',
           'Mythop Halite Base Depth XYZ Meters',
           'Mythop Halite Thickness XYZ Meters',
           'Mythop Halite Top Depth XYZ Meters',
           'Mythop Halite Top TWT XYZ Milliseconds',
           'Presall Halite Base Depth XYZ Meters-corrected',
           'Presall Halite Thickness - Zone Of Interest - XYZ Meters-corrected',
           'Presall Halite Thickness XYZ Meters-corrected',
           'Presall Halite Top Depth XYZ Meters-corrected',
           'Presall Halite Top TWT XYZ Milliseconds',
           'Rossall Halite Base Depth XYZ Meters',
           'Rossall Halite Thickness - Zone Of Interest - XYZ Meters',
           'Rossall Halite Thickness XYZ Meters',
           'Rossall Halite Top Depth XYZ Meters',
           'Rossall Halite Top TWT XYZ Milliseconds'],
          dtype='object', name='data'))

• Attributes: (0)

ds.rio.crs

CRS.from_epsg(23029)

ds.rio.resolution()

(200.0, -200.0)

ds.rio.bounds()

(696500.0, 5889100.0, 740100.0, 5936500.0)

xmin, ymin, xmax, ymax = extent.total_bounds

def plot_maps(plot_data):
    """
    Create facet contour plots of the Xarray dataset
    """

    fig = plot_data["Z"].plot.contourf(
        col="data",
        cmap="jet",
        col_wrap=2,
        robust=True,
        levels=15,
        subplot_kws={"projection": ccrs.epsg(CRS)},
        xlim=(xmin, xmax),
        ylim=(ymin, ymax),
        cbar_kwargs={"aspect": 20, "pad": 0.02},
    )

    # add a basemap
    basemap = cx.providers.CartoDB.PositronNoLabels
    for n, axis in enumerate(fig.axs.flat):
        cx.add_basemap(axis, crs=CRS, source=basemap, attribution=False)
        # add attribution for basemap tiles
        if n == 2:
            axis.text(xmin, ymin - 2500, basemap["attribution"], fontsize=8)

    # assign titles this way to prevent truncation/overflow
    for axis, title in zip(fig.axs.flat, plot_data["data"].values):
        axis.set_title("\n".join(wrap(title, 34)), fontsize=10)
    plt.show()

Halite thickness

plot_maps(ds.sel(data=[x for x in ds["data"].values if "Thickness XYZ" in x]))

Halite thickness - zones of interest

plot_maps(ds.sel(data=[x for x in ds["data"].values if "Zone" in x]))

Halite base depth

plot_maps(ds.sel(data=[x for x in ds["data"].values if "Base" in x]))

Halite top depth

plot_maps(ds.sel(data=[x for x in ds["data"].values if "Top Depth" in x]))

Halite top TWT (two-way thickness)

plot_maps(ds.sel(data=[x for x in ds["data"].values if "Millisecond" in x]))