# %% # Download static data # -------------------- # This is the first step of the Great Barrier Reef example workflow. It gathers # the static data used by the following scripts: coral-reef polygons, coastline # and land polygons, bathymetric rasters, and a small TPXO-like tidal dataset # with the expected file layout. # # The tidal files created here contain zero amplitudes and phases. They are only # meant to document and test the expected TPXO format; real tidal forcing can be # used by replacing these files with data downloaded from TPXO. import os import urllib.request import zipfile import numpy as np import xarray as xr from gbrtools import poly_to_line, shapefile_to_flat_gpkg from osgeo import gdal from slim4.netcdf_to_abin import tpxo2abin # %% # We define the working directories and the geographic bounding box used to clip # vector datasets to the model region. os.makedirs("tmp", exist_ok=True) os.makedirs("data", exist_ok=True) os.makedirs("data/netcdf/tpxo", exist_ok=True) bounds = [137.0, 159.0, -32.0, -6.0] # %% def extract_progress(file, outpath): """Extract a zip archive while printing a short progress message. Parameters ---------- file : str Path to the zip archive. outpath : str Directory where the archive contents are extracted. """ print(f"Extracting '{file}'... ", end="") with zipfile.ZipFile(file, "r") as zip_ref: zip_ref.extractall(outpath) print("done") # %% # We download coral-reef polygons and convert them to a flattened GeoPackage # clipped to the model bounding box. These polygons are later used both for mesh # refinement and for the reef-density roughness parameterization. if not os.path.exists("tmp/WCMC008_CoralReefs2018_v4_1.zip"): urllib.request.urlretrieve( "https://wcmc.io/WCMC_008", "tmp/WCMC008_CoralReefs2018_v4_1.zip" ) extract_progress("tmp/WCMC008_CoralReefs2018_v4_1.zip", "tmp/reefs") input_shp = "tmp/reefs/14_001_WCMC008_CoralReefs2018_v4_1/01_Data/WCMC008_CoralReef2018_Py_v4_1.shp" output_gpkg = "data/reefs.gpkg" output_layer_name = "reefs" shapefile_to_flat_gpkg(input_shp, output_gpkg, output_layer_name, bounds) # %% # We download the bathymetry datasets covering the Great Barrier Reef, Torres # Strait, and the Gulf of Papua. The original files are stored in # data/bathymetry and reused if already present. # # GBR 30 m bathymetry # Source: https://ecat.ga.gov.au/geonetwork/srv/eng/catalog.search#/metadata/115066 if not os.path.exists("tmp/Great Barrier Reef Bathymetry 2020 30m.zip"): urllib.request.urlretrieve( "https://files.ausseabed.gov.au/survey/Great%20Barrier%20Reef%20Bathymetry%202020%2030m.zip", "tmp/Great Barrier Reef Bathymetry 2020 30m.zip", ) extract_progress("tmp/Great Barrier Reef Bathymetry 2020 30m.zip", "data/bathymetry") # %% # Torres Strait 30 m bathymetry # Source: https://ecat.ga.gov.au/geonetwork/srv/eng/catalog.search#/metadata/144348 if not os.path.exists( "tmp/Australian Bathymetry Topography (Torres Strait) 2023 30m.zip" ): urllib.request.urlretrieve( "https://files.ausseabed.gov.au/survey/Australian%20Bathymetry%20Topography%20(Torres%20Strait)%202023%2030m.zip", "tmp/Australian Bathymetry Topography (Torres Strait) 2023 30m.zip", ) extract_progress( "tmp/Australian Bathymetry Topography (Torres Strait) 2023 30m.zip", "data/bathymetry", ) # %% # GBR 100 m bathymetry # Source: https://ecat.ga.gov.au/geonetwork/srv/eng/catalog.search#/metadata/133163 if not os.path.exists("tmp/Great Barrier Reef Bathymetry 2020 100m.zip"): urllib.request.urlretrieve( "https://files.ausseabed.gov.au/survey/Great%20Barrier%20Reef%20Bathymetry%202020%20100m.zip", "tmp/Great Barrier Reef Bathymetry 2020 100m.zip", ) extract_progress("tmp/Great Barrier Reef Bathymetry 2020 100m.zip", "data/bathymetry") # %% # Gulf of Papua 100 m bathymetry # Source: https://figshare.com/articles/dataset/Gulf_of_Papua_Bathymetry_raster_dataset/11986797/1 if not os.path.exists("data/bathymetry/GulfOfPapua_100m_bathymetry.tif"): urllib.request.urlretrieve( "https://ndownloader.figshare.com/files/22014357", "data/bathymetry/GulfOfPapua_100m_bathymetry.tif", ) # %% # We merge the bathymetry rasters into two products: a 0.001-degree raster used # for the mesh-size field and a coarser 0.01-degree raster useful for fast # prototyping. files = [ "GulfOfPapua_100m_bathymetry.tif", "Great_Barrier_Reef_2020_100m_MSL_cog.tif", "Australian_Bathymetry_Topography__Torres_Strait__2023_30m_MSL_cog.tif", "Great_Barrier_Reef_D_2020_30m_MSL_cog.tif", "Great_Barrier_Reef_C_2020_30m_MSL_cog.tif", "Great_Barrier_Reef_B_2020_30m_MSL_cog.tif", "Great_Barrier_Reef_A_2020_30m_MSL_cog.tif", ] files = ["data/bathymetry/" + f for f in files] warp_options = gdal.WarpOptions( format="GTiff", xRes=0.001, yRes=0.001, outputType=gdal.GDT_Float32, resampleAlg=gdal.GRA_Bilinear, warpMemoryLimit=4000, creationOptions=[ "COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=6", "BIGTIFF=YES", "NUM_THREADS=ALL_CPUS", ], multithread=True, ) os.makedirs("data/bathymetry/merged/", exist_ok=True) ds = gdal.Warp("data/bathymetry/merged/merge_d100.tif", files, options=warp_options) ds = gdal.Warp( "data/bathymetry/merged/merge_d1000.tif", "data/bathymetry/merged/merge_d100.tif", xRes=0.01, yRes=0.01, resampleAlg="average", ) ds = None # %% # We download global land polygons and convert them to coastlines. The coastline # layer is used by mesh.py to close the domain and to control mesh refinement # close to land. if not os.path.exists("tmp/land-polygons-split-4326.zip"): urllib.request.urlretrieve( "https://osmdata.openstreetmap.de/download/land-polygons-split-4326.zip", "tmp/land-polygons-split-4326.zip", ) extract_progress("tmp/land-polygons-split-4326.zip", "tmp/land_polys") input_shp = "tmp/land_polys/land-polygons-split-4326/land_polygons.shp" land_gpkg = "data/contour/land_polygons.gpkg" os.makedirs("data/contour", exist_ok=True) shapefile_to_flat_gpkg(input_shp, land_gpkg, "land", bounds) coastlines_gpkg = "data/contour/coastlines.gpkg" poly_to_line(land_gpkg, coastlines_gpkg, "coast") # %% # We create a dummy TPXO-like dataset with the expected variables, dimensions, # and file names. Elevation and velocity amplitudes/phases are prescribed as # simple functions of lon/lat so that their magnitudes remain roughly realistic. # # This allows the full preprocessing and runtime workflow to be tested without # downloading a real TPXO atlas, while keeping the exact same data path and file # structure as in a production setup. nx, ny, nc = 751, 1141, 2 con_labels = np.array(["m2 ", "s2 "], dtype="S4") lon = np.linspace(135.0, 160.0, nx, dtype=np.float32) lat = np.linspace(-33.0, -5.0, ny, dtype=np.float32) lon_u, lat_u = np.meshgrid(lon, lat, indexing="ij") lon_v, lat_v = np.meshgrid(lon, lat, indexing="ij") lon_z, lat_z = np.meshgrid(lon, lat, indexing="ij") zeros_data = np.zeros((nc, nx, ny), dtype="float32") ha = np.array( [-0.0252 * lon_z - 0.0216 * lat_z + 3.91, -0.0476 * lon_z - 0.0173 * lat_z + 7.18] ) ua = np.array([-3.32 * lon_u - 1.29 * lat_u + 488, -1.62 * lon_u - 0.592 * lat_u + 239]) va = np.array([-1.6 * lon_v - 0.704 * lat_v + 234, -1.77 * lon_v - 0.748 * lat_v + 258]) hp = np.array( [-11.5 * lon_z - 1.26 * lat_z + 1.7e03, -14.2 * lon_z - 5.94 * lat_z + 2.01e03] ) up = np.array( [-25.2 * lon_u - 8.75 * lat_u + 3.77e03, -25.2 * lon_u - 11.3 * lat_u + 3.71e03] ) vp = np.array( [-23.9 * lon_v - 20.9 * lat_v + 3.26e03, -13.5 * lon_v + 2.84 * lat_v + 2.3e03] ) hp = (hp + 180) % 360 - 180 up = (up + 180) % 360 - 180 vp = (vp + 180) % 360 - 180 ds_h = xr.Dataset( data_vars={ "lon_z": ( ("nx", "ny"), lon_z, {"units": "degree_east", "long_name": "longitude of Z nodes"}, ), "lat_z": ( ("nx", "ny"), lat_z, {"units": "degree_north", "long_name": "latitude of Z nodes"}, ), "con": (("nc"), con_labels, {"long_name": "tidal constituents"}), "ha": ( ("nc", "nx", "ny"), ha, { "units": "meter", "long_name": "Tidal elevation amplitude at Z-nodes", "field": "amplitude, scalar", }, ), "hp": ( ("nc", "nx", "ny"), hp, { "units": "degree GMT", "long_name": "Tidal elevation phase at Z-nodes", "field": "phase, scalar", }, ), } ) ds_h.to_netcdf( "data/netcdf/tpxo/h_tpxo_zone.nc", format="NETCDF4", encoding={var: {"_FillValue": np.nan} for var in ds_h.data_vars if var != "con"}, ) ds_uv = xr.Dataset( data_vars={ "lon_u": ( ("nx", "ny"), lon_u, {"units": "degree_east", "long_name": "longitude of U nodes"}, ), "lat_u": ( ("nx", "ny"), lat_u, {"units": "degree_north", "long_name": "latitude of U nodes"}, ), "lon_v": ( ("nx", "ny"), lon_v, {"units": "degree_east", "long_name": "longitude of V nodes"}, ), "lat_v": ( ("nx", "ny"), lat_v, {"units": "degree_north", "long_name": "latitude of V nodes"}, ), "con": (("nc"), con_labels, {"long_name": "tidal constituents"}), "ua": ( ("nc", "nx", "ny"), ua, { "units": "cm/s", "long_name": "Tidal WE velocity amplitude at U-nodes", "field": "vector, W->E", }, ), "Ua": ( ("nc", "nx", "ny"), zeros_data, { "units": "m^2/s", "long_name": "Tidal WE transport amplitude at U-nodes", "field": "vector, W->E", }, ), "up": ( ("nc", "nx", "ny"), up, { "units": "degree GMT", "long_name": "Tidal WE transport/velocity phase at U-nodes", "field": "phase, scalar", }, ), "va": ( ("nc", "nx", "ny"), va, { "units": "cm/s", "long_name": "Tidal SN velocity amplitude at V-nodes", "field": "vector, S->N", }, ), "Va": ( ("nc", "nx", "ny"), zeros_data, { "units": "m^2/s", "long_name": "Tidal SN transport amplitude at V-nodes", "field": "vector, S->N", }, ), "vp": ( ("nc", "nx", "ny"), vp, { "units": "degree GMT", "long_name": "Tidal SN transport/velocity phase at V-nodes", "field": "phase, scalar", }, ), } ) ds_uv.to_netcdf( "data/netcdf/tpxo/u_tpxo_zone.nc", format="NETCDF4", encoding={var: {"_FillValue": np.nan} for var in ds_uv.data_vars if var != "con"}, ) h_path = tpxo2abin("data/netcdf/tpxo/h_tpxo_zone.nc") u_path = tpxo2abin("data/netcdf/tpxo/u_tpxo_zone.nc") # For compatibility with case-insensitive filesystems, SLIM expects transport # variables to be named hua/hva instead of Ua/Va. for varname in ["Ua.abin", "Ua.ajson", "Va.abin", "Va.ajson"]: newname = "h" + varname.lower() os.replace(f"{u_path}/{varname}", f"{u_path}/{newname}")