.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/example_cross.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_example_cross.py: Plot cross-section ================== .. GENERATED FROM PYTHON SOURCE LINES 8-16 In this notebook, we will see how a cross-section can be created, and how we can extract temperature values from Norkyst along this section. **Python requirements** To access data from the model and extracting it into datasets we will make use of some Python packages. Xarray will be the main tool to opening the datasets, and allows us to display the contents nicely. Cartopy and matplotlib are the main plotting tools, in addition to Cmocean for colormaps. Additionally, we will use pyresample and pyproj to deal with geographical coordinates and their collocation with the model grid. .. GENERATED FROM PYTHON SOURCE LINES 16-31 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import BoundaryNorm, ListedColormap import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import cmocean.cm as cmo import xarray as xr from pyproj import Geod import pyresample .. GENERATED FROM PYTHON SOURCE LINES 32-33 We define the path for the target data set here. .. GENERATED FROM PYTHON SOURCE LINES 33-45 .. code-block:: Python YEAR = 2026 MT = str(4).zfill(2) # zero-padded DAY = str(25).zfill(2) path = 'https://thredds.met.no/thredds/dodsC/fou-hi/norkystv3_his_files/%i/%s/%s/norkyst800_his_zdepth_%i%s%sT00Z_m00_AN.nc'%(YEAR, MT, DAY, YEAR, MT, DAY) ds = xr.open_dataset(path) "" ds .. raw:: html 
<xarray.Dataset> Size: 28GB
    Dimensions:                  (time: 24, depth: 15, Y: 1148, X: 2747)
    Coordinates:
      * time                     (time) datetime64[ns] 192B 2026-04-25 ... 2026-0...
      * depth                    (depth) float64 120B 0.0 1.0 2.0 ... 200.0 300.0
      * Y                        (Y) int32 5kB 0 800 1600 ... 916000 916800 917600
      * X                        (X) int32 11kB 0 800 1600 ... 2196000 2196800
        lat                      (Y, X) float64 25MB ...
        lon                      (Y, X) float64 25MB ...
    Data variables:
        forecast_reference_time  datetime64[ns] 8B ...
        projection_stere         int32 4B ...
        AKs                      (time, depth, Y, X) float32 5GB ...
        Uwind_eastward           (time, Y, X) float32 303MB ...
        Vwind_northward          (time, Y, X) float32 303MB ...
        h                        (Y, X) float64 25MB ...
        salinity                 (time, depth, Y, X) float32 5GB ...
        temperature              (time, depth, Y, X) float32 5GB ...
        u_eastward               (time, depth, Y, X) float32 5GB ...
        v_northward              (time, depth, Y, X) float32 5GB ...
        w                        (time, depth, Y, X) float32 5GB ...
        zeta                     (time, Y, X) float32 303MB ...
    Attributes: (12/43)
        id:                              7c375a70-4ee2-4e3d-b24a-0877ac760a4a
        naming_authority:                no.met
        operational_status:              Operational
        iso_topic_category:              oceans
        activity_type:                   Numerical Simulation
        keywords_vocabulary:             GCMDSK:GCMD Science Keywords:https://gcm...
        ...                              ...
        contributor_name:                Magne Simonsen, Mateusz Matuszak
        contributor_role:                Technical contact, Metadata author
        contributor_email:               magnes@met.no, mateuszm@met.no
        history:                         Sat Apr 25 18:41:49 2026: ncks -O -7 -L ...
        NCO:                             netCDF Operators version 4.8.1 (Homepage...
        DODS_EXTRA.Unlimited_Dimension:  time


.. GENERATED FROM PYTHON SOURCE LINES 46-47 We select now the transect points we want to work on, and make a quick plot of it over the surface temperature values .. GENERATED FROM PYTHON SOURCE LINES 47-59 .. code-block:: Python # The transect we will be working with lat_sec = np.array([57.95, 57.10]) # start_lat, end_lat lon_sec = np.array([7.31, 8.41]) # start_lon, end_lon #lat_sec = np.array([58.50, 57.68]) # start_lat, end_lat #lon_sec = np.array([9.02, 10.32]) # start_lon, end_lon # Latitude and longitude of the dataset lat = ds.lat.values lon = ds.lon.values .. GENERATED FROM PYTHON SOURCE LINES 60-62 ## Make a map of the area with the transect ### ############################################### .. GENERATED FROM PYTHON SOURCE LINES 62-105 .. code-block:: Python # Choosing a map projection proj = ccrs.NorthPolarStereo() # Making figure and axes with the projection fig, ax = plt.subplots(subplot_kw={'projection': proj}, constrained_layout=True, dpi=200) # Setting the extent of the map to our model domain ax.set_extent([np.min(lon), np.max(lon), np.min(lat), np.max(lat)], crs=ccrs.PlateCarree()) # ccrs.PlateCarree() to tell the program our data is in coordinates lats/lons # Adding natural features to our map land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', edgecolor='gray', facecolor=cfeature.COLORS['land']) coastline = cfeature.NaturalEarthFeature(category='physical', name='coastline', scale='50m', edgecolor='gray', facecolor='none') borders = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_boundary_lines_land', edgecolor= 'gray', scale='50m', facecolor='none') ax.add_feature(land) ax.add_feature(coastline) ax.add_feature(borders) # Adding gridlines gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='lightgray', alpha=0.5, linestyle='--') gl.top_labels = False # Disable top labels gl.right_labels = False # Disable right labels gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER # Plotting boundaries of model ax.plot(lon[0,:], lat[0,:], '--', transform= ccrs.PlateCarree(), color = 'gray', linewidth =0.8) ax.plot(lon[-1,:], lat[-1,:], '--', transform= ccrs.PlateCarree(), color = 'gray', linewidth =0.8) ax.plot(lon[:,0], lat[:,0], '--', transform= ccrs.PlateCarree(), color = 'gray', linewidth =0.8) ax.plot(lon[:,-1], lat[:,-1], '--', transform= ccrs.PlateCarree(), color = 'gray', linewidth =0.8, label='Norkyst boundary') # Plotting sea surface height temp_norkyst = ds.temperature.isel(time=22, depth=0) cs = temp_norkyst.plot.pcolormesh(ax=ax, x='lon', y='lat', vmin=6, vmax=10, cmap=cmo.thermal, transform=ccrs.PlateCarree()) # Adding transect ax.plot([lon_sec[0], lon_sec[1]], [lat_sec[0], lat_sec[1]], transform = ccrs.PlateCarree(), color = 'r', label = 'Transect') # Adding legend ax.legend(loc = 'upper left') .. image-sg:: /gallery/images/sphx_glr_example_cross_001.png :alt: depth = 0.0 [meter], time = 2026-04-25T22:00:00 :srcset: /gallery/images/sphx_glr_example_cross_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 106-107 The next step now is the extraction of values along the defined transect. We want salinity and temperature values. For that, we have to define a function which will perform the job of `xroms.argsel2d`. .. GENERATED FROM PYTHON SOURCE LINES 107-173 .. code-block:: Python def generate_transect_points(lat_init, lon_init, lat_end, lon_end, n_points=100): """ Generate evenly spaced points along a geodesic transect. Parameters ---------- lat_init, lon_init : float Start point coordinates. lat_end, lon_end : float End point coordinates. n_points : int Total number of points INCLUDING endpoints. Returns ------- lats : ndarray lons : ndarray """ geod = Geod(ellps="WGS84") # npts excludes endpoints intermediate = geod.npts( lon_init, lat_init, lon_end, lat_end, n_points - 2 ) # Build full coordinate list lons = [lon_init] + [p[0] for p in intermediate] + [lon_end] lats = [lat_init] + [p[1] for p in intermediate] + [lat_end] return np.array(lats), np.array(lons) # Here we generate the transect points using the function defined above lat_t, lon_t = generate_transect_points(lat_sec[0], lon_sec[0], lat_sec[1], lon_sec[1], n_points=100) # Extract the 2D longitude and latitude arrays from the dataset lon2d = ds['lon'].values lat2d = ds['lat'].values ds_geo = pyresample.geometry.GridDefinition(lons=lon2d, lats=lat2d) pos_geo = pyresample.geometry.SwathDefinition(lons=lon_t, lats=lat_t) _, valid_output_index, index_array, distance_array = \ pyresample.kd_tree.get_neighbour_info( source_geo_def=ds_geo, target_geo_def=pos_geo, radius_of_influence=800, neighbours=1) index_array_2d = np.unravel_index(index_array, ds.lon.shape) # Ensure rows and cols are treated as 1D arrays rows_da = xr.DataArray(np.array(index_array_2d[0]), dims="points") cols_da = xr.DataArray(np.array(index_array_2d[1]), dims="points") # Extract the data for these unique indices transect_ds = ds.isel(Y=rows_da, X=cols_da) # Replace 'y' and 'x' with actual dimension names .. GENERATED FROM PYTHON SOURCE LINES 174-175 We use a function here to calculate the cumulative distance along the transect .. GENERATED FROM PYTHON SOURCE LINES 175-209 .. code-block:: Python def compute_transect_distance(lons, lats): """ Compute cumulative distance along transect. Parameters ---------- lons, lats : 1D arrays Returns ------- distance_km : 1D array Cumulative distance in km """ geod = Geod(ellps="WGS84") distance_km = np.zeros(len(lons)) for i in range(1, len(lons)): _, _, dist_m = geod.inv( lons[i-1], lats[i-1], lons[i], lats[i] ) distance_km[i] = distance_km[i-1] + dist_m / 1000 return distance_km distance_km = compute_transect_distance(lon_t, lat_t) "" .. rst-class:: sphx-glr-script-out .. code-block:: none '' .. GENERATED FROM PYTHON SOURCE LINES 210-212 We plot our transaction using a grid interpolator, so the bathymetry contour is properly displayed. .. GENERATED FROM PYTHON SOURCE LINES 212-288 .. code-block:: Python from scipy.interpolate import RegularGridInterpolator # ----------------------- # INPUT DATA # ----------------------- temp = transect_ds.temperature[0,:,:].values # (depth, points) z = transect_ds.depth.values # (depth,) dist = distance_km # (points,) h = transect_ds.h.values # (points,) # ----------------------- # 1. CREATE TARGET GRID # ----------------------- dist_grid = np.linspace(dist.min(), dist.max(), 300) z_grid = np.linspace(0, z.max(), 200) dist2d, z2d = np.meshgrid(dist_grid, z_grid) # ----------------------- # 2. INTERPOLATOR (structured grid: depth × distance) # ----------------------- f = RegularGridInterpolator( (z, dist), temp, bounds_error=False, fill_value=np.nan ) # evaluate on grid points = np.column_stack([z2d.ravel(), dist2d.ravel()]) temp_grid = f(points).reshape(z2d.shape) # ----------------------- # 3. BATHYMETRY MASK # ----------------------- # interpolate h onto grid horizontally h_interp = np.interp(dist_grid, dist, h) mask = z2d <= h_interp[None, :] temp_grid[~mask] = np.nan # ----------------------- # 4. PLOT # ----------------------- fig, ax = plt.subplots(figsize=(12, 5)) # seabed ax.fill_between( dist_grid, z_grid.min(), h_interp, color="w", alpha=0.3, ec='k', ) pcm = ax.pcolormesh( dist_grid, z_grid, temp_grid, vmin=6, vmax=8, cmap=cmo.thermal, shading="auto" ) cb = plt.colorbar(pcm, ax=ax) cb.set_label(r"Potential temperature [C$^{\degree}$]") ax.set_xlabel("Distance (km)") ax.set_ylabel("Depth (m)") ax.invert_yaxis() ax.set_title("Interpolated cross-section") .. image-sg:: /gallery/images/sphx_glr_example_cross_002.png :alt: Interpolated cross-section :srcset: /gallery/images/sphx_glr_example_cross_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Interpolated cross-section') .. GENERATED FROM PYTHON SOURCE LINES 289-298 Cross-section in S-coordinates ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We can also make a plot of the same cross-section using outputs defined in S-coords. The later represents a terrain-following coordinate system, therefore, we must convert the output variables there to true depth values using the following equation: $$Z_0 = \frac{h_c \, S + h \, C}{h_c + h}$$ $$z = Z_0 (\zeta + h) + \zeta$$ .. GENERATED FROM PYTHON SOURCE LINES 298-305 .. code-block:: Python path_s = './datasets/transect_scoord.nc' ds_s = xr.open_dataset(path_s, engine='netcdf4') "" .. rst-class:: sphx-glr-script-out .. code-block:: none '' .. GENERATED FROM PYTHON SOURCE LINES 306-307 Employ the conversion here: .. GENERATED FROM PYTHON SOURCE LINES 307-313 .. code-block:: Python Zo_rho = (ds_s.hc * ds_s.s_rho + ds_s.Cs_r * ds_s.h) / (ds_s.hc + ds_s.h) z_rho = ds_s.zeta + (ds_s.zeta + ds_s.h) * Zo_rho ds_s.coords["z_rho"] = z_rho.transpose() .. GENERATED FROM PYTHON SOURCE LINES 314-315 You can have a quick look at the results using: .. GENERATED FROM PYTHON SOURCE LINES 315-320 .. code-block:: Python section = ds_s.temp section.plot(x="lon_rho", y="z_rho", figsize=(15, 6), clim=(6, 8), cmap=cmo.thermal) plt.ylim([-500, 1]); .. image-sg:: /gallery/images/sphx_glr_example_cross_003.png :alt: ocean_time = 2026-04-25T22:00:00 :srcset: /gallery/images/sphx_glr_example_cross_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (-500.0, 1.0) .. GENERATED FROM PYTHON SOURCE LINES 321-322 Or, plot it using Matplotlib if you want more control over your options: .. GENERATED FROM PYTHON SOURCE LINES 322-440 .. code-block:: Python from scipy.interpolate import griddata # ----------------------------------- # SELECT TIME # ----------------------------------- temp = ds_s.temp.values # shape: (40, 100) # ----------------------------------- # DISTANCE # ----------------------------------- dist = distance_km # shape: (100,) # ----------------------------------- # TRUE DEPTHS # ----------------------------------- z_rho = ds_s.z_rho.values # ----------------------------------- # BUILD SCATTERED POINTS # ----------------------------------- dist2d = np.tile(dist, (temp.shape[0], 1)) points = np.column_stack([ dist2d.ravel(), z_rho.ravel() ]) values = temp.ravel() # remove NaNs mask = np.isfinite(values) & np.isfinite(points[:,1]) points = points[mask] values = values[mask] # ----------------------------------- # TARGET GRID # ----------------------------------- dist_grid = np.linspace(dist.min(), dist.max(), 400) z_grid = np.linspace( np.nanmin(z_rho), 0, 300 ) distg, zg = np.meshgrid(dist_grid, z_grid) # ----------------------------------- # INTERPOLATE # ----------------------------------- temp_grid = griddata( points, values, (distg, zg), method="linear" ) # ----------------------------------- # BATHYMETRY # ----------------------------------- h = -transect_ds.h.values h_grid = np.interp( dist_grid, dist, h ) # mask below seafloor temp_grid[zg < h_grid[None, :]] = np.nan # ----------------------------------- # PLOT # ----------------------------------- fig, ax = plt.subplots(figsize=(12,5)) pcm = ax.contourf( distg, zg, temp_grid, levels=np.linspace(6, 8, 41), shading="auto", cmap=cmo.thermal, extend="both" ) # Define colorbar cbar = plt.colorbar(pcm, ax=ax, label=r"Potential temperature [C$^{\degree}$]") cbar.ax.locator_params(nbins=5) # fill land ax.fill_between( dist_grid, h_grid, zg.min()-20, color="k", alpha=0.3 ) # Set labels ax.set_ylim(-480, 0) ax.set_xlabel("Distance (km)") ax.set_ylabel("Depth (m)") plt.tight_layout() plt.show() .. image-sg:: /gallery/images/sphx_glr_example_cross_004.png :alt: example cross :srcset: /gallery/images/sphx_glr_example_cross_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/norkyst.github.io/norkyst.github.io/examples/example_cross.py:409: UserWarning: The following kwargs were not used by contour: 'shading' pcm = ax.contourf( .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 39.863 seconds) .. _sphx_glr_download_gallery_example_cross.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_cross.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_cross.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_cross.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_