From 8745f85dcab07a5ac8e0818ac3150234fe4d3bd7 Mon Sep 17 00:00:00 2001 From: Alexander Rey Date: Thu, 29 Jun 2023 16:57:42 -0400 Subject: [PATCH] Final Baird Commit 2 --- Azure/hydroUsage.py | 220 ++++++ Azure/secret.py | 2 + EWR_D3D/.idea/EWR_D3D.iml | 2 +- EWR_D3D/.idea/misc.xml | 2 +- EWR_D3D/EWR_PlotHD_D3D.py | 5 +- EWR_D3D/EWR_PlotHD_D3D_Profile.py | 12 +- EWR_D3D/EWR_Plotting.py | 4 +- EWR_D3D/WSP_DB.py | 90 +++ EWR_Data/EWR_Hg_Interp.py | 2 + EWR_Data/EWR_Hg_Interp_RevB.py | 59 +- EWR_Data/EWR_Hg_Syn.py | 213 +----- EWR_Data/createSeg.py | 33 + EWR_Data/output_seg.seg | 1 + EWR_Data/output_seg10.bin | Bin 0 -> 240 bytes EWR_Data/output_seg2.seg | Bin 0 -> 200 bytes EWR_Data/output_seg3.seg | Bin 0 -> 240 bytes EWR_Data/output_seg4.seg | Bin 0 -> 240 bytes EWR_Data/output_seg5.seg | Bin 0 -> 240 bytes EWR_Data/output_seg6.seg | Bin 0 -> 240 bytes EWR_Data/output_seg7.seg | Bin 0 -> 240 bytes EWR_Data/output_seg8.bin | Bin 0 -> 220 bytes EWR_Data/output_seg9.bin | Bin 0 -> 240 bytes LSU_D3D/PenobscotExample.py | 67 ++ Mustique/WestCoastDataTemplate_V5.py | 704 +++++++++++++----- NTC_DFM/Delft3DBoundsCompare.py | 83 +++ NTC_DFM/NTC_PlottingD3D_2023.py | 2 +- NTC_DFM/NTC_SedTestPlot.py | 45 +- NTC_DFM/NTC_VelAnimation.py | 433 ++++------- ...odel_Temp_Sal_EFDC_Delft3D_TM_2014-2015.py | 153 ++-- .../inspectionProfiles/Project_Default.xml | 13 + .../inspectionProfiles/profiles_settings.xml | 6 + SatBathy/.idea/misc.xml | 4 + SatBathy/.idea/modules.xml | 8 + SatBathy/.idea/vcs.xml | 6 + SatBathy/.idea/workspace.xml | 62 ++ 35 files changed, 1466 insertions(+), 765 deletions(-) create mode 100644 Azure/hydroUsage.py create mode 100644 EWR_Data/output_seg.seg create mode 100644 EWR_Data/output_seg10.bin create mode 100644 EWR_Data/output_seg2.seg create mode 100644 EWR_Data/output_seg3.seg create mode 100644 EWR_Data/output_seg4.seg create mode 100644 EWR_Data/output_seg5.seg create mode 100644 EWR_Data/output_seg6.seg create mode 100644 EWR_Data/output_seg7.seg create mode 100644 EWR_Data/output_seg8.bin create mode 100644 EWR_Data/output_seg9.bin create mode 100644 LSU_D3D/PenobscotExample.py create mode 100644 SatBathy/.idea/inspectionProfiles/Project_Default.xml create mode 100644 SatBathy/.idea/inspectionProfiles/profiles_settings.xml create mode 100644 SatBathy/.idea/misc.xml create mode 100644 SatBathy/.idea/modules.xml create mode 100644 SatBathy/.idea/vcs.xml create mode 100644 SatBathy/.idea/workspace.xml diff --git a/Azure/hydroUsage.py b/Azure/hydroUsage.py new file mode 100644 index 0000000..341ff74 --- /dev/null +++ b/Azure/hydroUsage.py @@ -0,0 +1,220 @@ +import requests +import json +import time +import xmltodict +import datetime +from bs4 import BeautifulSoup +import secret + +loginUrl = ('https://api.loginradius.com/identity/v2/auth/login?apiKey=d842e884-2dfb-4c8f-a971-f9eacf8e9f54&loginUrl=&emailTemplate=Verification English&verificationUrl=https://account.hydroottawa.com/login') +authenticatorUrl = ('https://account.hydroottawa.com/ajax/authenticator') +secureUrl = ('https://account.hydroottawa.com/') +downloadDataUrl = ('https://secure.hydroottawa.com/Usage/Secure/TOU/DownloadMyData.aspx') +ssoMhlUrl = ('https://account.hydroottawa.com/ajax/sso_mhl') +ssoTargetUrl = ('https://secure.hydroottawa.com/Usage/SSO/SSOTarget.aspx') + +jsonPayload = { + 'email': secret.email, + 'password': secret.password +} +headersPayloadOptions = { + 'referrer': secureUrl, + 'origin': secureUrl, + 'Access-Control-Request-Method': 'POST', + 'Connection': 'keep-alive' +} +headersPayload = { + 'referrer': secureUrl, + 'origin': secureUrl, + 'Connection': 'keep-alive', + 'Content-Type': 'application/json' +} + +with requests.session() as s: + r = s.get(secureUrl) + print('get: ', r) + r = s.options(loginUrl, headers=headersPayloadOptions) + print('options: ', r) + r = s.post(loginUrl, json=jsonPayload, headers=headersPayload) + print('login post: ', r) + responseJson = r.json() + payload = { + 'source': 'login', + 'uid': responseJson['Profile']['Uid'], + 'token': responseJson['access_token'], + 'expires': responseJson['expires_in'], + 'profile[Identities]': responseJson['Profile']['Identities'], + 'profile[Password]': responseJson['Profile']['Password'], + 'profile[LastPasswordChangeDate]': responseJson['Profile']['LastPasswordChangeDate'], + 'profile[PasswordExpirationDate]': responseJson['Profile']['PasswordExpirationDate'], + 'profile[LastPasswordChangeToken]': responseJson['Profile']['LastPasswordChangeToken'], + 'profile[EmailVerified]': responseJson['Profile']['EmailVerified'], + 'profile[IsActive]': responseJson['Profile']['IsActive'], + 'profile[IsDeleted]': responseJson['Profile']['IsDeleted'], + 'profile[Uid]': responseJson['Profile']['Uid'], + 'profile[CustomFields][EBilling]': responseJson['Profile']['CustomFields']['EBilling'], + 'profile[CustomFields][AccountID]': responseJson['Profile']['CustomFields']['AccountID'], + 'profile[CustomFields][PreviousBillAmount]': responseJson['Profile']['CustomFields']['PreviousBillAmount'], + 'profile[CustomFields][PreviousDueDate]': responseJson['Profile']['CustomFields']['PreviousDueDate'], + 'profile[IsEmailSubscribed]': responseJson['Profile']['IsEmailSubscribed'], + 'profile[UserName]': responseJson['Profile']['UserName'], + 'profile[NoOfLogins]': responseJson['Profile']['NoOfLogins'], + 'profile[PhoneId]': responseJson['Profile']['PhoneId'], + 'profile[PhoneIdVerified]': responseJson['Profile']['PhoneIdVerified'], + 'profile[Roles]': responseJson['Profile']['Roles'], + 'profile[ExternalUserLoginId]': responseJson['Profile']['ExternalUserLoginId'], + 'profile[RegistrationProvider]': responseJson['Profile']['RegistrationProvider'], + 'profile[IsLoginLocked]': responseJson['Profile']['IsLoginLocked'], + 'profile[LoginLockedType]': responseJson['Profile']['LoginLockedType'], + 'profile[LastLoginLocation]': responseJson['Profile']['LastLoginLocation'], + 'profile[RegistrationSource]': responseJson['Profile']['RegistrationSource'], + 'profile[IsCustomUid]': responseJson['Profile']['IsCustomUid'], + 'profile[UnverifiedEmail]': responseJson['Profile']['UnverifiedEmail'], + 'profile[IsSecurePassword]': responseJson['Profile']['IsSecurePassword'], + 'profile[PrivacyPolicy]': responseJson['Profile']['PrivacyPolicy'], + 'profile[ExternalIds]': responseJson['Profile']['ExternalIds'], + 'profile[IsRequiredFieldsFilledOnce]': responseJson['Profile']['IsRequiredFieldsFilledOnce'], + 'profile[ConsentProfile]': responseJson['Profile']['ConsentProfile'], + 'profile[PIN]': responseJson['Profile']['PIN'], + 'profile[RegistrationData]': responseJson['Profile']['RegistrationData'], + 'profile[ID]': responseJson['Profile']['ID'], + 'profile[Provider]': responseJson['Profile']['Provider'], + 'profile[Prefix]': responseJson['Profile']['Prefix'], + 'profile[FirstName]': responseJson['Profile']['FirstName'], + 'profile[MiddleName]': responseJson['Profile']['MiddleName'], + 'profile[LastName]': responseJson['Profile']['LastName'], + 'profile[Suffix]': responseJson['Profile']['Suffix'], + 'profile[FullName]': responseJson['Profile']['FullName'], + 'profile[NickName]': responseJson['Profile']['NickName'], + 'profile[ProfileName]': responseJson['Profile']['ProfileName'], + 'profile[BirthDate]': responseJson['Profile']['BirthDate'], + 'profile[Gender]': responseJson['Profile']['Gender'], + 'profile[Website]': responseJson['Profile']['Website'], + 'profile[Email][0][Type]': responseJson['Profile']['Email'][0]['Type'], + 'profile[Email][0][Value]': responseJson['Profile']['Email'][0]['Value'], + 'profile[Country]': responseJson['Profile']['Country'], + 'profile[ThumbnailImageUrl]': responseJson['Profile']['ThumbnailImageUrl'], + 'profile[ImageUrl]': responseJson['Profile']['ImageUrl'], + 'profile[Favicon]': responseJson['Profile']['Favicon'], + 'profile[ProfileUrl]': responseJson['Profile']['ProfileUrl'], + 'profile[HomeTown]': responseJson['Profile']['HomeTown'], + 'profile[State]': responseJson['Profile']['State'], + 'profile[City]': responseJson['Profile']['City'], + 'profile[Industry]': responseJson['Profile']['Industry'], + 'profile[About]': responseJson['Profile']['About'], + 'profile[TimeZone]': responseJson['Profile']['TimeZone'], + 'profile[LocalLanguage]': responseJson['Profile']['LocalLanguage'], + 'profile[CoverPhoto]': responseJson['Profile']['CoverPhoto'], + 'profile[TagLine]': responseJson['Profile']['TagLine'], + 'profile[Language]': responseJson['Profile']['Language'], + 'profile[Verified]': responseJson['Profile']['Verified'], + 'profile[UpdatedTime]': responseJson['Profile']['UpdatedTime'], + 'profile[Positions]': responseJson['Profile']['Positions'], + 'profile[Educations]': responseJson['Profile']['Educations'], + 'profile[PhoneNumbers]': responseJson['Profile']['PhoneNumbers'], + 'profile[IMAccounts]': responseJson['Profile']['IMAccounts'], + 'profile[Addresses]': responseJson['Profile']['Addresses'], + 'profile[MainAddress]': responseJson['Profile']['MainAddress'], + 'profile[Created]': responseJson['Profile']['Created'], + 'profile[CreatedDate]': responseJson['Profile']['CreatedDate'], + 'profile[ModifiedDate]': responseJson['Profile']['ModifiedDate'], + 'profile[ProfileModifiedDate]': responseJson['Profile']['ProfileModifiedDate'], + 'profile[LocalCity]': responseJson['Profile']['LocalCity'], + 'profile[ProfileCity]': responseJson['Profile']['ProfileCity'], + 'profile[LocalCountry]': responseJson['Profile']['LocalCountry'], + 'profile[ProfileCountry]': responseJson['Profile']['ProfileCountry'], + 'profile[FirstLogin]': responseJson['Profile']['FirstLogin'], + 'profile[IsProtected]': responseJson['Profile']['IsProtected'], + 'profile[RelationshipStatus]': responseJson['Profile']['RelationshipStatus'], + 'profile[Quota]': responseJson['Profile']['Quota'], + 'profile[Quote]': responseJson['Profile']['Quote'], + 'profile[InterestedIn]': responseJson['Profile']['InterestedIn'], + 'profile[Interests]': responseJson['Profile']['Interests'], + 'profile[Religion]': responseJson['Profile']['Religion'], + 'profile[Political]': responseJson['Profile']['Political'], + 'profile[Sports]': responseJson['Profile']['Sports'], + 'profile[InspirationalPeople]': responseJson['Profile']['InspirationalPeople'], + 'profile[HttpsImageUrl]': responseJson['Profile']['HttpsImageUrl'], + 'profile[FollowersCount]': responseJson['Profile']['FollowersCount'], + 'profile[FriendsCount]': responseJson['Profile']['FriendsCount'], + 'profile[IsGeoEnabled]': responseJson['Profile']['IsGeoEnabled'], + 'profile[TotalStatusesCount]': responseJson['Profile']['TotalStatusesCount'], + 'profile[Associations]': responseJson['Profile']['Associations'], + 'profile[NumRecommenders]': responseJson['Profile']['NumRecommenders'], + 'profile[Honors]': responseJson['Profile']['Honors'], + 'profile[Awards]': responseJson['Profile']['Awards'], + 'profile[Skills]': responseJson['Profile']['Skills'], + 'profile[CurrentStatus]': responseJson['Profile']['CurrentStatus'], + 'profile[Certifications]': responseJson['Profile']['Certifications'], + 'profile[Courses]': responseJson['Profile']['Courses'], + 'profile[Volunteer]': responseJson['Profile']['Volunteer'], + 'profile[RecommendationsReceived]': responseJson['Profile']['RecommendationsReceived'], + 'profile[Languages]': responseJson['Profile']['Languages'], + 'profile[Projects]': responseJson['Profile']['Projects'], + 'profile[Games]': responseJson['Profile']['Games'], + 'profile[Family]': responseJson['Profile']['Family'], + 'profile[TeleVisionShow]': responseJson['Profile']['TeleVisionShow'], + 'profile[MutualFriends]': responseJson['Profile']['MutualFriends'], + 'profile[Movies]': responseJson['Profile']['Movies'], + 'profile[Books]': responseJson['Profile']['Books'], + 'profile[AgeRange]': responseJson['Profile']['AgeRange'], + 'profile[PublicRepository]': responseJson['Profile']['PublicRepository'], + 'profile[Hireable]': responseJson['Profile']['Hireable'], + 'profile[RepositoryUrl]': responseJson['Profile']['RepositoryUrl'], + 'profile[Age]': responseJson['Profile']['Age'], + 'profile[Patents]': responseJson['Profile']['Patents'], + 'profile[FavoriteThings]': responseJson['Profile']['FavoriteThings'], + 'profile[ProfessionalHeadline]': responseJson['Profile']['ProfessionalHeadline'], + 'profile[ProviderAccessCredential]': responseJson['Profile']['ProviderAccessCredential'], + 'profile[RelatedProfileViews]': responseJson['Profile']['RelatedProfileViews'], + 'profile[KloutScore]': responseJson['Profile']['KloutScore'], + 'profile[LRUserID]': responseJson['Profile']['LRUserID'], + 'profile[PlacesLived]': responseJson['Profile']['PlacesLived'], + 'profile[Publications]': responseJson['Profile']['Publications'], + 'profile[JobBookmarks]': responseJson['Profile']['JobBookmarks'], + 'profile[Suggestions]': responseJson['Profile']['Suggestions'], + 'profile[Badges]': responseJson['Profile']['Badges'], + 'profile[MemberUrlResources]': responseJson['Profile']['MemberUrlResources'], + 'profile[TotalPrivateRepository]': responseJson['Profile']['TotalPrivateRepository'], + 'profile[Currency]': responseJson['Profile']['Currency'], + 'profile[StarredUrl]': responseJson['Profile']['StarredUrl'], + 'profile[GistsUrl]': responseJson['Profile']['GistsUrl'], + 'profile[PublicGists]': responseJson['Profile']['PublicGists'], + 'profile[PrivateGists]': responseJson['Profile']['PrivateGists'], + 'profile[Subscription]': responseJson['Profile']['Subscription'], + 'profile[Company]': responseJson['Profile']['Company'], + 'profile[GravatarImageUrl]': responseJson['Profile']['GravatarImageUrl'], + 'profile[ProfileImageUrls]': responseJson['Profile']['ProfileImageUrls'], + 'profile[WebProfiles]': responseJson['Profile']['WebProfiles'], + 'profile[PinsCount]': responseJson['Profile']['PinsCount'], + 'profile[BoardsCount]': responseJson['Profile']['BoardsCount'], + 'profile[LikesCount]': responseJson['Profile']['LikesCount'], + 'profile[SignupDate]': responseJson['Profile']['SignupDate'], + 'profile[LastLoginDate]': responseJson['Profile']['LastLoginDate'], + 'profile[PreviousUids]': responseJson['Profile']['PreviousUids'] + } + print('Number of logins: ', responseJson['Profile']['NoOfLogins']) + r = s.post(authenticatorUrl, data=payload) + print('authenticator post: ', r) + r = s.get(secureUrl) + print('succesfully logged in') + + ssoMhlPayload = { + 'account': responseJson['Profile']['CustomFields']['AccountID'], + 'request_type': 'BILLING' + } + + r = s.post(ssoMhlUrl, data=ssoMhlPayload) + print('SSO: ', r) + responseSsoMhlJson = r.json() + + ssoTargetPayload = { + 'authssotoken': responseSsoMhlJson['token'] + } + r = s.post(ssoTargetUrl, data=ssoTargetPayload) + print('SSO Target: ', r) + + r = s.get(downloadDataUrl) + print('Download My Data: ', r) + + # TODO: Need to figure out how to read data diff --git a/Azure/secret.py b/Azure/secret.py index e69de29..54c4c60 100644 --- a/Azure/secret.py +++ b/Azure/secret.py @@ -0,0 +1,2 @@ +email = "alexander.rey1@gmail.com" +password = "7hW02#t03!" \ No newline at end of file diff --git a/EWR_D3D/.idea/EWR_D3D.iml b/EWR_D3D/.idea/EWR_D3D.iml index b9ba6dd..8f0515e 100644 --- a/EWR_D3D/.idea/EWR_D3D.iml +++ b/EWR_D3D/.idea/EWR_D3D.iml @@ -2,7 +2,7 @@ - + diff --git a/EWR_D3D/.idea/misc.xml b/EWR_D3D/.idea/misc.xml index e62796c..b189f70 100644 --- a/EWR_D3D/.idea/misc.xml +++ b/EWR_D3D/.idea/misc.xml @@ -1,4 +1,4 @@ - + \ No newline at end of file diff --git a/EWR_D3D/EWR_PlotHD_D3D.py b/EWR_D3D/EWR_PlotHD_D3D.py index 9b8e0dd..3ece528 100644 --- a/EWR_D3D/EWR_PlotHD_D3D.py +++ b/EWR_D3D/EWR_PlotHD_D3D.py @@ -15,8 +15,11 @@ import datetime as datetime from scipy.interpolate import LinearNDInterpolator, interp1d +import dfm_tools as dfmt from dfm_tools.get_nc import get_netdata, get_ncmodeldata, plot_netmapdata -from dfm_tools.get_nc_helpers import get_timesfromnc, get_ncfilelist, get_hisstationlist, get_ncvardimlist, get_timesfromnc +from dfm_tools.get_nc_helpers import get_ncvarproperties, get_stationid_fromstationlist, get_varnamefromattrs + + import cartopy.crs as ccrs import contextily as ctx from dfm_tools.regulargrid import scatter_to_regulargrid diff --git a/EWR_D3D/EWR_PlotHD_D3D_Profile.py b/EWR_D3D/EWR_PlotHD_D3D_Profile.py index cefd70a..10a6a64 100644 --- a/EWR_D3D/EWR_PlotHD_D3D_Profile.py +++ b/EWR_D3D/EWR_PlotHD_D3D_Profile.py @@ -64,7 +64,7 @@ dataPath = "FlowFM_map.nc" #%% Import using DFM functions #modelPlot = [35, 36, 37, 39, 40, 41, 42, 43, 44, 45] -modelPlot = [18, 35] +modelPlot = [18, 35, 41, 46, 47, 48, 49] dataIN_gp = [None] * (max(modelPlot)+1) dataOUT = [None] * (max(modelPlot)+1) @@ -181,14 +181,14 @@ ax = axes.flat modelPlot = [35, 37, 40, 41, 42, 43, 44] modelPlot = [40, 41, 42, 43, 44] modelPlot = [35, 37, 40, 45] -modelPlot = [18, 35] +modelPlot = [41, 47, 48, 49] for i in modelPlot: # S1 ax[0].set_title('Water Level') ax[0].plot(dataOUT[i].index, dataOUT[i]['s1'], linewidth=3, label=runLog['Run Title'][i]) - ax[0].set_ylim([330, 345]) + ax[0].set_ylim([325, 357]) # ax[0].set_xlim([0, 10000]) ax[0].set_xlim([0, 90000]) @@ -198,7 +198,7 @@ for i in modelPlot: # Shear ax[1].set_title('Shear Stress') ax[1].plot(dataOUT[i].index, dataOUT[i]['taus'], linewidth=3, label=runLog['Run Title'][i]) - ax[1].set_ylim([0, 2]) + ax[1].set_ylim([0, 0.25]) ax[1].legend(loc='upper right') ax[1].set_ylabel('Shear Stress [N/m^2]') @@ -233,8 +233,8 @@ plt.show() # fig.savefig('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/12_WAQ/FullDomainFigures/DomainHDCompare.png', # bbox_inches='tight', dpi=200) -fig.savefig('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/07_PostProcessing/HD_Profile_Compare.png', - bbox_inches='tight', dpi=200) +# fig.savefig('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/07_PostProcessing/HD_Profile_Compare.png', +# bbox_inches='tight', dpi=200) #%% Plot select data along centerline fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 6), sharex=True) diff --git a/EWR_D3D/EWR_Plotting.py b/EWR_D3D/EWR_Plotting.py index 83a256b..b287f82 100644 --- a/EWR_D3D/EWR_Plotting.py +++ b/EWR_D3D/EWR_Plotting.py @@ -594,8 +594,8 @@ ax[1].set_xlim([0, 100000]) #DetC in sediment ax[2].set_title('DetC (Fast) in Sediment') ax[2].plot(dataOUT.index, dataOUT['DetCS1'] - / (dataOUT['IM1S1'] + dataOUT['IM2S1'] + dataOUT['IM3S1'] + dataOUT['DetCS1'] + dataOUT['OOCS1']) - , linewidth=3, label='DetCS1') + / (dataOUT['IM1S1'] + dataOUT['IM2S1'] + dataOUT['IM3S1'] + dataOUT['DetCS1'] + dataOUT['OOCS1']), + linewidth=3, label='DetCS1') ax[2].legend() ax[2].set_ylabel('DetCS1 [g/g]') diff --git a/EWR_D3D/WSP_DB.py b/EWR_D3D/WSP_DB.py index e69de29..5f594fc 100644 --- a/EWR_D3D/WSP_DB.py +++ b/EWR_D3D/WSP_DB.py @@ -0,0 +1,90 @@ +#%% Script to read and work with Access database + +import pyodbc +import pandas as pd +import geopandas as gpd + +#%% Establish a connection to the Access database +conn_str = r'DRIVER={Microsoft Access Driver (*.mdb, *.accdb)};DBQ=C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/Wood English-Wabigoon Export_alldata_2017-current.accdb;' +conn = pyodbc.connect(conn_str) + +# Get the list of tables +tables = conn.cursor().tables() + +# Save a list of the table names +table_names = [table.table_name for table in tables] + + +#%% Print rows from a table +cursor = conn.cursor() +cursor.execute('SELECT * FROM [General_Query]') +for row in cursor.fetchall(): + print(row) + + +#%% Print the column names +cursor = conn.cursor() +cursor.execute('SELECT * FROM [General_Query]') +columns = [column[0] for column in cursor.description] +print(columns) + +#%% Save the table to a pandas dataframe +df = pd.read_sql_query('SELECT * FROM [General_Query]', conn) + +#%% Save the table to a pandas dataframe +df_Samples = pd.read_sql_query('SELECT * FROM [EW_Export_Samples]', conn) + +#%% Save the table to a pandas dataframe +df_Results = pd.read_sql_query('SELECT * FROM [EW_Export_Results]', conn) + +#%% Save the table to a pandas dataframe +df_Locations = pd.read_sql_query('SELECT * FROM [EW_Export_Locations]', conn) + +#%% Close the connection +conn.close() + + +#%% Select a subset of the data where the media is Surface Water ('SW') and the param_name is 'Mercury' and the fraction is 'T' +df_SW_Hg = df.loc[(df['media'] == 'SW') & (df['param_name'] == 'Mercury') & (df['fraction'] == 'T'), ['Long_coord', 'Lat_coord', 'ppm_result', 'ppm_uom', 'field_sample_date', 'param_name', 'tech_task_name']] + +# Convert to geo dataframe using Lat and Long +gdf_SW_Hg = gpd.GeoDataFrame(df_SW_Hg, geometry=gpd.points_from_xy(df_SW_Hg.Long_coord, df_SW_Hg.Lat_coord, crs="EPSG:4326")) + +# Convert result to numeric +gdf_SW_Hg['ppm_result'] = pd.to_numeric(gdf_SW_Hg['ppm_result']) + +# Where ppm_uom is 'MG/L', convert to NG/L +gdf_SED_Hg.loc[gdf_SED_Hg['ppm_uom'] == 'MG/L', 'ppm_result'] = gdf_SED_Hg['ppm_result'] * 1000000 + +# Drop where no result or location +gdf_SW_Hg = gdf_SW_Hg.dropna(subset=['ppm_result', 'Long_coord', 'Lat_coord']) + +# Convert datetime to string for shapefile +gdf_SW_Hg['field_sample_date'] = gdf_SW_Hg['field_sample_date'].dt.strftime('%Y-%m') + +# Save as a shapefile +gdf_SW_Hg.to_file('C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/EW_Export_Samples_SW_Hg_RevB.shp') + + + +#%% Select a subset of the data where: +# - media is Sediment ('SED') +# - the param_name is 'Mercury' +# - The top depth is 0 + +df_SED_Hg = df.loc[(df['media'] == 'SED') & (df['param_name'] == 'Mercury') & (df['top_depth'] == 0), ['Long_coord', 'Lat_coord', 'ppm_result', 'ppm_uom']] + +# Convert to geo dataframe using Lat and Long +gdf_SED_Hg = gpd.GeoDataFrame(df_SED_Hg, geometry=gpd.points_from_xy(df_SED_Hg.Long_coord, df_SED_Hg.Lat_coord, crs="EPSG:4326")) + +# Convert result to numeric +gdf_SED_Hg['ppm_result'] = pd.to_numeric(gdf_SED_Hg['ppm_result']) + +# Where ppm_uom is 'MG/KG', convert to NG/G +gdf_SED_Hg.loc[gdf_SED_Hg['ppm_uom'] == 'MG/KG', 'ppm_result'] = gdf_SED_Hg['ppm_result'] * 1000 + +# Drop rows with no locatuion +gdf_SED_Hg = gdf_SED_Hg.dropna(subset=['Long_coord', 'Lat_coord']) + +# Save as a shapefile +df_SW_Hg.to_file('C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/EW_Export_Samples_SW_Hg.shp') diff --git a/EWR_Data/EWR_Hg_Interp.py b/EWR_Data/EWR_Hg_Interp.py index 5697806..915ca0e 100644 --- a/EWR_Data/EWR_Hg_Interp.py +++ b/EWR_Data/EWR_Hg_Interp.py @@ -26,6 +26,8 @@ dataIN_gp = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Longitude, d # Convert to UTM dataIN_gp.to_crs(crs='EPSG:32615', inplace=True) + + #%% Read in centerline shapefile river_centerline = gp.read_file('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/16_Waterline/Centerline_for_Modelling_UTMZ15.shp') diff --git a/EWR_Data/EWR_Hg_Interp_RevB.py b/EWR_Data/EWR_Hg_Interp_RevB.py index b1dba11..11faf1d 100644 --- a/EWR_Data/EWR_Hg_Interp_RevB.py +++ b/EWR_Data/EWR_Hg_Interp_RevB.py @@ -42,9 +42,22 @@ river_centerline_merge_gpd2.iloc[0, 2] = 0 #%% Read in combined dataset -obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/output/combined dataset-SGJ.csv") +# obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/output/combined dataset-SGJ.csv") +# obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/combined dataset.csv") +# obs_IN = pd.read_csv("C:/Users/arey/Downloads/ExportDownloads_EwrrpDataExport_rharris_2023-05-12_104911/EwrrpDataExport_rharris_2023-05-12_104911.txt", sep='\t', encoding='UTF-16 LE') +df = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Results.txt", sep = '\t', encoding='UTF-16 LE') +locations = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Sites.txt", sep = '\t', encoding='UTF-16 LE') +samples = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Samples.txt", sep = '\t', encoding='UTF-16 LE') + +#%% Merge locations and results +obs_IN = pd.merge(left=df, right=locations, left_on= 'Samplesiteid', right_on= 'Id', how = 'left') # 1096 locations w/Hg +obs_IN = pd.merge(left=obs_IN, right=samples, left_on= 'Sampleid', right_on= 'Id', how = 'left') # 1096 locations w/Hg + + # Add Datetime -obs_IN['DateTime'] = pd.to_datetime(obs_IN.Sampledate_x) +obs_IN['Sampledate'] = pd.to_datetime(obs_IN.Sampledate) + + #%% Process combined datsaset # Convert to geodataframe @@ -54,11 +67,11 @@ obs_gdf = gp.GeoDataFrame(obs_IN, geometry=gp.points_from_xy( obs_gdf = gp.sjoin_nearest(river_centerline_merge_gpd2, obs_gdf, how='right', max_distance=250).reset_index() # Sediment Hg GeoDataFrame where media is Sediment -obsMask = (((obs_gdf.Media_x == 'SED') | (obs_gdf.Media_x == 'SOIL')) & +obsMask = (((obs_gdf.Media == 'SED') | (obs_gdf.Media == 'SOIL')) & ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & - (obs_gdf.DateTime > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna() & - ((obs_gdf.TopDepth_x == 0) | (obs_gdf.TopDepth_x.isna()))) + (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna() & + ((obs_gdf.TopDepth == 0) | (obs_gdf.TopDepth.isna()))) # obsMask = (((obs_gdf.Media_x == 'SED') | (obs_gdf.Media_x == 'SOIL')) & # ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') # | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & @@ -74,6 +87,42 @@ obs_HgSed.loc[obs_HgSed.Unit == 'ug/g', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Samplevalue'] * 0.001 obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Samplevalue'] * 0.001 + + +# Save subset of columns as a shapefile +obs_HgSed.loc[:, ['Sample_NG/G', 'Sampledate_x', 'RiverKM', 'Latitude', 'Longitude', 'geometry', + 'Surveyname','Sitecode','Parameter']].to_file( + 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgSed.shp') + + +#%% Extract Hg data from surface water samples + +# Water Hg GeoDataFrame where media is SurfaceWater (SW) +obsMask = (((obs_gdf.Media == 'SW')) & + ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)-Total') + | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & + (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna()) + +obs_HgWater = obs_gdf.loc[obsMask, :].copy() +# Adjust Units +obs_HgWater.loc[obs_HgWater.Unit == 'NG/L', 'Samplevalue'] = obs_HgWater.Samplevalue +obs_HgWater.loc[obs_HgWater.Unit == 'ng/L', 'Samplevalue'] = obs_HgWater.Samplevalue +obs_HgWater.loc[obs_HgWater.Unit == 'UG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 +obs_HgWater.loc[obs_HgWater.Unit == 'ug/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 +obs_HgWater.loc[obs_HgWater.Unit == 'MG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 +obs_HgWater.loc[obs_HgWater.Unit == 'mg/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 + +# Drop where no result +obs_HgWater = obs_HgWater.dropna(subset=['Samplevalue']) + +# Save subset of columns as a shapefile +obs_HgWater.loc[:, ['Samplevalue', 'Sampledate', 'RiverKM', 'geometry', + 'Contributor', 'Sitecode', 'Parameter']].to_file( + 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgWater_RevB.shp') + + + + #%% River Profile Plotting fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 6)) diff --git a/EWR_Data/EWR_Hg_Syn.py b/EWR_Data/EWR_Hg_Syn.py index 11faf1d..37d0754 100644 --- a/EWR_Data/EWR_Hg_Syn.py +++ b/EWR_Data/EWR_Hg_Syn.py @@ -40,198 +40,59 @@ river_centerline_merge_gpd2['RiverKM'] = river_centerline_merge_gpd2['DistanceFr river_centerline_merge_gpd2.iloc[0, 1] = 0 river_centerline_merge_gpd2.iloc[0, 2] = 0 +#%% Read in water depth to get wet and dry points -#%% Read in combined dataset -# obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/output/combined dataset-SGJ.csv") -# obs_IN = pd.read_csv("//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/05_Analyses/02_Data Analysis/combined dataset.csv") -# obs_IN = pd.read_csv("C:/Users/arey/Downloads/ExportDownloads_EwrrpDataExport_rharris_2023-05-12_104911/EwrrpDataExport_rharris_2023-05-12_104911.txt", sep='\t', encoding='UTF-16 LE') -df = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Results.txt", sep = '\t', encoding='UTF-16 LE') -locations = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Sites.txt", sep = '\t', encoding='UTF-16 LE') -samples = pd.read_csv(r"P:\12828.101 English Wabigoon River\03_Data\02_Physical\30_Site Sampling Database\EwrrpDataExport_rharris_2022-08-23_102403_Samples.txt", sep = '\t', encoding='UTF-16 LE') +depth_xyz = pd.read_csv('//srv-ott3/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/04_Bathymetry/Water depth at pressure points - mesh2d_nFaces_ mean.csv', + names=['x', 'y', 'z', 'depth'], header=1) -#%% Merge locations and results -obs_IN = pd.merge(left=df, right=locations, left_on= 'Samplesiteid', right_on= 'Id', how = 'left') # 1096 locations w/Hg -obs_IN = pd.merge(left=obs_IN, right=samples, left_on= 'Sampleid', right_on= 'Id', how = 'left') # 1096 locations w/Hg +# Make geodataframe +depth_xyz_gp = gp.GeoDataFrame(depth_xyz, geometry=gp.points_from_xy(depth_xyz.x, depth_xyz.y), crs="EPSG:32615") + +# Merge with centerline to get river km +depth_xyz_gp = gp.sjoin_nearest(river_centerline_merge_gpd2, depth_xyz_gp, how='right') + +#%% Create synthetic Hg data following Reed's equation +# These need to be multiplied by the total mass of dry matter in each cell to get the total Hg/m2 in each layer +CstartHg = 10 +CfinalHg = 1 +kHg = 0.08 + +depth_xyz_gp['Hg'] = (CstartHg - CfinalHg) * np.exp(-kHg * depth_xyz_gp['RiverKM']/1000) + CfinalHg + +CstartMeHg = 0.05 +CfinalMeHg = 0.005 +kMeHg = 0.08 + +depth_xyz_gp['MeHg'] = (CstartMeHg - CfinalMeHg) * np.exp(-kMeHg * depth_xyz_gp['RiverKM']/1000) + CfinalMeHg -# Add Datetime -obs_IN['Sampledate'] = pd.to_datetime(obs_IN.Sampledate) +# Set points where depth is zero (no water) to 0 +depth_xyz_gp.loc[depth_xyz_gp['depth'] == 0, 'Hg'] = 0 +depth_xyz_gp.loc[depth_xyz_gp['depth'] == 0, 'MeHg'] = 0 - - -#%% Process combined datsaset -# Convert to geodataframe -obs_gdf = gp.GeoDataFrame(obs_IN, geometry=gp.points_from_xy( - obs_IN.loc[:, 'Longitude'], obs_IN.loc[:, 'Latitude'], crs="EPSG:4326")).to_crs(crs="EPSG:32615") -# Add centerline and reset index -obs_gdf = gp.sjoin_nearest(river_centerline_merge_gpd2, obs_gdf, how='right', max_distance=250).reset_index() - -# Sediment Hg GeoDataFrame where media is Sediment -obsMask = (((obs_gdf.Media == 'SED') | (obs_gdf.Media == 'SOIL')) & - ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') - | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & - (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna() & - ((obs_gdf.TopDepth == 0) | (obs_gdf.TopDepth.isna()))) -# obsMask = (((obs_gdf.Media_x == 'SED') | (obs_gdf.Media_x == 'SOIL')) & -# ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg) Total') -# | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & -# (obs_gdf.DateTime > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna()) - -obs_HgSed = obs_gdf.loc[obsMask, :].copy() - -# Adjust Units -obs_HgSed.loc[obs_HgSed.Unit == 'NG/G', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'NG/G', 'Samplevalue'] -obs_HgSed.loc[obs_HgSed.Unit == 'ng/g', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'ng/g', 'Samplevalue'] -obs_HgSed.loc[obs_HgSed.Unit == 'UG/G', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'UG/G', 'Samplevalue'] * 1000 -obs_HgSed.loc[obs_HgSed.Unit == 'ug/g', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'ug/g', 'Samplevalue'] * 1000 -obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'MG/KG', 'Samplevalue'] * 0.001 -obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Sample_NG/G'] = obs_HgSed.loc[obs_HgSed.Unit == 'mg/kg', 'Samplevalue'] * 0.001 - - - -# Save subset of columns as a shapefile -obs_HgSed.loc[:, ['Sample_NG/G', 'Sampledate_x', 'RiverKM', 'Latitude', 'Longitude', 'geometry', - 'Surveyname','Sitecode','Parameter']].to_file( - 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgSed.shp') - - -#%% Extract Hg data from surface water samples - -# Water Hg GeoDataFrame where media is SurfaceWater (SW) -obsMask = (((obs_gdf.Media == 'SW')) & - ((obs_gdf.Parameter == 'Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)-Total') - | (obs_gdf.Parameter == 'Total Mercury') | (obs_gdf.Parameter == 'Mercury (Hg)')) & - (obs_gdf.Sampledate > datetime.datetime(2000, 1, 1)) & obs_gdf.RiverKM.notna()) - -obs_HgWater = obs_gdf.loc[obsMask, :].copy() -# Adjust Units -obs_HgWater.loc[obs_HgWater.Unit == 'NG/L', 'Samplevalue'] = obs_HgWater.Samplevalue -obs_HgWater.loc[obs_HgWater.Unit == 'ng/L', 'Samplevalue'] = obs_HgWater.Samplevalue -obs_HgWater.loc[obs_HgWater.Unit == 'UG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 -obs_HgWater.loc[obs_HgWater.Unit == 'ug/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000 -obs_HgWater.loc[obs_HgWater.Unit == 'MG/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 -obs_HgWater.loc[obs_HgWater.Unit == 'mg/L', 'Samplevalue'] = obs_HgWater.Samplevalue * 1000000 - -# Drop where no result -obs_HgWater = obs_HgWater.dropna(subset=['Samplevalue']) - -# Save subset of columns as a shapefile -obs_HgWater.loc[:, ['Samplevalue', 'Sampledate', 'RiverKM', 'geometry', - 'Contributor', 'Sitecode', 'Parameter']].to_file( - 'C:/Users/arey/files/Projects/Grassy Narrows/Wood English-Wabigoon Export_alldata_2017-current/BairdDB_HgWater_RevB.shp') - - - - -#%% River Profile Plotting -fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 6)) - -axes.scatter(obs_HgSed.RiverKM, obs_HgSed.Samplevalue) -axes.scatter(obs_HgSed.RiverKM, obs_HgSed.Samplevalue) -axes.set_ylim(0, 20000) -axes.set_xlim(0, 100000) -axes.set_xlabel('Distance Along River [m]') -axes.set_ylabel('Surface Sediment Hg [ng/g]') - -plt.show() - -fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Hg_SurfaceProfile.png', - bbox_inches='tight', dpi=300) - -#%% Plot Vertical Profiles -fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(7, 7)) -obs_HgSed.plot.scatter('RiverKM', 'TopDepth_x', 12, 'Sample_NG/G', - ax=axes, vmax=20000, cmap='viridis') - -axes.invert_yaxis() -axes.set_xlabel('Distance Along River [m]') -axes.set_ylabel('Sample Depth [cm]') -axes.set_xlim([0, 100000]) - -# Get colorbar axis to set a legend -cax = fig.get_axes()[1] -#and we can modify it, i.e.: -cax.set_ylabel('Sediment Hg (ng/g)') - -plt.show() - -fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Hg_Vertical.png', - bbox_inches='tight', dpi=300) - -#%% Hg Map Plotting +#%% Hg Test plot mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 7)) fig.patch.set_facecolor('white') -# Full system limits -# axes.set_xlim(350000, 520000) -# axes.set_ylim(5510000, 5580000) - -# # Clay Lake Limits -# axes.set_xlim(466418, 515846) -# axes.set_ylim(5513437, 5545362) - -# Dryden Limits -axes.set_xlim(497697, 515846) -axes.set_ylim(5513437, 5525166) - - # Add Basemap -ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615') +# ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615') -# Plot Hg At surface -obs_HgSed.plot(ax=axes, column='Samplevalue', legend=True, - vmin=0, vmax=5000, markersize=2, cmap='viridis', - legend_kwds={'label': 'Surface Sediment Hg (ng/g)'}) +# Plot river distance using geopandas +# depth_xyz_gp.plot(ax=axes, column='RiverKM', cmap='viridis') +depth_xyz_gp.loc[depth_xyz_gp['depth'] != 0, :].plot(ax=axes, column='Hg', cmap='viridis') plt.show() -fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceHgMap.png', - bbox_inches='tight', dpi=300) +#%% Save to csv +np.savetxt('//srv-ott3/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/99_Setups/SynHg_ug_g.xyz', + np.array([depth_xyz_gp.geometry.x.values, depth_xyz_gp.geometry.y.values, depth_xyz_gp['Hg'].values]).T, + delimiter=" ") - -#%% Interpolate Hg -bathy_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Bed_Level.xyz', - names=['x', 'y', 'z'], header=0, delim_whitespace=True) - -rough_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Roughness.xyz', - names=['x', 'y', 'z'], header=0, delim_whitespace=True) - -HgInterp = griddata(np.array([dataOUTmax.geometry.x, dataOUTmax.geometry.y]).T, - np.array(dataOUTmax.Samplevalue).T, np.array([bathy_xyz.x, bathy_xyz.y]).T, - method='nearest') - -HgInterpOut = np.vstack((bathy_xyz.x, bathy_xyz.y, HgInterp)) -HgInterpOut = np.transpose(HgInterpOut) - -# Set Wetlands to zero -# Interpolate Roughness -RoughInterp = sp.interpolate.griddata(np.transpose(np.array([rough_xyz.x, rough_xyz.y])), rough_xyz.z, np.transpose(np.array([bathy_xyz.x, bathy_xyz.y])), - method='nearest') - -HgInterpOut[RoughInterp==0.05, 2] = 0 - -np.savetxt('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceInterFiltDB.xyz', - HgInterpOut, delimiter=" ") - - -#%% Hg Interp Plotting -mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' - -fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 7)) -fig.patch.set_facecolor('white') - -# Clay Lake Limits -axes.set_xlim(466418, 515846) -axes.set_ylim(5513437, 5545362) -# Add Basemap -ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615') -# Plot Hg -axes.scatter(HgInterpOut[:, 0], HgInterpOut[:, 1], 5, HgInterpOut[:, 2], - vmin=0, vmax=5000) - -plt.show() +np.savetxt('//srv-ott3/Projects/12828.101 English Wabigoon River/06_Models/00_Delft3D/99_Setups/SynMeHg_ug_g.xyz', + np.array([depth_xyz_gp.geometry.x.values, depth_xyz_gp.geometry.y.values, depth_xyz_gp['MeHg'].values]).T, + delimiter=" ") diff --git a/EWR_Data/createSeg.py b/EWR_Data/createSeg.py index e69de29..6ab54da 100644 --- a/EWR_Data/createSeg.py +++ b/EWR_Data/createSeg.py @@ -0,0 +1,33 @@ +#%% Create a Delft3D WAQ segment input file +# Alexander Rey, 2023-06-27 + +# These files can be used to define a value for every cell at every time step +# The files are binary and have a specific format +# The format is described in the Delft3D WAQ manual +# Each lauer is a different segmnent +# Typically, segments increase sequentially from 1 to N, then the second layer will be N+1 to 2N, etc. +# The number of layers is defined in the GUI +# Use QUICKPLOT to extract segment numbers in layer 1, then adjust from there +# Note: it will go to zero at the end of the time series, so ensure that the last time step is covered + +#%% Import libraries +import struct + +NOSEG = 5 # Number of segments +NOTIME = 10 # Number of time steps + +VALUE = [[0.0] * NOSEG for _ in range(NOTIME)] # Initialize VALUE as a 2D matrix + +for i in range(NOTIME): + for j in range(NOSEG): + # Fill the matrix with the desired values. + VALUE[i][j] = i + j + +with open('output_seg10.bin', 'wb') as file: + for i in range(NOTIME): + # TIMER is the number of seconds from the reference date (set in the GUI) + TIMER = (i - 1) * 31536000 + file.write(struct.pack('i', TIMER)) # Write TIMER as a 4-byte integer + + for j in range(NOSEG): + file.write(struct.pack('f', VALUE[i][j])) # Write VALUE[i][j] as a 4-byte float diff --git a/EWR_Data/output_seg.seg b/EWR_Data/output_seg.seg new file mode 100644 index 0000000..177e962 --- /dev/null +++ b/EWR_Data/output_seg.seg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/EWR_Data/output_seg10.bin b/EWR_Data/output_seg10.bin new file mode 100644 index 0000000000000000000000000000000000000000..90a85c074386ce9e60dbf855696afe492f239145 GIT binary patch literal 240 zcmZoTBlnL13>xf#lmid}aRU&80J_)$hX&(^j6ewogt`Sl_5lZm^h3-TVhj%)8fGnK zg_s9b2NVP9cVsx%!p^{eCdMG(*ua>^2{sR=4k*UJpy0@0>dwu8F2-Qs*pQ*Y3pNj? U4k!k6has z%8Cog*hCO{)s`&oO}?fxSttzv literal 0 HcmV?d00001 diff --git a/LSU_D3D/PenobscotExample.py b/LSU_D3D/PenobscotExample.py new file mode 100644 index 0000000..b8bb90f --- /dev/null +++ b/LSU_D3D/PenobscotExample.py @@ -0,0 +1,67 @@ +#%% Example Penobscot Plotting Script +# Alexander Rey, May 18, 2023 + +#%% Import modules +import xarray as xr +import matplotlib.pyplot as plt +import contextily as ctx +import numpy as np +import dfm_tools as dfmt +import pandas as pd +import pathlib as pl + +#%% Read in Model Log +pth = pl.Path("N:/Penobscot/Modelling/12345.678.W.AJMR.Rev0_PenobscotModelLog.xlsx") + +runLog = pd.read_excel(pth.as_posix(), "ModelLog") + +dataPath = "FlowFM_map.nc" + + +#%% Read in Delft3D Data +i = 3 + +# Get hist file path +file_nc_hist = pl.Path(runLog['Run Location'][i], + 'FlowFM', 'output', 'FlowFM_his.nc') +file_nc_hist = "N:/Penobscot/Modelling/Runs/04_SalinityRun.dsproj_data/FlowFM/output/FlowFM_his.nc" +# Open hist file dataset +data_xr = xr.open_mfdataset(file_nc_hist, preprocess=dfmt.preprocess_hisnc) + +# Get stations +stations_pd = data_xr['stations'].to_dataframe() + +#%% Extract and plot water level + +# Extract water level +waterlevel_xr = data_xr['waterlevel'].sel(stations=stations_pd.index[1]) + +# Plot water level +fig, ax = plt.subplots(figsize=(12, 6)) +waterlevel_xr.plot(ax=ax, label='Water Level') + +plt.show() + +#%% Extract and plot water level map +file_nc_map = "N:/Penobscot/Modelling/Runs/04_SalinityRun.dsproj_data/FlowFM/output/FlowFM_map.nc" + +# Read in map file +d3dfm_DataArray = dfmt.open_partitioned_dataset(file_nc_map) + +# Get Var info +vars_pd = dfmt.get_ncvarproperties(d3dfm_DataArray) + +# Get map of water level +waterlevel_map_xr = d3dfm_DataArray.mesh2d_s1[3, :].compute() + + +#%% Plot water level map using delft3d map plotting function and contextily +fig, ax = plt.subplots(figsize=(12, 6)) + +pc = waterlevel_map_xr.ugrid.plot(ax=ax, edgecolors='face', cmap='jet') + +# Define crs +crs = 'EPSG:32619' +ctx.add_basemap(ax=ax, source=ctx.providers.Esri.WorldImagery, crs=crs, attribution=False) + +plt.show() diff --git a/Mustique/WestCoastDataTemplate_V5.py b/Mustique/WestCoastDataTemplate_V5.py index 271cfb8..f88dae1 100644 --- a/Mustique/WestCoastDataTemplate_V5.py +++ b/Mustique/WestCoastDataTemplate_V5.py @@ -10,6 +10,8 @@ import math from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar, AnchoredDirectionArrows import matplotlib.pyplot as plt import matplotlib.font_manager as fm +import matplotlib.dates as mdates + import matplotlib as mpl import cartopy.crs as ccrs @@ -81,7 +83,8 @@ importPaths = ['//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Resear None, '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/20230212 WQ Sampling/Old Queens Fort', '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/20230212 WQ Sampling/Crane', - '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/20230212 WQ Sampling/Caribee'] + '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/20230212 WQ Sampling/Caribee', + '//srv-ott3.baird.com/Projects/13711.101 Caribee (Indigo) Hotel, Barbados/03_Data/03_Field/20220613_EnvData'] siteNames = ['Great House', 'Greensleeves', @@ -121,6 +124,7 @@ siteNames = ['Great House', None, 'Old Queens Fort', 'Crane', + 'Caribee', 'Caribee'] timeLabels= ['Before Construction', @@ -161,7 +165,8 @@ timeLabels= ['Before Construction', None, 'February', 'February', - 'February'] + 'February', + 'June 2022'] wave_bts_file = [ 'T:/Alexander/WestCoast/Barbados Nowcast 2021-09-15 to 2021-11-15/spawnee_mid_27m.bts', @@ -202,33 +207,69 @@ wave_bts_file = [ None, None, None, + None, None] # %% Read in site shapefile siteShp = gp.read_file('//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/SitePolygons_Crane_Caribee.shp') siteShp.geometry = siteShp.geometry.to_crs("EPSG:32621") -for s in [34, 36, 37, 38]: #range(15, 27): +# Switch to save figures/ sheets +figureSave = True + +for s in range(22, 40): #range(15, 27): + ## Define master import path importPath = importPaths[s] siteName = siteNames[s] timeLabel = timeLabels[s] importFiles = os.listdir(importPath) + # Skip if no importPath + if importPath is None: + print('Skipping site ' + str(s + 1)) + continue + + # Monitor progress + print('Processing site ' + str(s+1) + ' of ' + str(len(importPaths)) + ' (' + siteName + '_' + timeLabel + ')') + # Initialize import variables OBS_File = None + GPS_File = None + RBR_File = None + JFE_File = None for i in range(0, len(importFiles)): - if '.csv' in importFiles[i] and 'Summary' not in importFiles[i] and 'export' not in importFiles[i]: + if '_A.csv' in importFiles[i] and 'Summary' not in importFiles[i]: + JFE_File = importFiles[i] + elif '.csv' in importFiles[i] and 'Summary' not in importFiles[i] and 'export' not in importFiles[i]: OBS_File = importFiles[i] elif '.csv' in importFiles[i] and 'export' in importFiles[i]: GPS_File = importFiles[i] GPS_Type = 'CSV' + elif '.xlsx' in importFiles[i] and 'Summary' not in importFiles[i]: + RBR_File = importFiles[i] elif '.xls' in importFiles[i] and 'Summary' not in importFiles[i]: GPS_File = importFiles[i] GPS_Type = 'xls' + #%% RBR Import Data + if RBR_File is not None: + RBR_Obs = pd.read_excel(importPath + '/' + RBR_File, + sheet_name='Data', skiprows=0, header=1) + RBR_Obs['DateTime'] = pd.to_datetime(RBR_Obs['Time']) + RBR_Obs.set_index('DateTime', inplace=True) + + #%% JFE Import Data + if JFE_File is not None: + JFE_Obs = pd.read_csv(importPath + '/' + JFE_File, skiprows=30) + JFE_Obs['DateTime'] = pd.to_datetime(JFE_Obs['Date']) + # Add seconds to time based on mod 60 of index + JFE_Obs['DateTime'] = JFE_Obs['DateTime'] + pd.to_timedelta(JFE_Obs.index % 60, unit='s') + JFE_Obs.set_index('DateTime', inplace=True) + + #%% Obs Import Data if OBS_File is not None: Obs_dat = pd.read_csv(importPath + '/' + OBS_File, skiprows=0, header=0) @@ -245,6 +286,26 @@ for s in [34, 36, 37, 38]: #range(15, 27): # Set Index Obs_dat.set_index('DateTime', inplace=True) + else: + # Merge RBR and JFS obs into one dataframe + if JFE_File is not None: + Obs_dat = pd.merge(RBR_Obs, JFE_Obs, how='outer', on='DateTime') + else: + Obs_dat = RBR_Obs + + # Sort index by DateTime + Obs_dat.sort_index(inplace=True) + + # Rename columns to match Obs_dat + Obs_dat = Obs_dat.rename(columns={'Pressure ': 'CH0:Pressure(dbar)', 'Temperature ': 'CH1:Temperature(degC)', + 'Chl-a[ug/l]': 'CH2:Chlorophyll_a(ppb)', 'Conductivity ': 'CH3:Conductivity(mS/cm)', + 'Temp.[deg C]': 'CH4:Temperature(degC)', 'Turb. -M[FTU]': 'CH6:Turbidity(NTU)', + 'Dissolved O₂ saturation ': 'CH8:Oxygen_Saturation(%)'}) + + # Set Time Zone for Sensor. Sometimes it is set to UTC, sometimes it is set to EST + Obs_dat.index = Obs_dat.index.tz_localize('UTC') + + #%% GPS Import Data if GPS_File is not None: if GPS_Type == 'CSV': @@ -267,22 +328,29 @@ for s in [34, 36, 37, 38]: #range(15, 27): GPS_gdf['DateTime'] = pd.to_datetime(GPS_gdf['time']).dt.tz_localize('UTC') - # Set Datetime as index - GPS_gdf.set_index('DateTime', inplace=True) + else: + # Synthesize GPS data for great house + GPS_gdf = gp.read_file('C:/Users/arey/files/Projects/West Coast/GreatHousePath_R3.shp') + GPS_gdf['DateTime'] = pd.date_range(pd.to_datetime(RBR_Obs['Time'].iloc[0]), + pd.to_datetime(RBR_Obs['Time'].iloc[-1]), + periods=len(GPS_gdf)).values - # Sort by time - GPS_gdf.sort_index(inplace=True) + # Set Datetime as index + GPS_gdf.set_index('DateTime', inplace=True) - # Convert to UTM - GPS_gdf.geometry = GPS_gdf.geometry.to_crs("EPSG:32621") + # Sort by time + GPS_gdf.sort_index(inplace=True) - #%% Merge GPS to RBR - # Process RBR into datetime - if OBS_File is not None: - # Merge with GPS as dataframe - Obs_geo = pd.merge_asof(Obs_dat, GPS_gdf, - left_index=True, right_index=True, direction='nearest', tolerance=pd.Timedelta('300s')) - Obs_geo = gp.GeoDataFrame(Obs_geo, geometry=Obs_geo.geometry, crs="EPSG:32621") + # Convert to UTM + GPS_gdf.geometry = GPS_gdf.geometry.to_crs("EPSG:32621") + + + #%% Merge Obs to RBR + # Process Obs into datetime + # Merge with GPS as dataframe + Obs_geo = pd.merge_asof(Obs_dat, GPS_gdf, + left_index=True, right_index=True, direction='nearest', tolerance=pd.Timedelta('300s')) + Obs_geo = gp.GeoDataFrame(Obs_geo, geometry=Obs_geo.geometry, crs="EPSG:32621") #%% Find and setup key plotting details # Shaded Water @@ -322,7 +390,7 @@ for s in [34, 36, 37, 38]: #range(15, 27): GFS_Lat = 13.1075 Obs_geo['inArea'] = Obs_geo.within(siteShp.iloc[4, 1]) - # Set min and max times using conductivity + # Set min and max times using shape if Obs_geo['inArea'].any(): # First and last times from area in shapefile minTime = Obs_geo[Obs_geo['inArea']==True].index[0] @@ -455,240 +523,319 @@ for s in [34, 36, 37, 38]: #range(15, 27): met_tp = -999 met_mwd = -999 - #%% Add PSU and depth units - Obs_geo['PSU'] = gsw.conversions.SP_from_C(Obs_geo['CH3:Conductivity(mS/cm)'].values, - Obs_geo['CH4:Temperature(degC)'].values, - Obs_geo['CH0:Pressure(dbar)'].values) + #%% Add PSU and depth units if not present + + if ('Salinity ' not in Obs_geo.columns) and (JFE_File is None): + Obs_geo['PSU'] = gsw.conversions.SP_from_C(Obs_geo['CH3:Conductivity(mS/cm)'].values, + Obs_geo['CH4:Temperature(degC)'].values, + Obs_geo['CH0:Pressure(dbar)'].values) - Obs_geo['Depth'] = (Obs_geo['CH0:Pressure(dbar)'].astype(float)*10000)/(gsw.rho(Obs_geo['PSU'].values, - Obs_geo['CH1:Temperature(degC)'].values, - Obs_geo['CH0:Pressure(dbar)'].values) * 9.81) + Obs_geo['Depth'] = (Obs_geo['CH0:Pressure(dbar)'].astype(float)*10000)/(gsw.rho(Obs_geo['PSU'].values, + Obs_geo['CH1:Temperature(degC)'].values, + Obs_geo['CH0:Pressure(dbar)'].values) * 9.81) + elif (JFE_File is not None): + # Rename Salinity to PSU + Obs_geo.rename(columns={'Salinity ': 'PSU'}, inplace=True) + # Fix depth name + Obs_geo.rename(columns={'Depth ': 'Depth'}, inplace=True) + + else: # Fix turbidity name + Obs_geo.rename(columns={'Turbidity ': 'CH6:Turbidity(NTU)'}, inplace=True) + + + #%% Save data as geojson + Obs_geo.loc[minTime:maxTime, :].to_file( + importPaths[s] + '/MergedGeoDataFrame.geojson', driver='GeoJSON') #%% Plot time series for Geo data - fontprops = fm.FontProperties(size=25) + if figureSave: + fontprops = fm.FontProperties(size=25) + if (JFE_File is None) and (RBR_File is not None): + fig, axesTMP = plt.subplots(nrows=1, ncols=1, figsize=(19, 5), constrained_layout=True) + axes = [] + axes.append(axesTMP) - fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(19, 25), constrained_layout=True) - dataTable = np.zeros([14, 4]) - roundIDX = [1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1] + dataTable = np.zeros([9, 4]) + roundIDX = [1, 1, 0, 1, 1, 1, 1, 1] - params = ['Depth', 'PSU', 'CH8:Oxygen_Saturation(%)', 'CH1:Temperature(degC)', - 'CH6:Turbidity(NTU)', 'CH2:Chlorophyll_a(ppb)'] - paramName = ['Depth [m]', 'Salinity [PSU]', 'Dissolved O₂ sat [%]', 'Temperature [degC]', - 'Turbidity [FTU]', 'Chl-a [ppb]'] + params = ['CH6:Turbidity(NTU)'] + paramName = ['Turbidity [FTU]'] + parmCmap = [cmo.turbid] + RBRparamMin = [0.0] + RBRparamMax = [60.0] + paramMin = [0] + paramMax = [60] + order = 12 + smoothIDX = [60] + smoothPeriod = '1 Minute Average' - parmCmap = [cmo.deep, 'cividis', cmo.dense, cmo.thermal, cmo.turbid, cmo.algae] - # paramMin = [0.0, 34.0, 32.5, 25.0, 0, 0] - # paramMax = [1.0, 36.0, 34.0, 31.0, 20.0, 1.0] - paramMin = [0.0, 32.0, 100, 24.0, 0, 0] - paramMax = [6, 36.5, 130, 34.0, 150.0, 12000] - - - fig.patch.set_facecolor('white') - # fig.tight_layout(pad=1.05) - - fontprops = fm.FontProperties(size=25) - - dataTable[0, 0] = met_temp.m.item(0)-272.15 - dataTable[1, 0] = met_wind - dataTable[2, 0] = met_wdir - dataTable[3, 0] = met_prep.m.item(0) - dataTable[4, 0] = met_hmo - dataTable[5, 0] = met_tp - dataTable[6, 0] = met_mwd - - ilocs_max = [] - ilocs_max_pts = [] - OBS_mask = [] - - for paramIDX, param in enumerate(params): - Obs_geo.loc[minTime:maxTime, param].plot( - ax=axes[paramIDX], label='1 Second Observations', color='lightgrey') # Note the space in the col name - - # Create mask for RBR data based on time and parameter minimum - OBS_mask.append(((Obs_geo.index>minTime) & - (Obs_geo.indexparamMin[paramIDX]))) - - OBS_smoothed = Obs_geo.loc[OBS_mask[paramIDX], param].rolling( - 12, win_type='nuttall',center=True).mean() - - OBS_smoothed.plot( - ax=axes[paramIDX], label='10 Second Average', color='black', - linewidth=3) - - # Find the local maximums for Turbidity - if param == 'CH6:Turbidity(NTU)': - ##### Adjust "order" to control the number of turbidity points that are selected - ilocs_max.append(argrelextrema(OBS_smoothed.values, - np.greater_equal, order=12, mode='wrap')[0]) - - # Add start and end points? - # ilocs_max = np.insert(ilocs_max, 0, 10) - # ilocs_max[-1] = len(OBS_smoothed.values) - - ilocs_max_pts.append(OBS_smoothed.iloc[ilocs_max[paramIDX]].index.values) - - # Add labels if GPS data is available - if GPS_File is not None: - # axes[paramIDX].scatter(OBS_smoothed.iloc[ - # ilocs_max[paramIDX]].index, np.ones(len(ilocs_max[paramIDX])) * 30, 75, - # color='blue') - for a in range(0, len(ilocs_max[paramIDX])): - axes[paramIDX].annotate(str(a+1), (ilocs_max_pts[paramIDX][a], 12), fontsize=30) else: - ilocs_max.append(None) - ilocs_max_pts.append(None) + fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(19, 25), constrained_layout=True) + dataTable = np.zeros([14, 4]) + params = ['Depth', 'PSU', 'CH8:Oxygen_Saturation(%)', 'CH1:Temperature(degC)', + 'CH6:Turbidity(NTU)', 'CH2:Chlorophyll_a(ppb)'] + paramName = ['Depth [m]', 'Salinity [PSU]', 'Dissolved O₂ sat [%]', 'Temperature [degC]', + 'Turbidity [FTU]', 'Chl-a [ppb]'] - dataTable[paramIDX+7, 0] = OBS_smoothed.mean(skipna=True) - dataTable[paramIDX+7, 1] = OBS_smoothed.std(skipna=True) - dataTable[paramIDX+7, 2] = max(OBS_smoothed.min(skipna=True), 0) - dataTable[paramIDX+7, 3] = OBS_smoothed.max(skipna=True) + parmCmap = [cmo.deep, 'cividis', cmo.dense, cmo.thermal, cmo.turbid, cmo.algae] + roundIDX = [1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1] + if s<12: # Original sensor + # paramMin = [0.0, 34.0, 32.5, 25.0, 0, 0] + # paramMax = [6.0, 36.0, 34.0, 31.0, 20.0, 1.0] + paramMin = [0.0, 32.0, 32.5, 24.0, 0, 0] + paramMax = [2.0, 36.0, 34.0, 34.0, 20.0, 1.0] + order = 12 + smoothIDX = [60, 60, 60, 20, 20, 20] + smoothPeriod = '1 Minute Average' + else: # New sensor + paramMin = [0.0, 32.0, 100, 24.0, 0, 0] + paramMax = [6, 36.5, 130, 34.0, 150.0, 12000] + order = 60 + smoothIDX = [12, 12, 12, 12, 12, 12] + smoothPeriod = '10 Second Average' - axes[paramIDX].set_ylabel(paramName[paramIDX]) - axes[paramIDX].set_title(paramName[paramIDX]) - axes[paramIDX].set_xlabel('') - axes[paramIDX].set_ylim(paramMin[paramIDX], paramMax[paramIDX]) - axes[paramIDX].legend(loc='upper right') - # axes[paramIDX].set_xlabel(minTime.strftime('%B %#d, %Y')) - # Formate Data Table - dataTableFormat_mean = [] - dataTableFormat_maxmin = [] - for d in range(0, len(roundIDX)): - dataTableFormat_mean.append(round(dataTable[d, 0], roundIDX[d])) - if dataTable[d, 3] == 0: - dataTableFormat_maxmin.append('--') + + fig.patch.set_facecolor('white') + # fig.tight_layout(pad=1.05) + + fontprops = fm.FontProperties(size=25) + + dataTable[0, 0] = met_temp.m.item(0)-272.15 + dataTable[1, 0] = met_wind + dataTable[2, 0] = met_wdir + dataTable[3, 0] = met_prep.m.item(0) + dataTable[4, 0] = met_hmo + dataTable[5, 0] = met_tp + dataTable[6, 0] = met_mwd + + ilocs_max = [] + ilocs_max_pts = [] + OBS_mask = [] + + for paramIDX, param in enumerate(params): + Obs_geo.loc[minTime:maxTime, param].plot( + ax=axes[paramIDX], label='1 Second Observations', color='lightgrey') # Note the space in the col name + + # Create mask for RBR data based on time and parameter minimum + OBS_mask.append(((Obs_geo.index>minTime) & + (Obs_geo.indexparamMin[paramIDX]))) + + OBS_smoothed = Obs_geo.loc[OBS_mask[paramIDX], param].rolling( + 12, win_type='nuttall',center=True).mean() + + OBS_smoothed.plot( + ax=axes[paramIDX], label=smoothPeriod, color='black', + linewidth=3) + + # Find the local maximums for Turbidity + if param == 'CH6:Turbidity(NTU)': + ##### Adjust "order" to control the number of turbidity points that are selected + + # Find an order that creates between 5-10 points + while (len(argrelextrema(OBS_smoothed.values, np.greater_equal, order=order, mode='wrap')[0]) < 5) or\ + (len(argrelextrema(OBS_smoothed.values, np.greater_equal, order=order, mode='wrap')[0]) > 10) or\ + (order == 0): + + if len(argrelextrema(OBS_smoothed.values, np.greater_equal, order=order, mode='wrap')[0]) < 5: + order = order - 5 + elif len(argrelextrema(OBS_smoothed.values, np.greater_equal, order=order, mode='wrap')[0]) > 10: + order = order + 6 + + if order == 0: + order = 1 + + ilocs_max.append(argrelextrema(OBS_smoothed.values, + np.greater_equal, order=order, mode='wrap')[0]) + + + # Add start and end points? + # ilocs_max = np.insert(ilocs_max, 0, 10) + # ilocs_max[-1] = len(OBS_smoothed.values) + + ilocs_max_pts.append(OBS_smoothed.iloc[ilocs_max[paramIDX]].index.values) + + # Add labels if GPS data is available + if GPS_File is not None: + # axes[paramIDX].scatter(OBS_smoothed.iloc[ + # ilocs_max[paramIDX]].index, np.ones(len(ilocs_max[paramIDX])) * 30, 75, + # color='blue') + for a in range(0, len(ilocs_max[paramIDX])): + axes[paramIDX].annotate(str(a+1), (ilocs_max_pts[paramIDX][a], 12), fontsize=30) + else: + ilocs_max.append(None) + ilocs_max_pts.append(None) + + dataTable[paramIDX+7, 0] = OBS_smoothed.mean(skipna=True) + dataTable[paramIDX+7, 1] = OBS_smoothed.std(skipna=True) + dataTable[paramIDX+7, 2] = max(OBS_smoothed.min(skipna=True), 0) + dataTable[paramIDX+7, 3] = OBS_smoothed.max(skipna=True) + + axes[paramIDX].set_ylabel(paramName[paramIDX]) + axes[paramIDX].set_title(paramName[paramIDX]) + axes[paramIDX].set_xlabel('') + axes[paramIDX].set_ylim(paramMin[paramIDX], paramMax[paramIDX]) + axes[paramIDX].legend(loc='upper right') + + # axes[paramIDX].set_xlabel(minTime.strftime('%B %#d, %Y')) + + # Formate Data Table + dataTableFormat_mean = [] + dataTableFormat_maxmin = [] + for d in range(0, len(roundIDX)): + dataTableFormat_mean.append(round(dataTable[d, 0], roundIDX[d])) + if dataTable[d, 3] == 0: + dataTableFormat_maxmin.append('--') + else: + dataTableFormat_maxmin.append(str(round(dataTable[d, 2], roundIDX[d])) + ' / ' + str(round(dataTable[d, 3], roundIDX[d]))) + + + dfOut = pd.DataFrame(dataTable[:, :]) + dfOutFormat = pd.DataFrame([dataTableFormat_mean, dataTableFormat_maxmin]).transpose() + + # Change the column names + dfOut.columns =['Mean', 'Standard Deviation', 'Min', 'Max'] + dfOutFormat.columns = [np.array([minTime.strftime('%B %#d, %Y'), minTime.strftime('%B %#d, %Y')]), + np.array(['Mean', 'Min / Max'])] + # Change the row indexes + + if (JFE_File is None) and (RBR_File is not None): + dfOut.iloc[8, 0] = minTime.strftime('%B %#d, %Y') + rowNames = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]', + 'Significant Wave Height Offshore [m]', 'Peak Wave Period Offshore [s]', + 'Mean Wave Direction Offshore [deg]', + 'Turbidity [NTU]', 'Observation Date'] + dfOut.index = rowNames + dfOutFormat.index = rowNames[0:-1] + else: - dataTableFormat_maxmin.append(str(round(dataTable[d, 2], roundIDX[d])) + ' / ' + str(round(dataTable[d, 3], roundIDX[d]))) + dfOut.iloc[13, 0] = minTime.strftime('%B %#d, %Y') + rowNames = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]', + 'Significant Wave Height Offshore [m]' ,'Peak Wave Period Offshore [s]', + 'Mean Wave Direction Offshore [deg]', 'Depth [m]', 'Salinity [PSU]', + 'Dissolved O₂ saturation [%]', 'Temperature [degC]', + 'Turbidity [FTU]', 'Chl-a [ppb]', 'Observation Date'] - dfOut = pd.DataFrame(dataTable[:, :]) - dfOutFormat = pd.DataFrame([dataTableFormat_mean, dataTableFormat_maxmin]).transpose() + dfOut.index = rowNames + dfOutFormat.index = rowNames[0:-1] - # Change the column names - dfOut.columns =['Mean', 'Standard Deviation', 'Min', 'Max'] - dfOutFormat.columns = [np.array([minTime.strftime('%B %#d, %Y'), minTime.strftime('%B %#d, %Y')]), - np.array(['Mean', 'Min / Max'])] - # Change the row indexes - dfOut.iloc[13, 0] = minTime.strftime('%B %#d, %Y') - rowNames = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]', - 'Significant Wave Height Offshore [m]' ,'Peak Wave Period Offshore [s]', - 'Mean Wave Direction Offshore [deg]', 'Depth [m]', 'Salinity [PSU]', - 'Dissolved O₂ saturation [%]', 'Temperature [degC]', - 'Turbidity [FTU]', 'Chl-a [ppb]', 'Observation Date'] + fig.suptitle(siteName + ', ' + minTime.strftime('%B %#d, %Y') + ' (' + timeLabel + ')', fontsize=30) + plt.show() - dfOut.index = rowNames - dfOutFormat.index = rowNames[0:-1] + # outPath = '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + timeLabels[s] + outPath = importPaths[s] - fig.suptitle(siteName + ', ' + minTime.strftime('%B %#d, %Y') + ' (' + timeLabel + ')', fontsize=30) + if not os.path.exists(outPath): + os.mkdir(outPath) - plt.show() + dfOut.to_excel(outPath + '/Summary_Stats_' + siteName + '_RevB.xlsx') + dfOutFormat.to_excel(outPath + '/Summary_StatsFormat_' + siteName + '_RevB.xlsx') - # outPath = '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + timeLabels[s] - outPath = importPaths[s] + dfOut.to_csv(outPath + '/Summary_Stats_' + siteName + '_RevB.csv') - if not os.path.exists(outPath): - os.mkdir(outPath) + if not os.path.exists(outPath + '/Figures'): + os.mkdir(outPath + '/Figures') - dfOut.to_excel(outPath + '/Summary_Stats_' + siteName + '.xlsx') - dfOutFormat.to_excel(outPath + '/Summary_StatsFormat_' + siteName + '.xlsx') - - dfOut.to_csv(outPath + '/Summary_Stats_' + siteName + '.csv') - - if not os.path.exists(outPath + '/Figures'): - os.mkdir(outPath + '/Figures') - - fig.savefig(outPath + '/Figures/SummaryTimeSeries_' + siteName + '.pdf', - bbox_inches='tight') - fig.savefig(outPath + '/Figures/SummaryTimeSeries_' + siteName + '.png', - bbox_inches='tight', dpi=500) + fig.savefig(outPath + '/Figures/SummaryTimeSeries_' + siteName + '_RevB.pdf', + bbox_inches='tight') + fig.savefig(outPath + '/Figures/SummaryTimeSeries_' + siteName + '_RevB.png', + bbox_inches='tight', dpi=500) #%% Plot Maps - fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(19, 25), constrained_layout=True) - ax = axes.flat + if figureSave: + if (JFE_File is None) and (RBR_File is not None): + # Only Turbidity Data + fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 9), constrained_layout=True) + ax = [] + ax.append(axes) + else: + fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(19, 25), constrained_layout=True) + ax = axes.flat - fig.patch.set_facecolor('white') - # fig.tight_layout(pad=1.05) + fig.patch.set_facecolor('white') + # fig.tight_layout(pad=1.05) - fontprops = fm.FontProperties(size=25) - x, y, arrow_length = 0.95, 0.93, 0.20 - plt.rcParams.update({'font.size': 22}) + fontprops = fm.FontProperties(size=25) + x, y, arrow_length = 0.95, 0.93, 0.20 + plt.rcParams.update({'font.size': 22}) - axXlimTT = (Obs_geo.loc[minTime:maxTime].geometry.x.min()-100, - Obs_geo.loc[minTime:maxTime].geometry.x.max()+100) - axYlimTT = (Obs_geo.loc[minTime:maxTime].geometry.y.min()-100, - Obs_geo.loc[minTime:maxTime].geometry.y.max()+100) + axXlimTT = (Obs_geo.loc[minTime:maxTime].geometry.x.min()-100, + Obs_geo.loc[minTime:maxTime].geometry.x.max()+100) + axYlimTT = (Obs_geo.loc[minTime:maxTime].geometry.y.min()-100, + Obs_geo.loc[minTime:maxTime].geometry.y.max()+100) - plt.setp(axes, xlim=axXlim, ylim=axYlim) - # plt.setp(axes, xlim=axXlimTT, ylim=axYlimTT) + plt.setp(axes, xlim=axXlim, ylim=axYlim) + # plt.setp(axes, xlim=axXlimTT, ylim=axYlimTT) - # Plot the RBR observations - # Salinity - for paramIDX, param in enumerate(params): - Obs_geo.loc[minTime:maxTime].plot( - column=param, ax=ax[paramIDX], vmin=paramMin[paramIDX], vmax=paramMax[paramIDX], - legend=True, legend_kwds={'label': paramName[paramIDX]}, - cmap=parmCmap[paramIDX], markersize=20) + # Plot the RBR observations + # Salinity + for paramIDX, param in enumerate(params): + Obs_geo.loc[minTime:maxTime].plot( + column=param, ax=ax[paramIDX], vmin=paramMin[paramIDX], vmax=paramMax[paramIDX], + legend=True, legend_kwds={'label': paramName[paramIDX]}, + cmap=parmCmap[paramIDX], markersize=20) - ctx.add_basemap(ax[paramIDX], source=mapbox, crs='EPSG:32621') + ctx.add_basemap(ax[paramIDX], source=mapbox, crs='EPSG:32621') - # Add time labels - # plt.scatter(RBR_Obs_geo.loc[RBR_mask].iloc[ - # ilocs_max, :].geometry.x, - # RBR_Obs_geo.loc[RBR_mask].iloc[ - # ilocs_max, :].geometry.y, 75, marker='o', color='black') + # Add time labels + # plt.scatter(RBR_Obs_geo.loc[RBR_mask].iloc[ + # ilocs_max, :].geometry.x, + # RBR_Obs_geo.loc[RBR_mask].iloc[ + # ilocs_max, :].geometry.y, 75, marker='o', color='black') - if (not Obs_geo.geometry.isnull().all()) &\ - (GPS_File is not None) & (param == 'CH6:Turbidity(NTU)'): - for a in range(0, len(ilocs_max[paramIDX])): - ax[paramIDX].annotate(str(a + 1), (Obs_geo.loc[OBS_mask[paramIDX]].iloc[ - ilocs_max[paramIDX][a], :].geometry.x, - Obs_geo.loc[OBS_mask[paramIDX]].iloc[ - ilocs_max[paramIDX][a], :].geometry.y), fontsize=30) + if (not Obs_geo.geometry.isnull().all()) &\ + (GPS_File is not None) & (param == 'CH6:Turbidity(NTU)'): + for a in range(0, len(ilocs_max[paramIDX])): + ax[paramIDX].annotate(str(a + 1), (Obs_geo.loc[OBS_mask[paramIDX]].iloc[ + ilocs_max[paramIDX][a], :].geometry.x, + Obs_geo.loc[OBS_mask[paramIDX]].iloc[ + ilocs_max[paramIDX][a], :].geometry.y), fontsize=30) - ax[paramIDX].set_title(paramName[paramIDX]) - # ax[paramIDX].set_ylabel('UTM 21N [m]') - # ax[paramIDX].set_xlabel('UTM 21N [m]') - ax[paramIDX].locator_params(axis='y', nbins=3) - ax[paramIDX].ticklabel_format(useOffset=False, style='plain', axis='both') + ax[paramIDX].set_title(paramName[paramIDX]) + # ax[paramIDX].set_ylabel('UTM 21N [m]') + # ax[paramIDX].set_xlabel('UTM 21N [m]') + ax[paramIDX].locator_params(axis='y', nbins=3) + ax[paramIDX].ticklabel_format(useOffset=False, style='plain', axis='both') - ax[paramIDX].get_xaxis().set_ticks([]) - ax[paramIDX].get_yaxis().set_ticks([]) + ax[paramIDX].get_xaxis().set_ticks([]) + ax[paramIDX].get_yaxis().set_ticks([]) - #Add scale-bar - scalebar = AnchoredSizeBar(ax[paramIDX].transData, - 100, '100 m', 'lower right', pad=0.5, size_vertical=10, fontproperties=fontprops) - ax[paramIDX].add_artist(scalebar) - ax[paramIDX].annotate('N', xy=(x, y), xytext=(x, y-arrow_length), - arrowprops=dict(facecolor='black', width=6, headwidth=30), - ha='center', va='center', fontsize=35, - xycoords=ax[paramIDX].transAxes) + #Add scale-bar + scalebar = AnchoredSizeBar(ax[paramIDX].transData, + 100, '100 m', 'lower right', pad=0.5, size_vertical=10, fontproperties=fontprops) + ax[paramIDX].add_artist(scalebar) + ax[paramIDX].annotate('N', xy=(x, y), xytext=(x, y-arrow_length), + arrowprops=dict(facecolor='black', width=6, headwidth=30), + ha='center', va='center', fontsize=35, + xycoords=ax[paramIDX].transAxes) - fig.suptitle(siteName + ', ' + minTime.strftime('%b %#d, %Y') + ' (' + timeLabel + ')', fontsize=30) + fig.suptitle(siteName + ', ' + minTime.strftime('%b %#d, %Y') + ' (' + timeLabel + ')', fontsize=30) - plt.show() + plt.show() - #outPath = '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + timeLabels[s] - outPath = importPaths[s] + #outPath = '//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + timeLabels[s] + outPath = importPaths[s] - if not os.path.exists(outPath): - os.mkdir(outPath) + if not os.path.exists(outPath): + os.mkdir(outPath) - if not os.path.exists(outPath + '/Figures'): - os.mkdir(outPath + '/Figures') + if not os.path.exists(outPath + '/Figures'): + os.mkdir(outPath + '/Figures') - fig.savefig(outPath + '/Figures/SummaryMap_' + timeLabels[s] + '_' + siteName + '.pdf', - bbox_inches='tight') + fig.savefig(outPath + '/Figures/SummaryMap_' + timeLabels[s] + '_' + siteName + '_RevB.pdf', + bbox_inches='tight') - fig.savefig(outPath + '/Figures/SummaryMap_' + timeLabels[s] + '_' + siteName + '.png', - bbox_inches='tight', dpi=500) + fig.savefig(outPath + '/Figures/SummaryMap_' + timeLabels[s] + '_' + siteName + '_RevB.png', + bbox_inches='tight', dpi=500) #%% Summary Sheet plotIDXsLoop = [] @@ -788,3 +935,150 @@ for i in range(2, 3): bbox_inches='tight') fig.savefig('//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + siteName + '.png', bbox_inches='tight', dpi=500) + +#%% Save tracks for one location to shapefile + +mergedLines = [] + +# Loop through all samples +for s in range(0, len(siteNames)): + siteName = siteNames[s] + + if siteName == 'Crane': + importPath = importPaths[s] + + # Load in track data + locationGDF = gp.read_file(importPath + '/MergedGeoDataFrame.geojson') + + # Convert points to linestring and add to list + mergedLines.append(LineString([[a.x, a.y] for a in locationGDF.geometry.values])) + + # Track progress + print(s) + +# Convert list of tracks to geodataframe +mergedLines_gdf = gp.GeoDataFrame(geometry=mergedLines, crs="EPSG:32621") + +# Save to shapefile +mergedLines_gdf.to_file('//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/Crane_Tracks.shp') + + + +#%% Read GeoJSON file for each site and extract data at two points + +# Plot trubidity at each point through time +fig, axes = plt.subplots(1, 2, figsize=(15, 6)) + +# Define points of interest as geodataframe 5 m from point + +meanTurbA = [] +meanTurbB = [] +meanTime = [] + +# Set plotting site +# plotSite = 'Great House' +# plotSite = 'Caribee' +# plotSite = 'Greensleeves' +# plotSite = 'Old Queens Fort' +plotSite = 'Crane' + + +# Mapping parameter +if plotSite == 'Great House': + axXlim = (213210.7529575412, 213562.64172686986) + axYlim = (1464769.2243017585, 1465135.2219089477) + A_pt = [213316.46, 1464972.00] + B_pt = [213329.09, 1464898.30] + +elif plotSite == 'Greensleeves': + axXlim = (213269.99233348924, 213648.1643157148) + axYlim = (1463378.1020314451, 1463950.5442048472) + A_pt = [213391.92, 1463659.98] + B_pt = [213518.95, 1463517.72] + +elif plotSite == 'Old Queens Fort': + axXlim = (213368.59866770002, 213745.6997016811) + axYlim = (1460192.707288096, 1460672.371780407) + A_pt = [213527.75, 1460354.13] + B_pt = [213510.13, 1460565.55] + +elif plotSite == 'Crane': + axXlim = (234835, 235507) + axYlim = (1449966, 1450580) + A_pt = [235287.01, 1450171.33] + B_pt = [235300.81, 1450233.21] + +elif plotSite == 'Caribee': + axXlim = (217929, 218345) + axYlim = (1446528, 1446818) + A_pt = [218120.54, 1446683.057] + B_pt = [218250.99, 1446676.20] + +# Buffer to 20 m +A_GDF = gp.GeoDataFrame(geometry=gp.points_from_xy([A_pt[0]], [A_pt[1]]), crs="EPSG:32621").buffer(22) +B_GDF = gp.GeoDataFrame(geometry=gp.points_from_xy([B_pt[0]], [B_pt[1]]), crs="EPSG:32621").buffer(22) + +# Set Map Bounds +axes[0].set_xlim(axXlim) +axes[0].set_ylim(axYlim) + +# Add Base Map +ctx.add_basemap(axes[0], source=mapbox, crs='EPSG:32621') + +# Loop through all samples +for s in range(0, len(siteNames)): + + siteName = siteNames[s] + + if siteName == plotSite: + importPath = importPaths[s] + + # Load in track data + locationGDF = gp.read_file(importPath + '/MergedGeoDataFrame.geojson') + + # Find points within the point of interest buffer + # Average turbidity + meanTurbA.append(locationGDF.loc[locationGDF.within(A_GDF.iloc[0]), 'CH6:Turbidity(NTU)'].mean()) + meanTurbB.append(locationGDF.loc[locationGDF.within(B_GDF.iloc[0]), 'CH6:Turbidity(NTU)'].mean()) + meanTime.append(locationGDF.loc[0, 'DateTime']) + + + # Plot Track + locationGDF.plot(ax=axes[0]) + +# Plot turbidity locations as empty circles +A_GDF.plot(ax=axes[0], facecolor="none", edgecolor='k', linewidth=5) +B_GDF.plot(ax=axes[0], facecolor="none", edgecolor='k', linewidth=5) + +# Add text labels to turbidity locations +axes[0].text(A_pt[0]+30, A_pt[1]-10, 'A', fontsize=25) +axes[0].text(B_pt[0]+30, B_pt[1]-10, 'B', fontsize=25) + +# Remove ticks +axes[0].get_xaxis().set_ticks([]) +axes[0].get_yaxis().set_ticks([]) + +# Plot turbidity +axes[1].plot(pd.to_datetime(meanTime), meanTurbA, label='Turbidity Point A') +axes[1].plot(pd.to_datetime(meanTime), meanTurbB, label='Turbidity Point B') + +# Add Legend +axes[1].legend() + +# Add axis labels +axes[1].set_xlabel('Date') +axes[1].set_ylabel('Turbidity (NTU)') + +# Limit x axis to 4 ticks +axes[1].xaxis.set_major_locator(plt.MaxNLocator(4)) + +# Format x axis dates +axes[1].xaxis.set_major_formatter(mdates.DateFormatter('%b %y')) + +# Add plot title +plt.suptitle(plotSite) + +plt.show() + +# Save figure +fig.savefig('//srv-ott3.baird.com/Projects/INNOVATION/Phase97 Barbados Research/05_Analyses/' + plotSite + '_Turbidity_Time.png', dpi=300, bbox_inches='tight') \ No newline at end of file diff --git a/NTC_DFM/Delft3DBoundsCompare.py b/NTC_DFM/Delft3DBoundsCompare.py index e69de29..7a912f4 100644 --- a/NTC_DFM/Delft3DBoundsCompare.py +++ b/NTC_DFM/Delft3DBoundsCompare.py @@ -0,0 +1,83 @@ +#%% Scratch file to compare Delft3D input bounds +# Alexander Rey, Nay 16, 2023 + +#%% Import modules +import os +import numpy as np +import matplotlib.pyplot as plt +plt.close('all') +import dfm_tools as dfmt +import hydrolib.core.dflowfm as hcdfm +import pandas as pd +import datetime +from matplotlib import dates as mpl_dates + +#%% Read in Salinity data +file_bc = r'\\ott-diomedes\Projects\11934.201_Newtown_Creek_TPP\Delft3D\Run169_2014HD_Sed_OrthoB_pliz_WorkingNoteNoMin_Wind.dsproj_data\FlowFM\input\Salinity.bc' + +#Load .bc-file using HydroLib object ForcingModel. +m = hcdfm.ForcingModel(file_bc) +ForcingModel_object_out = hcdfm.ForcingModel() + +# Create xarray object for SouthBound +salinityN_xr = dfmt.forcinglike_to_Dataset(m.forcing[0], convertnan=True) +salinityS_xr = dfmt.forcinglike_to_Dataset(m.forcing[2], convertnan=True) + +#%% Read in Temperature data +file_bc = r'\\ott-diomedes\Projects\11934.201_Newtown_Creek_TPP\Delft3D\Run169_2014HD_Sed_OrthoB_pliz_WorkingNoteNoMin_Wind.dsproj_data\FlowFM\input\Temperature.bc' + +#Load .bc-file using HydroLib object ForcingModel. +m = hcdfm.ForcingModel(file_bc) +ForcingModel_object_out = hcdfm.ForcingModel() + +# Create xarray object for NorthBound and SouthBound +temperatureN_xr = dfmt.forcinglike_to_Dataset(m.forcing[0], convertnan=True) +temperatureS_xr = dfmt.forcinglike_to_Dataset(m.forcing[2], convertnan=True) + + +#%% Read in Water Level data +file_bc = r'\\ott-diomedes\Projects\11934.201_Newtown_Creek_TPP\Delft3D\Run169_2014HD_Sed_OrthoB_pliz_WorkingNoteNoMin_Wind.dsproj_data\FlowFM\input\WaterLevel.bc' + +#Load .bc-file using HydroLib object ForcingModel. +m = hcdfm.ForcingModel(file_bc) +ForcingModel_object_out = hcdfm.ForcingModel() + +# Create xarray object for NorthBound +waterlevelN_xr = dfmt.forcinglike_to_Dataset(m.forcing[0], convertnan=True) +waterlevelS_xr = dfmt.forcinglike_to_Dataset(m.forcing[2], convertnan=True) + + +#%% Plot Bound Data + +# Set Time Bounds +dateMin = pd.to_datetime(datetime.datetime(2015, 7, 18, 0, 0, 0, 0)) +dateMax = pd.to_datetime(datetime.datetime(2015, 7, 22, 0, 0, 0, 0)) + +fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(12, 6), sharex=True) + +# Plot a line for salinity for all times for midpoint layer +salinityN_xr.salinitybnd.isel(depth=9).plot.line(ax=ax[0], linestyle='--', color='r', label='North Top') +salinityS_xr.salinitybnd.isel(depth=9).plot.line(ax=ax[0], linestyle='-', color='r', label='South Top') +salinityN_xr.salinitybnd.isel(depth=0).plot.line(ax=ax[0], linestyle='--', color='b', label='North Bottom') +salinityS_xr.salinitybnd.isel(depth=0).plot.line(ax=ax[0], linestyle='-', color='b', label='South Bottom') + +temperatureN_xr.temperaturebnd.isel(depth=9).plot.line(ax=ax[1], linestyle='--', color='r') +temperatureS_xr.temperaturebnd.isel(depth=9).plot.line(ax=ax[1], linestyle='-', color='r') +temperatureN_xr.temperaturebnd.isel(depth=0).plot.line(ax=ax[1], linestyle='--', color='b') +temperatureS_xr.temperaturebnd.isel(depth=0).plot.line(ax=ax[1], linestyle='-', color='b') + +waterlevelN_xr.waterlevelbnd.plot.line(ax=ax[2], linestyle='--', color='r') +waterlevelS_xr.waterlevelbnd.plot.line(ax=ax[2], linestyle='-', color='r') + +# Set y axis limits +ax[0].set_ylim(12, 27) +ax[1].set_ylim(15, 25) +ax[2].set_ylim(-2, 2) + +ax[2].set_xlim(dateMin, dateMax) +ax[2].xaxis.set_major_formatter( + mpl_dates.ConciseDateFormatter(ax[2].xaxis.get_major_locator())) + +ax[0].legend() + +plt.show() diff --git a/NTC_DFM/NTC_PlottingD3D_2023.py b/NTC_DFM/NTC_PlottingD3D_2023.py index 4b47776..fc751f6 100644 --- a/NTC_DFM/NTC_PlottingD3D_2023.py +++ b/NTC_DFM/NTC_PlottingD3D_2023.py @@ -97,7 +97,7 @@ for i in modelPlot: Delft_WL[i] = data_xr.waterlevel.sel(stations=['West_WL', 'East_WL']).to_pandas() # Put in UTC if needed - if i == 104: + if i == 104 or i == 105 or i >= 135: Delft_WL[i] = Delft_WL[i].tz_localize( pytz.timezone('EST')).tz_convert(pytz.utc) else: diff --git a/NTC_DFM/NTC_SedTestPlot.py b/NTC_DFM/NTC_SedTestPlot.py index 2cf232c..4efd0b6 100644 --- a/NTC_DFM/NTC_SedTestPlot.py +++ b/NTC_DFM/NTC_SedTestPlot.py @@ -20,8 +20,18 @@ import matplotlib.animation as animation import dfm_tools as dfmt +# Geotiff +import cartopy.crs as ccrs +import contextily as ctx +import pathlib as pl + +from shapely.geometry import Point, MultiPoint +import alphashape + +import rioxarray import xarray as xr + import cartopy.crs as ccrs import contextily as ctx from dfm_tools.regulargrid import scatter_to_regulargrid @@ -41,13 +51,14 @@ histPath = "FlowFM/dflowfm/output/FlowFM_his.nc" #%% Import Model results for static images # Define which models to import -modelPlot = [120, 121, 129, 130, 131, 132, 133, 134] +modelPlot = [144] d3dfm_DataArray = [None] * (max(modelPlot)+1) sedMass = [None] * (max(modelPlot)+1) # Available Mass of Sediment erodeSed = [None] * (max(modelPlot)+1) # Erosion and deposition initialBed = [None] * (max(modelPlot)+1) # Erosion and deposition modelCRS = [None] * (max(modelPlot)+1) # Model CRS +sedConc = [None] * (max(modelPlot)+1) # Model CRS for i in modelPlot: file_nc_map = Path(runLog['Run Location'][i]) / dataPath @@ -70,6 +81,9 @@ for i in modelPlot: erodeSed[i] = erodeSed[i].compute() initialBed[i] = initialBed[i].compute() + # Import concentration of sediment 2 at timestep 95 and layer 10 + sedConc[i] = d3dfm_DataArray[i].mesh2d_sedfrac_concentration[94, 2, :, 9].compute() + modelCRS[i] = vars_pd['EPSG_code']['projected_coordinate_system'] #%% Define axis limits @@ -99,7 +113,7 @@ axLim_Xmin.append(305500) axLim_Xmax.append(306800) axLim_Ymin.append(60000) -axLim_Ymax.append(62200) +axLim_Ymax.append(62500) # Zoom far trib axLim_Xmin.append(305800) @@ -161,8 +175,8 @@ for i in modelPlot: # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/ErodeDep_' + runLog['Run'][i] + '.png', # bbox_inches='tight', dpi=300) - fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/InitialBed_' + runLog['Run'][i] + '.png', - bbox_inches='tight', dpi=300) + # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/InitialBed_' + runLog['Run'][i] + '.png', + # bbox_inches='tight', dpi=300) #%% Plot using DFM functions #Cmap Path @@ -425,3 +439,26 @@ for m in modelPlot: # Save video writer.finish() + + +#%% Save Sediment Geotiff + +# Create mask within bounds +boundSelect = 2 +croppedGridMask = (d3dfm_DataArray[i].mesh2d_face_x >= axLim_Xmin[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_x <= axLim_Xmax[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_y >= axLim_Ymin[boundSelect]) &\ + (d3dfm_DataArray[i].mesh2d_face_y <= axLim_Ymax[boundSelect]) +croppedGridMask = croppedGridMask.compute() + +# Step through time steps +for sedTstep in range(105, 120): + # Rasterize sediment concentration within bounds + sedConcTif = dfmt.rasterize_ugrid(d3dfm_DataArray[i]['mesh2d_sedfrac_concentration'].isel( + time=sedTstep, mesh2d_nLayers=4, nSedSus=0).where(croppedGridMask, drop=True).to_dataset(), + resolution=2) + + # Add CRS + sedConcTif.rio.write_crs("epsg:32118", inplace=True) + sedConcTif.rio.to_raster('C:/Users/arey/files/Projects/Newtown/SedFigs2/GeoTiff_RevB_' + + 'Sed_Conc4_Time_' + str(sedTstep) + '.tif') diff --git a/NTC_DFM/NTC_VelAnimation.py b/NTC_DFM/NTC_VelAnimation.py index 4efd0b6..f78983a 100644 --- a/NTC_DFM/NTC_VelAnimation.py +++ b/NTC_DFM/NTC_VelAnimation.py @@ -41,6 +41,8 @@ from pathlib import Path # Colormaps import cm_xml_to_matplotlib as cm +import pytz + #%% Read Model Log runLog = pd.read_excel('C:/Users/arey/files/Projects/Newtown/Model Log NTC.xlsx', 'ModelLog') @@ -48,303 +50,172 @@ runLog = pd.read_excel('C:/Users/arey/files/Projects/Newtown/Model Log NTC.xlsx' dataPath = "FlowFM/dflowfm/output/FlowFM_map.nc" histPath = "FlowFM/dflowfm/output/FlowFM_his.nc" -#%% Import Model results for static images - -# Define which models to import -modelPlot = [144] - -d3dfm_DataArray = [None] * (max(modelPlot)+1) -sedMass = [None] * (max(modelPlot)+1) # Available Mass of Sediment -erodeSed = [None] * (max(modelPlot)+1) # Erosion and deposition -initialBed = [None] * (max(modelPlot)+1) # Erosion and deposition -modelCRS = [None] * (max(modelPlot)+1) # Model CRS -sedConc = [None] * (max(modelPlot)+1) # Model CRS - -for i in modelPlot: - file_nc_map = Path(runLog['Run Location'][i]) / dataPath - - # Open Map file as xarrray Dataset - d3dfm_DataArray[i] = dfmt.open_partitioned_dataset(file_nc_map.as_posix()) - - # Get Var info - vars_pd = dfmt.get_ncvarproperties(d3dfm_DataArray[i]) - - # Get last timestep - tStep = d3dfm_DataArray[i].time[-1].values - - # Import sand sediment class - sedMass[i] = d3dfm_DataArray[i].mesh2d_bodsed[-1, :, 2].compute() - initialBed[i]= d3dfm_DataArray[i].mesh2d_s1[0, :] - d3dfm_DataArray[i].mesh2d_waterdepth[0, :] - finalBed = d3dfm_DataArray[i].mesh2d_s1[-1, :] - d3dfm_DataArray[i].mesh2d_waterdepth[-1, :] - - erodeSed[i] = finalBed - initialBed[i] - erodeSed[i] = erodeSed[i].compute() - initialBed[i] = initialBed[i].compute() - - # Import concentration of sediment 2 at timestep 95 and layer 10 - sedConc[i] = d3dfm_DataArray[i].mesh2d_sedfrac_concentration[94, 2, :, 9].compute() - - modelCRS[i] = vars_pd['EPSG_code']['projected_coordinate_system'] - -#%% Define axis limits -axLim_Xmin = [] -axLim_Xmax = [] - -axLim_Ymin = [] -axLim_Ymax = [] - -# Set axis limits -# Full domain -axLim_Xmin.append(298500) -axLim_Xmax.append(306800) - -axLim_Ymin.append(58500) -axLim_Ymax.append(68000) - -# Zoom North River -axLim_Xmin.append(303500) -axLim_Xmax.append(305500) - -axLim_Ymin.append(66000) -axLim_Ymax.append(68000) - -# Zoom Tribs -axLim_Xmin.append(305500) -axLim_Xmax.append(306800) - -axLim_Ymin.append(60000) -axLim_Ymax.append(62500) - -# Zoom far trib -axLim_Xmin.append(305800) -axLim_Xmax.append(305900) - -axLim_Ymin.append(60150) -axLim_Ymax.append(60400) - - -#%% Plot using DFM functions - -for i in modelPlot: - # Create figure - fig, axesFig = plt.subplots(1, 4, figsize=(8, 2)) - - axes = axesFig.flatten() - - # Loop through axis and plot with different limits - for subIDX in range(0, 4): - - # # Plot Sediment Mass - # smp = sedMass[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='turbo', - # vmin=0, vmax=2000, add_colorbar=False) - - # # Plot Erosion - # smp = erodeSed[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='coolwarm', - # vmin=-1, vmax=1, add_colorbar=False) - - # Plot Bed - smp = initialBed[i].ugrid.plot(ax=axes[subIDX], edgecolors='face', cmap='nipy_spectral', - vmin=-25, vmax=0, add_colorbar=False) - - - # Set axis labels - axes[subIDX].set_xlabel('') - axes[subIDX].set_ylabel('') - axes[subIDX].set_title('') - axes[subIDX].set_xticks([]) - axes[subIDX].set_yticks([]) - - axes[subIDX].set_xlim(axLim_Xmin[subIDX], axLim_Xmax[subIDX]) - axes[subIDX].set_ylim(axLim_Ymin[subIDX], axLim_Ymax[subIDX]) - - # Add basemap - # ctx.add_basemap(ax=axes, source=ctx.providers.Esri.WorldImagery, crs=modelCRS[i], attribution=False) - mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' - ctx.add_basemap(axes[subIDX], source=mapbox, crs=modelCRS[i]) - - - # cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Sand Mass [kg/m2]') - # cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Erosion/ deposition [m]') - cbar = plt.colorbar(smp, ax=axesFig.ravel().tolist(), label='Initial Bed [mNAVD88]') - - plt.show() - - # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/SedMass_' + runLog['Run'][i] + '.png', - # bbox_inches='tight', dpi=300) - - # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/ErodeDep_' + runLog['Run'][i] + '.png', - # bbox_inches='tight', dpi=300) - - # fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs2/InitialBed_' + runLog['Run'][i] + '.png', - # bbox_inches='tight', dpi=300) - -#%% Plot using DFM functions -#Cmap Path -cmap_path = 'C:/Users/arey/Repo/MATLAB_Q/Downloads/KeyColormaps/' - -for i in modelPlot: - fig, axes = plt.subplots(nrows=1, ncols=1, subplot_kw={'projection': ccrs.epsg(32118)}, figsize=(6, 6)) - fig.patch.set_facecolor('white') - - fig.tight_layout(pad=3.0) - - # Load cmaps - cmap = 'AJMR_Sed5_RevE.xml' - plotcmap = cm.make_cmap(cmap_path + cmap) - - pc = plot_netmapdata(ugrid_all[i].verts, values=modelMaxShear[i][0, :], - ax=axes, linewidth=0.5, cmap=plotcmap) - - axes.set_xlim(305900, 306400) - axes.set_ylim(61500, 62200) - - # axes.quiver(modelX[i], modelX[i], U[i], V[i], color='w', scale=1.00) - - mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' - ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118') - - extent = axes.get_extent() - axes.set_xticks(np.linspace(extent[0], extent[1], 4)) - axes.set_yticks(np.linspace(extent[2], extent[3], 4)) - - pc.set_clim([0, 0.145]) - cbarTicks = [0, 0.0378, 0.063, 0.0826, 0.11, 0.145] - cbar = fig.colorbar(pc, ax=axes, shrink=0.95, - ticks=cbarTicks) - cbar.set_label('Total Max Shear Stress [N/m^2]') - - # Add label on colorbar - cbar.ax.text(0.08, (cbarTicks[0]+cbarTicks[1])/2, 'CLAY', ha='center', va='center', - rotation=90, fontweight='bold') - cbar.ax.text(0.08, (cbarTicks[1]+cbarTicks[2])/2, 'FSILT', ha='center', va='center', - rotation=90, fontweight='bold') - cbar.ax.text(0.08, (cbarTicks[2]+cbarTicks[3])/2, 'MSILT', ha='center', va='center', - rotation=90, fontweight='bold') - cbar.ax.text(0.08, (cbarTicks[3]+cbarTicks[4])/2, 'CSILT', ha='center', va='center', - rotation=90, fontweight='bold') - cbar.ax.text(0.08, (cbarTicks[4]+cbarTicks[5])/2, 'FSAND', ha='center', va='center', - rotation=90, fontweight='bold') - - plt.show() - - fig.savefig('C:/Users/arey/files/Projects/Newtown/SedFigs/Shear_Class_' + runLog['Run'][i] + '.png', - bbox_inches='tight', dpi=300) - #%% Import Model results for animations # modelPlot = [20, 22] # modelPlot = [18, 20, 21, 22] # modelPlot = [18, 21] -modelPlot = [25, 26] +modelPlot = [146] -modelSedConc = [None] * (max(modelPlot)+1) -modelMaxShear_animate = [None] * (max(modelPlot)+1) -ugrid_all = [None] * (max(modelPlot)+1) +d3dfm_DataArray = [None] * (max(modelPlot)+1) +modelVelX = [None] * (max(modelPlot)+1) +modelVelY = [None] * (max(modelPlot)+1) tSteps = [None] * (max(modelPlot)+1) +# Time steps to import +tStepsImport = range(1000, 7000) + for i in modelPlot: + # Create variables + modelVelX[i] = [None] * (max(tStepsImport)+1) + modelVelY[i] = [None] * (max(tStepsImport)+1) + file_nc_map = Path(runLog['Run Location'][i]) / dataPath + d3dfm_DataArray[i] = dfmt.open_partitioned_dataset(file_nc_map.as_posix()) + # Get Var info - vars_pd, dims_pd = get_ncvardimlist(file_nc=file_nc_map.as_posix()) + vars_pd = dfmt.get_ncvarproperties(d3dfm_DataArray[i]) - # Import - tSteps[i] = get_timesfromnc(file_nc=file_nc_map.as_posix(), varname='time') + # Get time steps + tSteps[i] = d3dfm_DataArray[i].time.values - ugrid_all[i] = get_netdata(file_nc=file_nc_map.as_posix()) - modelSedConc[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_sedfrac_concentration', - timestep='all', station=0, layer=0) #timestep='all' OR timestep=range(0, 30) - modelMaxShear_animate[i] = get_ncmodeldata(file_nc=file_nc_map.as_posix(), varname='mesh2d_tausmax', - timestep='all') + # Loop through time steps and import data + for t in tStepsImport: + # Import 3D velocity data + modelVelXin = d3dfm_DataArray[i].mesh2d_ucx[t, :, :].compute() + modelVelYin = d3dfm_DataArray[i].mesh2d_ucy[t, :, :].compute() + # Average over layers + modelVelX[i][t] = modelVelXin.mean(axis=1) + modelVelY[i][t] = modelVelYin.mean(axis=1) -# %% Import water levels from hist + print(t) + + # %% Import water levels from hist modelHistWL = [None] * (max(modelPlot)+1) -modelHistTime = [None] * (max(modelPlot)+1) -# tSteps = [None] * (max(modelPlot)+1) for i in modelPlot: + # Hist file path file_nc_hist = Path(runLog['Run Location'][i]) / histPath - vars_pd, dims_pd = get_ncvardimlist(file_nc=file_nc_hist.as_posix()) + # Open hist file dataset + data_xr = xr.open_mfdataset(file_nc_hist, preprocess=dfmt.preprocess_hisnc) - # Get station names - histStations = get_hisstationlist(file_nc=file_nc_hist.as_posix(), varname='waterlevel') + # Get stations + stations_pd = data_xr['stations'].to_dataframe() - # Import - modelHistTime[i] = get_timesfromnc(file_nc=file_nc_hist.as_posix(), varname='time') - modelHistWL[i] = get_ncmodeldata(file_nc=file_nc_hist.as_posix(), varname='waterlevel', - timestep='all', station=1) + # Convert to Pandas + modelHistWL[i] = data_xr.waterlevel.sel(stations=['West_WL', 'East_WL']).to_pandas() + + # Put in UTC if needed + if i == 104 or i == 105 or i >= 135: + modelHistWL[i] = modelHistWL[i].tz_localize( + pytz.timezone('EST')).tz_convert(pytz.utc) + else: + modelHistWL[i] = modelHistWL[i].tz_localize( + pytz.timezone('UTC')) -#%% Sediment animations +#%% Velocity Quiver animations plt.rcParams['animation.ffmpeg_path'] = \ - 'C:/Users/arey/Local/ffmpeg-2022-02-14-git-59c647bcf3-full_build/bin/ffmpeg.exe' + 'C:/Users/arey/files/Local/ffmpeg-2022-02-14-git-59c647bcf3-full_build/bin/ffmpeg.exe' mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckex9vtri0o6619p55sl5qiyv/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' x, y, arrow_length = 0.93, 0.95, 0.12 # fontprops = fm.FontProperties(size=12) # Set model to plot -# modelPlot = [20, 22] -# modelPlot = [18, 20, 21, 22] -modelPlot = [25, 26] -cmap_path = 'C:/Users/arey/Repo/MATLAB_Q/Downloads/KeyColormaps/' +modelPlot = [84] +tStepsPlot = tStepsImport +tStepsPlot = range(970, 971) +tStepsPlot = range(955, 990) + + +vmin = 0 # Legend min +vmax = 0.05 # Legend max +scale = 30 # Size of arrows, smaller is bigger #25 +qmax = 0.002# Scaling threshold. Values above this are set to 1 +amin = 0.1 # Minimum arrow length. A dot is used below this value for m in modelPlot: - vmax = 0.75 - vmin = 0 - # vmax = 0.145 - # vmin = 0 - - cbarTicks = [0, 0.0378, 0.063, 0.0826, 0.11, 0.145] + # Find water level plot range + wlPlotMinIDX = wlHistTimeIDX = np.argmin(abs(tSteps[m][tStepsPlot[0]] - modelHistWL[m].index.values)) + wlPlotMaxIDX = wlHistTimeIDX = np.argmin(abs(tSteps[m][tStepsPlot[-1]] - modelHistWL[m].index.values)) cmapName = 'turbo' cmap = mpl.cm.turbo - # cmapName = 'AJMR_Sed5_RevE.xml' - # cmap = cm.make_cmap(cmap_path + cmapName) - - # Zoom limits + # # Zoom limits # xLimits = [305900, 306500] # yLimits = [61400, 62200] - # Wide Limits - xLimits = [302719.4, 306934.2] - yLimits = [60263.7, 64194.8] + # Corner limits + xLimits = [306088, 306444] + yLimits = [61370, 61741] # Setup Video - metadata = dict(title='NTC Sediment Animation', artist='Matplotlib', - comment='AJMR June 28, 2022') - writer = animation.FFMpegWriter(fps=2, metadata=metadata, codec='h264', bitrate=5000) - fig, axes = plt.subplots(figsize=(6, 6)) + metadata = dict(title='NTC Velocity', artist='Matplotlib', + comment='AJMR June 2, 2023') + writer = animation.FFMpegWriter(fps=1, metadata=metadata, codec='h264', bitrate=5000) + fig, axes = plt.subplots(figsize=(10, 10)) - writer.setup(fig, '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + - '06_Models/04_Delft3D/Figures/SedConc/HQ3_Wide_RevC_Long_Sed_Cond_' - + runLog['Run'][m] + '.mp4') + writer.setup(fig, 'C:/Users/arey/files/Projects/Newtown/Vel_Animation4.mp4') mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw' - for i in np.arange(1, len(modelSedConc[m])-1, 1): #289 + frameCount = 1 + + for t in tStepsPlot: #289 # Clear axes for video axes.cla() # Create new figure for stills - # fig, axes = plt.subplots(figsize=(6, 6)) + # fig, axes = plt.subplots(figsize=(12, 12)) # Clear inset WL Plot - # if i > 1: - # axin1.cla() + if frameCount > 1: + axin1.cla() - # Plot sediment - pc = plot_netmapdata(ugrid_all[m].verts, values=modelSedConc[m][i, 0, :, 0], - ax=axes, cmap=cmapName, - antialiaseds=False) - # Plot shear - # pc = plot_netmapdata(ugrid_all[m].verts, values=modelMaxShear_animate[m][i,:], - # ax=axes, cmap=cmap, - # antialiaseds=False) + # Cropt to zoom limits + croppedGridMask = (d3dfm_DataArray[m].mesh2d_face_x >= xLimits[0]) & \ + (d3dfm_DataArray[m].mesh2d_face_x <= xLimits[1]) & \ + (d3dfm_DataArray[m].mesh2d_face_y >= yLimits[0]) & \ + (d3dfm_DataArray[m].mesh2d_face_y <= yLimits[1]) + croppedGridMask = croppedGridMask.compute() + + # Rasterize Ugrid + rasterU = dfmt.rasterize_ugrid(d3dfm_DataArray[m]['mesh2d_ucx'].isel( + time=t, mesh2d_nLayers=4).where(croppedGridMask, drop=True).to_dataset(), resolution=10) + rasterV = dfmt.rasterize_ugrid(d3dfm_DataArray[m]['mesh2d_ucy'].isel( + time=t, mesh2d_nLayers=4).where(croppedGridMask, drop=True).to_dataset(), resolution=10) + + # Using fixed grid-> select desiered time step/layer/variable, convert to dataset, rasterize, then to numpy + overGrid_U_Scale = rasterU.to_array().to_numpy().squeeze() + overGrid_V_Scale = rasterV.to_array().to_numpy().squeeze() + + r = np.power(np.add(np.power(overGrid_U_Scale, 2), np.power(overGrid_V_Scale, 2)), 0.5) + + curPlotMask = np.sqrt((overGrid_U_Scale) ** 2 + (overGrid_V_Scale) ** 2) > qmax + + overGrid_U_Plot = np.zeros(overGrid_U_Scale.shape) + overGrid_V_Plot = np.zeros(overGrid_V_Scale.shape) + + # Normalize arrows to 1 where not zero + overGrid_U_Plot[r != 0] = overGrid_U_Scale[r != 0] / r[r != 0] + overGrid_V_Plot[r != 0] = overGrid_V_Scale[r != 0] / r[r != 0] + + # Scale arrows above qmax + overGrid_U_Plot[~curPlotMask] = overGrid_U_Plot[~curPlotMask] * (r[~curPlotMask] / qmax) + overGrid_V_Plot[~curPlotMask] = overGrid_V_Plot[~curPlotMask] * (r[~curPlotMask] / qmax) + + # Calculate velocity magnitude + velMag = np.sqrt((overGrid_U_Scale) ** 2 + (overGrid_V_Scale) ** 2) + + # Using fixed grid + qv = axes.quiver(rasterU['x'].values, rasterU['y'].values, overGrid_U_Plot, overGrid_V_Plot, velMag, scale=scale, + color='k', headlength=5, headwidth=5, headaxislength=5, width=0.002, minlength=amin) #200 # Set map limits axes.set_xlim(xLimits[0], xLimits[1]) @@ -353,33 +224,41 @@ for m in modelPlot: # Add basemap ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118') - # Set map ticks - axes.set_xticks(np.linspace(xLimits[0], xLimits[1], 4)) - axes.set_yticks(np.linspace(yLimits[0], yLimits[1], 4)) - # Color Limits - pc.set_clim([vmin, vmax]) + qv.set_clim([vmin, vmax]) + + # Set colormap for quiver + qv.set_cmap(cmap) # Add title - axes.title.set_text('Bottom Sediment Concentration at: ' + tSteps[m][i].strftime("%Y/%m/%d %H:%M")) - # axes.title.set_text('Total Max Shear Stress at: ' + str(tSteps[m][i])) + axes.title.set_text('Mid-depth velocity: ' + np.datetime_as_string(tSteps[m][t], unit='m')) # Add inset wl plot axin1 = axes.inset_axes([0.1, 0.1, 0.4, 0.15]) - axin1.plot(modelHistTime[m], modelHistWL[m]) - # Add current point - # Find closest time - abs_deltas_from_target_date = np.absolute(modelHistTime[m] - tSteps[m][i]) - axin1.scatter(modelHistTime[m][np.argmin(abs_deltas_from_target_date)], modelHistWL[m][np.argmin(abs_deltas_from_target_date)], 10, 'r') - # axin1.set_xticks([modelHistTime[m][0], modelHistTime[m][len(modelHistTime[m])-1]]) + # Find matching time step in model water level history + wlHistTimeIDX = np.argmin(abs(tSteps[m][t] - modelHistWL[m].index.values)) + + # Plot model water level history + axin1.plot(modelHistWL[m].index, modelHistWL[m]) + # Add current point + axin1.scatter(modelHistWL[m].index[wlHistTimeIDX], + modelHistWL[m]['East_WL'][wlHistTimeIDX], 10, 'r') + axin1.set_xticks([]) axin1.set_yticks([]) + # Set time axis limits + axin1.set_xlim([modelHistWL[m].index[wlPlotMinIDX], + modelHistWL[m].index[wlPlotMaxIDX]]) + # Set inset axis y limits + axin1.set_ylim([-0.5, 2]) + + plt.setp(axin1.get_xticklabels(), backgroundcolor="white") plt.setp(axin1.get_yticklabels(), backgroundcolor="white") - if i == 1: #i < 1000 + if frameCount == 1: #i < 1000 fig.subplots_adjust(right=0.8) norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) @@ -387,55 +266,31 @@ for m in modelPlot: # cax = plt.axes([0.82, 0.25, 0.04, 0.5]) cax = plt.axes([0.85, 0.15, 0.04, 0.7]) - # Add colorbar with designated tick marks - # cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap, - # norm=norm, - # orientation='vertical', - # ticks=cbarTicks) - # Add colorbar without designated tick marks cb1 = mpl.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm, orientation='vertical') # Colorbar lavel - cb1.set_label('Sediment Concentration [$kg/m^3$]') - # cb1.set_label('Total Max Shear Stress [N/m^2]') - - # Add label on colorbar - # cb1.ax.text(0.08, (cbarTicks[0] + cbarTicks[1]) / 2, 'CLAY', ha='center', va='center', - # rotation=90, fontweight='bold') - # cb1.ax.text(0.08, (cbarTicks[1] + cbarTicks[2]) / 2, 'FSILT', ha='center', va='center', - # rotation=90, fontweight='bold') - # cb1.ax.text(0.08, (cbarTicks[2] + cbarTicks[3]) / 2, 'MSILT', ha='center', va='center', - # rotation=90, fontweight='bold') - # cb1.ax.text(0.08, (cbarTicks[3] + cbarTicks[4]) / 2, 'CSILT', ha='center', va='center', - # rotation=90, fontweight='bold') - # cb1.ax.text(0.08, (cbarTicks[4] + cbarTicks[5]) / 2, 'FSAND', ha='center', va='center', - # rotation=90, fontweight='bold') - - - # Change label font size for only the primary axis - for item in ([axes.xaxis.label, axes.yaxis.label] + - axes.get_xticklabels() + axes.get_yticklabels()): - item.set_fontsize(4) - + cb1.set_label('Water velocity [m/s]') # Save video frame writer.grab_frame() + # plt.show() + print(t) plt.show() - print(i) + + # Increase frame count + frameCount += 1 # Save stills - fig.savefig('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + - '06_Models/04_Delft3D/Figures/SedConc/Wide_RevC_SedConc_' + - runLog['Run'][m] + '_Frame_' + str(i) + '.png', - bbox_inches='tight', dpi=300) # fig.savefig('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/' + - # '06_Models/04_Delft3D/Figures/ShearStress/Wide_ShearStress_' + + # '06_Models/04_Delft3D/Figures/SedConc/Wide_RevC_SedConc_' + # runLog['Run'][m] + '_Frame_' + str(i) + '.png', # bbox_inches='tight', dpi=300) + + # Save video writer.finish() diff --git a/NTC_DFM/Plot_Meas_versus_Model_Temp_Sal_EFDC_Delft3D_TM_2014-2015.py b/NTC_DFM/Plot_Meas_versus_Model_Temp_Sal_EFDC_Delft3D_TM_2014-2015.py index 030eec0..8d4256b 100644 --- a/NTC_DFM/Plot_Meas_versus_Model_Temp_Sal_EFDC_Delft3D_TM_2014-2015.py +++ b/NTC_DFM/Plot_Meas_versus_Model_Temp_Sal_EFDC_Delft3D_TM_2014-2015.py @@ -51,34 +51,19 @@ import datetime as datetime #%% Load in Measured data -dfIN=pd.read_csv('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/21_YSI Data/D1-D9_Combined/Combined_with_2015/All_Combine.csv') +YSI_df = pd.read_csv('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/21_YSI Data/D1-D9_Combined/Combined_with_2015/All_Combine.csv') # Filter out invalid data -df1 = dfIN -df1.loc[df1['Temp'] < -30, 'Temp'] = np.nan -df1.loc[df1['Sal'] < 0, 'Sal'] = np.nan -df1['DateTime'] = pd.to_datetime(df1['Date Time Local']) -df1.set_index('DateTime', inplace=True) +YSI_df.loc[YSI_df['Temp'] < -30, 'Temp'] = np.nan +YSI_df.loc[YSI_df['Sal'] < 0, 'Sal'] = np.nan +YSI_df['DateTime'] = pd.to_datetime(YSI_df['Date Time Local']) +YSI_df.set_index('DateTime', inplace=True) -df1.index = df1.index.tz_localize(pytz.timezone('America/New_York'), +# Set index as DateTime +YSI_df.index = YSI_df.index.tz_localize(pytz.timezone('America/New_York'), ambiguous='NaT', nonexistent='NaT') -df1.index = df1.index.tz_convert(pytz.utc) - -# Convert to Dataframe -YSI_df = pd.DataFrame(df1, columns=['Station', 'Temp', 'Sal', 'Depth']) -# Set index as DateTime - -#Change Station of interest (EB043, EK108, NC310, NC313, NC316, NC318) -Station=['EB043','EK108','NC310','NC313','NC316','NC318'] -x=4 -i=Station[x] - - -Station_interest=i -Top=str(Station_interest)+str('A') -Bot=str(Station_interest)+str('C') - +YSI_df.index = YSI_df.index.tz_convert(pytz.utc) #%% Load in EFDC Output FILE=r"//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/06_Models/02_EFDC/06_EFDC_Outputs/2014-2015/EFDC_Compiled_2014_2015.xlsx" @@ -135,16 +120,16 @@ NGG2_pd = NGG2_pd.tz_localize( # Drop ambiguous times NGG2_pd = NGG2_pd.dropna() - -#%% Load in Delft Results +#%% Read in model log runLog = pd.read_excel('C:/Users/arey/files/Projects/Newtown/Model Log NTC.xlsx', 'ModelLog') +#%% Load in Delft Results dataPath = "FlowFM_map.nc" histPath = "FlowFM_his.nc" #read in time series -modelPlot = [105]#, 106] +modelPlot = [105, 139, 141, 142, 143, 144]#, 106] Delft_WL = [None] * (max(modelPlot)+1) Delft_T_B = [None] * (max(modelPlot)+1) @@ -157,11 +142,6 @@ Delft_S_T = [None] * (max(modelPlot)+1) for i in modelPlot: file_nc_hist = pl.Path(runLog['Run Location'][i], 'FlowFM', 'dflowfm', 'output', 'FlowFM_his.nc') - # file_nc_hist = "C:/Users/mrobinson/OneDrive - W.F. Baird & Associates Coastal Engineers Ltd/Documents/11934.202 NTC/06 Models/06 Delft3DFM/Results/Run129/FlowFM_his.nc" - - - # Hist file path - #file_nc_hist = file_nc_hist.as_posix() # Open hist file dataset data_xr = xr.open_mfdataset(file_nc_hist, preprocess=dfmt.preprocess_hisnc) @@ -169,7 +149,7 @@ for i in modelPlot: # Get stations stations_pd = data_xr['stations'].to_dataframe() - # Get observations + # Get observation stations #observations_pd=data_xr['stations'].to_dataframe() # Convert to Pandas @@ -184,7 +164,7 @@ for i in modelPlot: Delft_S_T[i] = data_xr.salinity.sel(stations=['NC310', 'NC313','NC316','NC318','EB043','EK108'], laydim=9).to_pandas() # Put in UTC if needed - if i == 104 or i == 105: + if i == 104 or i == 105 or i >= 135: Delft_WL[i] = Delft_WL[i].tz_localize( pytz.timezone('EST')).tz_convert(pytz.utc) @@ -199,7 +179,7 @@ for i in modelPlot: #%% Plot Water Levels -i = 105 +i = 142 # Create Figure fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(18, 8), sharex=False, sharey=True) @@ -212,19 +192,21 @@ NGG2_pd['water_surface_elevation'].multiply(0.3048).plot(label='FFG', color='k', Delft_WL[i]['West_WL'].plot(label='Delft', color='r', linewidth=2) # Plot YSI -YSI_df.loc[YSI_df['Station'] == 'NC318A', 'Depth'].add(-1.3).plot( - label='YSI', color='b', linewidth=2, ax=ax) +# YSI_df.loc[YSI_df['Station'] == 'NC318A', 'Depth'].add(-1.3).plot( +# label='YSI', color='b', linewidth=2, ax=ax) # Plot EFDC -EFDC_dat['NC318A'].loc[:, 'WSE_m'].plot(label='EFDC', color='m', linewidth=2, ax=ax) +# EFDC_dat['NC318A'].loc[:, 'WSE_m'].plot(label='EFDC', color='m', linewidth=2, ax=ax) # Plot Telemac -TM_dat['NC318'].loc[:, 'WSE_m'].plot(label='EFDC', color='m', linewidth=2, ax=ax) +# TM_dat['NC318'].loc[:, 'WSE_m'].plot(label='EFDC', color='m', linewidth=2, ax=ax) # Set Time Bounds -dateMin = pd.to_datetime(datetime.datetime(2014, 11, 1, 0, 0, 0, 0), utc=True) -dateMax = pd.to_datetime(datetime.datetime(2014, 11, 13, 0, 0, 0, 0), utc=True) +# dateMin = pd.to_datetime(datetime.datetime(2014, 11, 1, 0, 0, 0, 0), utc=True) +# dateMax = pd.to_datetime(datetime.datetime(2014, 11, 13, 0, 0, 0, 0), utc=True) +dateMin = pd.to_datetime(datetime.datetime(2015, 7, 15, 0, 0, 0, 0), utc=True) +dateMax = pd.to_datetime(datetime.datetime(2015, 9, 1, 0, 0, 0, 0), utc=True) # Set X Axis ax.set_xlim(dateMin, dateMax) @@ -255,12 +237,16 @@ fig, ax = plt.subplots(nrows=2, ncols=6, figsize=(16, 7), sharex=True, sharey=Tr # fig.suptitle('Comparison of Predicted and Measured Temperature at NC316', fontsize=18) # Set Time Bounds for all axis since linked -# dateMin = pd.to_datetime(datetime.datetime(2015, 7, 15, 0, 0, 0, 0), utc=True) -# dateMax = pd.to_datetime(datetime.datetime(2015, 9, 1, 0, 0, 0, 0), utc=True) +dateMin = pd.to_datetime(datetime.datetime(2015, 7, 15, 0, 0, 0, 0), utc=True) +dateMax = pd.to_datetime(datetime.datetime(2015, 9, 1, 0, 0, 0, 0), utc=True) # dateMin = pd.to_datetime(datetime.datetime(2015, 7, 28, 0, 0, 0, 0), utc=True) # dateMax = pd.to_datetime(datetime.datetime(2015, 8, 7, 0, 0, 0, 0), utc=True) -dateMin = pd.to_datetime(datetime.datetime(2014, 7, 1, 0, 0, 0, 0), utc=True) -dateMax = pd.to_datetime(datetime.datetime(2015, 10, 1, 0, 0, 0, 0), utc=True) +# dateMin = pd.to_datetime(datetime.datetime(2014, 7, 1, 0, 0, 0, 0), utc=True) +# dateMax = pd.to_datetime(datetime.datetime(2015, 10, 1, 0, 0, 0, 0), utc=True) + +# Set up Delft3D Variables +DelftInterp = [None] * (max(modelPlot)+1) +RMSE_D = [None] * (max(modelPlot)+1) # Loop through stations and plot @@ -271,17 +257,22 @@ for StatIDX, Stat in enumerate(Stations): ObsDat = YSI_df.loc[YSI_df['Station'] == Stat + 'A', 'Temp'] ObsDat = ObsDat.loc[(ObsDat.index > dateMin) & (ObsDat.index <= dateMax)] - DelftInterp = resampleToObs(Delft_T_T[i][Stat], ObsDat, '5T') - EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'TEM10'].resample('10T').nearest(), ObsDat, '5T') - # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'TEM10'].resample('15T').nearest(), ObsDat, '5T') + # Loop through Delft Runs + for i in modelPlot: + DelftInterp[i] = resampleToObs(Delft_T_T[i][Stat], ObsDat, '5T') + RMSE_D[i] = syncModelRMSE(DelftInterp[i], ObsDat) - RMSE_D = syncModelRMSE(DelftInterp, ObsDat) + EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'TEM10'].resample('10T').nearest(), ObsDat, '5T') RMSE_E = syncModelRMSE(EFDCInterp, ObsDat) + + # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'TEM10'].resample('15T').nearest(), ObsDat, '5T') # RMSE_T = syncModelRMSE(TMInterp, ObsDat) # Plot Top Layer Temperature ObsDat.plot(label='YSI', color='k', linewidth=2, ax=ax[0, StatIDX]) - ax[0, StatIDX].plot(Delft_WL[i].index, Delft_T_T[i][Stat], linestyle='-', label='Delft3D: ' + str(round(RMSE_D, 2)) + degree_sign + 'C') + for i in modelPlot: + ax[0, StatIDX].plot(Delft_WL[i].index, Delft_T_T[i][Stat], linestyle='-', label='Delft3D_' + str(i) + ':' + str(round(RMSE_D[i], 2)) + degree_sign + 'C') + EFDC_dat[Stat].loc[:, 'TEM10'].plot(linewidth=1, ax=ax[0, StatIDX], label='EFDC: ' + str(round(RMSE_E, 2)) + degree_sign + 'C') # TM_dat[Stat].loc[:, 'TEM10'].plot(linewidth=1, ax=ax[0, StatIDX], label='Telemac: ' + str(round(RMSE_T, 2)) + degree_sign + 'C') @@ -293,16 +284,20 @@ for StatIDX, Stat in enumerate(Stations): ObsDat = YSI_df.loc[YSI_df['Station'] == Stat + 'C', 'Temp'] ObsDat = ObsDat.loc[(ObsDat.index > dateMin) & (ObsDat.index <= dateMax)] - DelftInterp = resampleToObs(Delft_T_B[i][Stat], ObsDat, '5T') - EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'TEM1'].resample('10T').nearest(), ObsDat, '5T') - # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'TEM1'].resample('15T').nearest(), ObsDat, '5T') + # Loop through Delft Runs + for i in modelPlot: + DelftInterp[i] = resampleToObs(Delft_T_B[i][Stat], ObsDat, '5T') + RMSE_D[i] = syncModelRMSE(DelftInterp[i], ObsDat) - RMSE_D = syncModelRMSE(DelftInterp, ObsDat) + EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'TEM1'].resample('10T').nearest(), ObsDat, '5T') RMSE_E = syncModelRMSE(EFDCInterp, ObsDat) + + # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'TEM1'].resample('15T').nearest(), ObsDat, '5T') # RMSE_T = syncModelRMSE(TMInterp, ObsDat) ObsDat.plot(label='YSI', color='k', linewidth=2, ax=ax[1, StatIDX]) - ax[1, StatIDX].plot(Delft_WL[i].index, Delft_T_B[i][Stat], linestyle='-', label='Delft3D: ' + str(round(RMSE_D, 2)) + degree_sign + 'C') + for i in modelPlot: + ax[1, StatIDX].plot(Delft_WL[i].index, Delft_T_B[i][Stat], linestyle='-', label='Delft3D_' + str(i) + ':' + str(round(RMSE_D[i], 2)) + degree_sign + 'C') EFDC_dat[Stat].loc[:, 'TEM1'].plot(linewidth=1, ax=ax[1, StatIDX], label='EFDC: ' + str(round(RMSE_E, 2)) + degree_sign + 'C') # TM_dat[Stat].loc[:, 'TEM1'].plot(linewidth=1, ax=ax[1, StatIDX], label='Telemac: ' + str(round(RMSE_T, 2)) + degree_sign + 'C') @@ -327,23 +322,26 @@ plt.tight_layout() plt.show() # Save figure -fig.savefig('C:/Users/arey/files/Projects/Newtown/TempSalFigs/Temperature_All.png', - bbox_inches='tight', dpi=300) +# fig.savefig('C:/Users/arey/files/Projects/Newtown/TempSalFigs/Temperature_All_SalSens.png', +# bbox_inches='tight', dpi=300) #%% Plot Salinity AJMR Stations=['NC310', 'NC313', 'NC316', 'NC318', 'EB043', 'EK108'] # Create Figure -fig, ax = plt.subplots(nrows=2, ncols=6, figsize=(16, 7), sharex=True, sharey=True) +fig, ax = plt.subplots(nrows=2, ncols=6, figsize=(16, 10), sharex=True, sharey=True) # Set Time Bounds for all axis since linked -# dateMin = pd.to_datetime(datetime.datetime(2015, 7, 15, 0, 0, 0, 0), utc=True) -# dateMax = pd.to_datetime(datetime.datetime(2015, 9, 1, 0, 0, 0, 0), utc=True) -dateMin = pd.to_datetime(datetime.datetime(2015, 7, 28, 0, 0, 0, 0), utc=True) -dateMax = pd.to_datetime(datetime.datetime(2015, 8, 7, 0, 0, 0, 0), utc=True) +dateMin = pd.to_datetime(datetime.datetime(2015, 7, 15, 0, 0, 0, 0), utc=True) +dateMax = pd.to_datetime(datetime.datetime(2015, 9, 1, 0, 0, 0, 0), utc=True) +# dateMin = pd.to_datetime(datetime.datetime(2015, 7, 28, 0, 0, 0, 0), utc=True) +# dateMax = pd.to_datetime(datetime.datetime(2015, 8, 7, 0, 0, 0, 0), utc=True) +# Set up Delft3D Variables +DelftInterp = [None] * (max(modelPlot)+1) +RMSE_D = [None] * (max(modelPlot)+1) # Loop through stations and plot for StatIDX, Stat in enumerate(Stations): @@ -353,17 +351,22 @@ for StatIDX, Stat in enumerate(Stations): ObsDat = YSI_df.loc[YSI_df['Station'] == Stat + 'A', 'Sal'] ObsDat = ObsDat.loc[(ObsDat.index > dateMin) & (ObsDat.index <= dateMax)] - DelftInterp = resampleToObs(Delft_S_T[i][Stat], ObsDat, '5T') - EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'SAL10'].resample('10T').nearest(), ObsDat, '5T') - # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'SAL10'].resample('15T').nearest(), ObsDat, '5T') + # Loop through Delft Runs + for i in modelPlot: + DelftInterp[i] = resampleToObs(Delft_S_T[i][Stat], ObsDat, '5T') + RMSE_D[i] = syncModelRMSE(DelftInterp[i], ObsDat) - RMSE_D = syncModelRMSE(DelftInterp, ObsDat) + + EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'SAL10'].resample('10T').nearest(), ObsDat, '5T') RMSE_E = syncModelRMSE(EFDCInterp, ObsDat) + # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'SAL10'].resample('15T').nearest(), ObsDat, '5T') # RMSE_T = syncModelRMSE(TMInterp, ObsDat) # Plot Top Layer Salinity ObsDat.plot(label='YSI', color='k', linewidth=2, ax=ax[0, StatIDX]) - ax[0, StatIDX].plot(Delft_S_T[i].index, Delft_S_T[i][Stat], linestyle='-', label='Delft3D: ' + str(round(RMSE_D, 2)) + ' psu') + for i in modelPlot: + ax[0, StatIDX].plot(Delft_S_T[i].index, Delft_S_T[i][Stat], linestyle='-', label='Delft3D_' + str(i) + ':' + str(round(RMSE_D[i], 2)) + ' psu') + EFDC_dat[Stat].loc[:, 'SAL10'].plot(linewidth=1, ax=ax[0, StatIDX], label='EFDC: ' + str(round(RMSE_E, 2)) + ' psu') # TM_dat[Stat].loc[:, 'SAL10'].plot(linewidth=1, ax=ax[0, StatIDX], label='Telemac: ' + str(round(RMSE_T, 2)) + ' psu') @@ -373,17 +376,19 @@ for StatIDX, Stat in enumerate(Stations): ObsDat = YSI_df.loc[YSI_df['Station'] == Stat + 'C', 'Sal'] ObsDat = ObsDat.loc[(ObsDat.index > dateMin) & (ObsDat.index <= dateMax)] - DelftInterp = resampleToObs(Delft_S_B[i][Stat], ObsDat, '5T') - EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'SAL1'].resample('10T').nearest(), ObsDat, '5T') - # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'SAL1'].resample('15T').nearest(), ObsDat, '5T') + for i in modelPlot: + DelftInterp[i] = resampleToObs(Delft_S_B[i][Stat], ObsDat, '5T') + RMSE_D[i] = syncModelRMSE(DelftInterp[i], ObsDat) - RMSE_D = syncModelRMSE(DelftInterp, ObsDat) + EFDCInterp = resampleToObs(EFDC_dat[Stat].loc[:, 'SAL1'].resample('10T').nearest(), ObsDat, '5T') RMSE_E = syncModelRMSE(EFDCInterp, ObsDat) + # TMInterp = resampleToObs(TM_dat[Stat].loc[:, 'SAL1'].resample('15T').nearest(), ObsDat, '5T') # RMSE_T = syncModelRMSE(TMInterp, ObsDat) # Plot Bottom Layer Salinity ObsDat.plot(label='YSI', color='k', linewidth=2, ax=ax[1, StatIDX]) - ax[1, StatIDX].plot(Delft_S_B[i].index, Delft_S_B[i][Stat], linestyle='-', label='Delft3D: ' + str(round(RMSE_D, 2)) + ' psu') + for i in modelPlot: + ax[1, StatIDX].plot(Delft_S_B[i].index, Delft_S_B[i][Stat], linestyle='-', label='Delft3D_' + str(i) + ':' + str(round(RMSE_D[i], 2)) + ' psu') EFDC_dat[Stat].loc[:, 'SAL1'].plot(linewidth=1, ax=ax[1, StatIDX], label='EFDC: ' + str(round(RMSE_E, 2)) + ' psu') # TM_dat[Stat].loc[:, 'SAL1'].plot(linewidth=1, ax=ax[1, StatIDX], label='Telemac: ' + str(round(RMSE_T, 2)) + 'psu') @@ -397,7 +402,7 @@ ax[0, 0].xaxis.set_major_formatter( mpl_dates.ConciseDateFormatter(ax[0, 0].xaxis.get_major_locator())) # Set Y Axis -ax[0, 0].set_ylim(12, 27) +ax[0, 0].set_ylim(4, 27) ax[0, 0].set_ylabel('Salinity Top (psu)') ax[1, 0].set_ylabel('Salinity Bottom (psu)') @@ -406,6 +411,6 @@ plt.tight_layout() plt.show() # Save figure -fig.savefig('C:/Users/arey/files/Projects/Newtown/TempSalFigs/Salinity_2015_Zoom.png', - bbox_inches='tight', dpi=300) +# fig.savefig('C:/Users/arey/files/Projects/Newtown/TempSalFigs/Salinity_2015_Zoom_SensTests.png', +# bbox_inches='tight', dpi=300) diff --git a/SatBathy/.idea/inspectionProfiles/Project_Default.xml b/SatBathy/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 0000000..7c85937 --- /dev/null +++ b/SatBathy/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,13 @@ + + + + \ No newline at end of file diff --git a/SatBathy/.idea/inspectionProfiles/profiles_settings.xml b/SatBathy/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/SatBathy/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/SatBathy/.idea/misc.xml b/SatBathy/.idea/misc.xml new file mode 100644 index 0000000..d1bb16c --- /dev/null +++ b/SatBathy/.idea/misc.xml @@ -0,0 +1,4 @@ + + + + \ No newline at end of file diff --git a/SatBathy/.idea/modules.xml b/SatBathy/.idea/modules.xml new file mode 100644 index 0000000..bfafffd --- /dev/null +++ b/SatBathy/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/SatBathy/.idea/vcs.xml b/SatBathy/.idea/vcs.xml new file mode 100644 index 0000000..6c0b863 --- /dev/null +++ b/SatBathy/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/SatBathy/.idea/workspace.xml b/SatBathy/.idea/workspace.xml new file mode 100644 index 0000000..69040f7 --- /dev/null +++ b/SatBathy/.idea/workspace.xml @@ -0,0 +1,62 @@ + + + + + + + + + + + + + + + + + + + + + + 1683656322651 + + + + + + + + + \ No newline at end of file