From 537620cf7bcec9c5ad15d75cb146ae0a6da7c01b Mon Sep 17 00:00:00 2001 From: "Bauer, Leonie" <leonie.bauer-1@studium.uni-hamburg.de> Date: Tue, 9 Jan 2024 12:21:52 +0000 Subject: [PATCH] Add new file --- Messmethoden_U6/U6_leonie | 341 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 341 insertions(+) create mode 100644 Messmethoden_U6/U6_leonie diff --git a/Messmethoden_U6/U6_leonie b/Messmethoden_U6/U6_leonie new file mode 100644 index 0000000..9ffff99 --- /dev/null +++ b/Messmethoden_U6/U6_leonie @@ -0,0 +1,341 @@ +# Importing the required module for plotting data +import gsw # See https://teos-10.github.io/GSW-Python/ +import pycnv # See https://pypi.org/project/pycnv/ +import matplotlib.pyplot as plt +import pylab as pl + +# Defining the file path where the CNV file is located +# Choose either Seepraktikum 2023 data: https://bit.ly/3ut6NtV +#file_path = 'Messmethoden/Seepraktikum_2023/SBE19plus_01907321_2023_05_12_Cast36_loop_window.cnv' +# Or MSM121 data: https://bit.ly/3Ri19E2 +#file_path = 'Messmethoden/msm121_profiles/MSM121_015_1db.cnv' +# Or MSM121 to-yo profiles: https://bit.ly/3GlQh1D +file_path = 'Messmethoden/MSM121_060/MSM121_060_001_05db.cnv' + +# Load the data from the given file_path (may need to update to match your file location) +cnv = pycnv.pycnv(file_path) + +# Print some info to the screen +print('Test if we are in the Baltic Sea (usage of different equation of state): ' + str(cnv.baltic)) +print('Position of cast is: Longitude:', cnv.lon,'Latitude:',cnv.lat) +print('Time of cast was:', cnv.date) +print('Number of sensor entries (len(cnv.data.keys())):',len(cnv.data.keys())) +print('Names of sensor entries (cnv.data.keys()):',cnv.data.keys()) + +# Get data of entry +key0 = list(cnv.data.keys())[0] +data0 = cnv.data[key0] + +# Get derived data: +keyd0 = list(cnv.cdata.keys())[0] +datad0 = cnv.cdata[keyd0] +# Get unit of derived data +datad0_unit = cnv.cunits[keyd0] + +# Standard names are mapped to +# cnv.p,cnv.CT,cnv.T,cnv.SP,cnv.oxy +# units are _unit, e.g. cnv.p_unit + +print(list(cnv.data.keys())) +print(list(cnv.cdata.keys())) + +# Plot standard parameters +pl.figure(1) +pl.clf() +pl.subplot(1,2,1) +pl.plot(cnv.SA,cnv.p) # Note that pycnv has done the TEOS-10 conversion to SA for you using gsw +pl.xlabel('Absolute salinity [' + cnv.SA_unit + ']') +pl.ylabel('Pressure [' + cnv.p_unit + ']') +pl.gca().invert_yaxis() # Question: What happens if you comment out this line? +pl.grid(color='grey', linestyle='-', linewidth=0.5) + + +# Step 1a: Add grid lines +pl.grid(color='grey', linestyle='-', linewidth=0.5) + +# Step 1b: Limit the top of the plot to the surface (p=0) +pl.ylim(1600, 0) + +# Step 1c: Add a second plot for temperature, to the right of the salinity plot +pl.subplot(1,2,2) +pl.plot(cnv.T,cnv.p, color='red') # Note that pycnv has done the TEOS-10 conversion to SA for you using gsw +pl.xlabel('Temperature [' + cnv.T_unit + ']') +pl.ylabel('Pressure [' + cnv.p_unit + ']') +pl.gca().invert_yaxis() +pl.grid(color='grey', linestyle='-', linewidth=0.5) +pl.subplots_adjust(right=1) +pl.subplots_adjust(left=0.1, + bottom=0.1, + right=0.9, + top=0.9, + wspace=0.4, + hspace=0.4) #https://www.geeksforgeeks.org/how-to-set-the-spacing-between-subplots-in-matplotlib-in-python/ + +# Step 1d: Add a title to your figure, perhaps the station number, latitude and longitude +latitude = round(cnv.lat, 2) +longitude = round(cnv.lon, 2) #https://stackoverflow.com/questions/20457038/how-to-round-to-2-decimals-with-python +pl.suptitle('Station: 64,' +' Latitude: '+str(latitude)+', Longitude: '+str(longitude)) #title function needs strings +#one title for all subplots: https://stackoverflow.com/questions/7066121/how-to-set-a-single-main-title-above-all-the-subplots + + +#see which attributes are in the file (https://code.activestate.com/pypm/cnv/) +#file = 'Messmethoden/MSM121_060/MSM121_060_001_05db.cnv' +#cnvdump file.cnv +#cnv2cdf file.cnv +#from cnv.fCNV import fCNV +#profile = fCNV('file') +#profile.attributes + + +# Step 1e: Print the figure to a *png file +import os +desktop_pfad = os.path.join(os.path.expanduser('~'), 'Desktop') #zum Desktop +ordner_pfad = os.path.join(desktop_pfad, 'UNI', 'Semester5', 'Messmethoden' , 'Excercise_6') #zum Ordner dort +bild_pfad = os.path.join(ordner_pfad, '1_stat64_salt_temperature.png') +plt.savefig(bild_pfad, bbox_inches='tight') + + +# Plotting profile data on the same axes +pl.subplot(1,2,2) +pl.plot(cnv.oxy,cnv.p) +pl.plot(cnv.cdata['oxy0'],cnv.p, label='Oxy0') +pl.plot(cnv.cdata['oxy1'],cnv.p, label='Oxy1') +pl.legend() +pl.xlabel('Oxygen [' + cnv.oxy_unit + ']') +pl.ylabel('Pressure [' + cnv.p_unit + ']') +pl.gca().invert_yaxis() +pl.title('Station: 64,' +' Latitude: '+str(latitude)+', Longitude: '+str(longitude)) #title function needs strings +bild_pfad = os.path.join(ordner_pfad, '2_stat64_oxygen_oneplot.png') +plt.savefig(bild_pfad, bbox_inches='tight') + + +# Step 2a-e: Clean up the figure as above +pl.subplot(1,2,1) +pl.plot(cnv.oxy,cnv.p) +pl.plot(cnv.cdata['oxy0'],cnv.p, color='orange', label='Oxy0') +pl.xlabel('Oxygen [' + cnv.oxy_unit + ']') +pl.ylabel('Pressure [' + cnv.p_unit + ']') +pl.grid(color='grey', linestyle='-', linewidth=0.5) #grid +pl.gca().invert_yaxis() + + +pl.subplot(1,2,2) +pl.plot(cnv.cdata['oxy1'],cnv.p, color='green', label='Oxy1') +pl.xlabel('Oxygen [' + cnv.oxy_unit + ']') +pl.ylabel('Pressure [' + cnv.p_unit + ']') +pl.grid(color='grey', linestyle='-', linewidth=0.5) +pl.gca().invert_yaxis() +pl.subplots_adjust(right=1) +pl.subplots_adjust(left=0.1, + bottom=0.1, + right=0.9, + top=0.9, + wspace=0.4, + hspace=0.4) + +pl.suptitle('Station: 64,' +' Latitude: '+str(latitude)+', Longitude: '+str(longitude)) #title function needs strings + + +#savefig +bild_pfad = os.path.join(ordner_pfad, '3_stat64_oxygen_subplots.png') +plt.savefig(bild_pfad, bbox_inches='tight') + + + +#Step 2f (optional): Compute the apparent oxygen utilisation See: https://en.wikipedia.org/wiki/#Apparent_oxygen_utilisation you will need to calculate oxygen solubility (see the gsw documentation for O2_sol) + +#AOU = oxygen solubility - oxygen observed + +#solubility (Löslichkeit) +#observed (gemessen) + + + +#calculate solubility: +#https://www.teos-10.org/pubs/gsw/html/gsw_O2sol.html +import gsw +#from gsw import gsw_O2sol +#help(gsw_O2sol) + +SA = cnv.SA +CT = cnv.CT +p = cnv.p +long = cnv.lon +lat = cnv.lat + + +print('length SA:',len(SA)) +print('shape CT:', CT.shape) + +print() +print('Absolute Salinity:', SA) +print() +print('Conservative temperature:', CT) +print() +print('Pressure:', p) +print() +#print(long) +#print(lat) +print() + +#lat und lon into numpy arrays: https://stackoverflow.com/questions/5891410/numpy-array-initialization-fill-with-identical-values +import numpy as np +lat = np.full((5022,), str(lat)) +long = np.full((5022,), str(long)) +print('Shape Lat_array:', lat.shape) +print('Shape Lon_array:', long.shape) +print() + +def y(SA,CT,p,long,lat): + return gsw.O2sol(SA,CT,p,long,lat) + + +oxygen_solubility = gsw.O2sol(SA,CT,p,long,lat) + +print('Oxygen solobility [umol/kg]:', oxygen_solubility) #Einheit: https://teos-10.github.io/GSW-Python/gsw_flat.html# +print() + + +#AOU = oxygen solubility - oxygen observed + +oxygen_observed = cnv.oxy +print('Observed Oxygen:', '[',cnv.oxy_unit,']' , oxygen_observed) #Einheit nachgelesen eigentlich ml/l +print() + +AOU = oxygen_solubility - oxygen_observed +print('Apparent oxygen utilisation:', AOU) + + +# Make a T-S diagram +pl.clf() +pl.plot(cnv.SA,cnv.CT) +pl.xlabel('Absolute salinity [' + cnv.SA_unit + ']') +pl.ylabel('Conservative temperature [' + cnv.CT_unit + ']') + +# Step 3a: Add contours of potential density +# - Hint, you will need to create a vector of temperature and salinity spanning at +# least the ranges in the profile +# - You will then need to use these gridded vectors to calculate potential density using the gsw toolbox +print(cnv.SA.shape) +print() +#source: https://jakevdp.github.io/PythonDataScienceHandbook/04.04-density-and-contour-plots.html +x = cnv.SA +y = cnv.CT +X, Y = np.meshgrid(x, y)#builds two-dimensional grids from one-dimensional arrays (T+S to T&S) + +#Z as potentail desity +#define function for calculating potential density from SA & CT (pre-defined gsw function) +def f(x,y): + return gsw.density.sigma0(x, y) + +#Z as potential density calculated from temperature and salinity two dimensional grids, type(Z):numpy ndarray, Z.shape:(5022, 5022) +Z = f(X,Y) +plot = plt.contour(X, Y, Z, colors='black') + +#potential density in new variable, type(pot_dens): numpy ndarray, pot_dens.shape: p22,) +pot_dens = gsw.density.sigma0(x, y) +print('potential density from gsw function:', pot_dens) +print() + + +# Step 3b: Add contour labels for the density + +#countour adding (source: https://matplotlib.org/stable/gallery/images_contours_and_fields/contour_label_demo.html) +format = {} +labels = pot_dens.tolist() +for lev, lab in zip(plot.levels, labels): + format[lev]=lab + +plt.clabel(plot, plot.levels, inline=True, fmt=format, fontsize=10) + + + + +# Step 3c: Add a title +pl.suptitle('Station: 64,' +' Latitude: '+str(latitude)+', Longitude: '+str(longitude)) #title function needs strings + + + +# Step 3d: Print the figure as a *png +bild_pfad = os.path.join(ordner_pfad, '4_stat64_T_S_diagram.png') +plt.savefig(bild_pfad, bbox_inches='tight') + + + +# Challenge 4: Load multiple cnv files into python. +cnv1 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_001_1db.cnv') +cnv2 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_002_1db.cnv') +cnv3 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_003_1db.cnv') +cnv4 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_014_1db.cnv') +cnv5 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_015_1db.cnv') +cnv6 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_016_1db.cnv') +cnv7 = pycnv.pycnv('Messmethoden/msm121_profiles/MSM121_017_1db.cnv') + +# Step 4a: Plot a map of the stations + +#first: values for station, Lon and lat (https://stackoverflow.com/questions/53233228/plot-latitude-longitude-from-csv-in-python-3-6) +stations = ['1', '2', '3', '14', '15', '16', '17'] +latitudes = [cnv1.lat, cnv2.lat, cnv3.lat, cnv4.lat, cnv5.lat, cnv6.lat, cnv7.lat] +longitudes = [cnv1.lon, cnv2.lon, cnv3.lon, cnv4.lon, cnv5.lon, cnv6.lon, cnv7.lon] + +#second: data to csv (https://www.delftstack.com/de/howto/python/write-list-to-csv-python/) +import pandas as pd +table = {"station": stations, "latitude": latitudes, "longitude": longitudes} +df = pd.DataFrame(table) +df.to_csv("lon_lat.csv") + +#third: plot coordinates on geopandas map (https://stackoverflow.com/questions/53233228/plot-latitude-longitude-from-csv-in-python-3-6) +import geopandas as gpd +from geopandas import GeoDataFrame +import shapely +from shapely.geometry import Point +#define dataframe +df = pd.read_csv('lon_lat.csv', delimiter=',', skiprows=0, low_memory=False) +#define datapoints for map +geometry = [Point(xy) for xy in zip(df['longitude'], df['latitude'])] +#plot as geodataframe? +gdf = GeoDataFrame(df, geometry=geometry) +#load worldmap +world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres')) +#use gdf (defined geodatframe with coordinates) abd plot onto map +gdf.plot(ax=world.plot(figsize=(10, 6)), marker='o', color='red', markersize=15); +bild_pfad = os.path.join(ordner_pfad, '5_station_map.png') +plt.savefig(bild_pfad, bbox_inches='tight') +plt.show() + +#again, but zoomed in +#use gdf (defined geodatframe with coordinates) abd plot onto map +gdf.plot(ax=world.plot(figsize=(10, 6)), marker='o', color='red', markersize=15); +plt.ylim(25,75) +plt.xlim(-100,0) +bild_pfad = os.path.join(ordner_pfad, '6_station_map_zoomed.png') +plt.savefig(bild_pfad, bbox_inches='tight') + +#again, but zoomed in more +#use gdf (defined geodatframe with coordinates) abd plot onto map +gdf.plot(ax=world.plot(figsize=(10, 6)), marker='o', color='red', markersize=15); +plt.ylim(35,55) +plt.xlim(-60,-30) +bild_pfad = os.path.join(ordner_pfad, '7_station_map_zoomed_more.png') +plt.savefig(bild_pfad, bbox_inches='tight') + + +# Step 4b (optional): Add bathymetry +# - see e.g. GEBCO https://www.gebco.net/data_and_products/gridded_bathymetry_data/ +# - or ETOPO1. Note these files are big when full resolution and global. +# You will probably want to slice/prepare the data in a separate python notebook/file + +# Step 4c: Load multiple stations into python +# - Option one: Keep them all as separate variables (cnv1, cnv2, cnv3, etc) +# - Option two: Combine them into a multidimensaionl np.ndarray +# In order to combine them, they will need to be interpolated onto the same +# pressure grid (suggested 1 dbar grid). + +# Step 4d: Plot multiple temperature profiles on a single set of axes (see Step 1) + +# Step 4e: Repeat for salinity + +# Step 4f: Repeat for T-S diagrams (see Step 3) + +# Noting the data type that cnv is using may help: +type(cnv.SA) -- GitLab