# -*- coding: utf-8 -*- """ @author: aforsythe Copied from "P:/11934.201 Newtown Creek TPP – Privileged and Confidential/05_Analyses/07 ADCP/NC_CurrentMeter_All_Phase1_all_data_2012_05_20/ADCP_Plot_v4.py" _AJMR """ import pandas as pd import matplotlib.pyplot as plt import matplotlib.dates as mdates import numpy as np import geopandas as gp gp.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw' import scipy as sp import scipy.ndimage from astropy.convolution import Gaussian2DKernel import contextily as ctx import os from shapely.geometry import Point #%% read in data dataPath = '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/05_Analyses/07 ADCP/NC_CurrentMeter_All_Phase1_all_data_2012_05_20/' df = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NC_CurrentMeter_All_Phase1_all_data_2013_07_16/NC_CurrentMeter_All_Phase1_all_mobile_2013_05_02.xlsx', sheet_name='mag_all_mobile') #convert data to geodataframe gdf = gp.GeoDataFrame(df, geometry=gp.points_from_xy(df.Longitude, df.Latitude, crs="EPSG:4326")) #convert data to CRS of D3D gdf.geometry = gdf.geometry.to_crs("EPSG:32118") # Add new cols to geodataframe gdf['Distance'] = np.zeros([len(gdf), 1]) gdf['DistanceSmth'] = np.zeros([len(gdf), 1]) gdf['dist'] = np.zeros([len(gdf), 1]) gdf['date'] = np.zeros([len(gdf), 1]) # Set Date gdf['date'] = gdf['year'].astype(str) + '-' + gdf['mon'].astype(str).str.zfill(2) + '-' + gdf['day'].astype( str).str.zfill(2) # Column names to distances data for plotting depths1 = gdf.columns[11:43] # ADCP sensor depth reading column names to list depths = np.array([float(d.replace('m', '')) for d in depths1]) * -1 #loop through all transects transectCount = 0 transects = gdf['transect'].unique() kernel = Gaussian2DKernel(x_stddev=0.25) for transect_id in transects[range(221, 241)]: # Select a given trasect tMask = (gdf['transect']==transect_id) # Remove rows without locations tMask[pd.isna(gdf['Latitude'])] = False # Find distance betwen points in m gdf.loc[tMask, 'Distance'] = gdf[tMask].distance(gdf[tMask].shift()) # Many of the points are recorded as being at the same location... # Where the distance is zero, set it to half of the following distance # Find zeros and iloc after zeros zeroDist = (gdf['Distance'] == 0) & tMask zeroDistShift = zeroDist.shift() zeroDistShift.iloc[0] = False distShift = gdf['Distance'].shift(-1) distShift.iloc[-1] = False gdf.loc[tMask, 'DistanceSmth'] = gdf['Distance'][tMask] # Set zeros to half of following gdf.loc[zeroDist, 'DistanceSmth'] = distShift[zeroDist]/2 # Set following to half gdf.loc[zeroDistShift, 'DistanceSmth'] = gdf['Distance']/2 # Set initial to zero gdf.loc[gdf[tMask].index[0], 'DistanceSmth'] = 0 # If final location is duplicate, set to previous value if gdf.loc[gdf[tMask].index[-1], 'Distance']==0: gdf.loc[gdf[tMask].index[-1], 'DistanceSmth'] = gdf.loc[gdf[tMask].index[-2], 'DistanceSmth'] # Get cumulative sum of distances gdf.loc[tMask, 'dist'] = gdf.loc[tMask, 'DistanceSmth'].cumsum() # get velocity data all rows columns 11-43 V = gdf.values[tMask, 11:43].astype(float) # #__________________________________________________________________________________________________________________________ # %% Plotting if transectCount == 0: fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(16, 8)) fig.tight_layout(pad=2.5) ax = axes.flat colormap = 'jet' vmin=0 vmax=0.5 VS = sp.ndimage.filters.gaussian_filter(V, [1, 1], mode='nearest') # vDatEnd = (~np.isnan(V)).cumsum(1).argmax(1).astype(int) + 1 # VS = convolve(V, kernel) # for i in range(0, len(VS)): # VS[i, vDatEnd[i]:-1] = np.nan pltDat = ax[transectCount].pcolormesh(gdf.loc[tMask, 'dist'], depths, np.transpose(VS), shading='auto', vmin=vmin, vmax=vmax, cmap=colormap) cbar = fig.colorbar(pltDat, ax=ax[transectCount], shrink=0.95,aspect=5) cbar.set_label('Magnitude [m/s]') ax[transectCount].set_xlabel('Distance Along Transect [m]') ax[transectCount].set_ylabel('Depth below WSL [m]') stationStart = next((i for i, j in enumerate(tMask) if j), None) ax[transectCount].set_title(gdf.loc[stationStart, 'station']) ax[transectCount].set_ylim(-6, 0) transectCount = transectCount + 1 plt.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Transects221-241.png', bbox_inches='tight', dpi=300) # %% Load in moored data df_moored_data = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NC_CurrentMeter_All_Phase1_all_data_2013_07_16/NC_CurrentMeter_All_Phase1_all_moored_2013_05_20.xlsx', sheet_name='mag_all_moored') # Shift col names to put second back in colIN = list(df_moored_data.columns) colOUT = colIN[0:6] + ['second'] + colIN[6:-1] df_moored_data.columns = colOUT ### Moored current meter locations from here: ### file://srv-ott3/Projects/11934.201%20Newtown%20Creek%20TPP%20%E2%80%93%20Privileged%20and%20Confidential/03_Data/01_BkgrdReports/NC_DRAFT_DSR_Submittal_No_3_2013-07-03.pdf ### Table 3-1, PDF page 65 ### Locations change between deployment 1-9 and 10/11 df_moored = pd.DataFrame( {'depOrig': ['NC083CM', 'NC081CM', 'NC082CM', 'NC086CM', 'EK023CM'], 'depCurr': ['NC086CM', 'NC087CM', 'NC088CM', 'NC089CM', 'EK023CM'], 'Northing': [208519.95, 206083.24, 203835.10, 201381.55, 200827.33], 'Easting': [996198.97, 1000959.45, 1004387.42, 1005283.07, 1004644.90], 'bedElev': [-16, -17, -20, -18, -20]}) gdf_moored_loc = gp.GeoDataFrame( df_moored, geometry=gp.points_from_xy(df_moored.Easting, df_moored.Northing), crs="EPSG:2263") gdf_moored_loc.geometry = gdf_moored_loc.geometry.to_crs("EPSG:32118") # Loop through Station IDs for deployments 1-9 and assign locations # for d in range(1,10): # depMask = df_moored_data['deployment'] == ('dep' + str(d)) # for stationIDX, station in enumerate(df_moored['depOrig']): # stationMask = (df_moored_data['station'] == station) & depMask # # # Assign geography to station plus deployment # df_moored_data.loc[stationMask, 'Northing'] = df_moored.loc[stationIDX, 'Northing'] # df_moored_data.loc[stationMask, 'Easting'] = df_moored.loc[stationIDX, 'Easting'] # Loop through Station IDs for deployments 10-11 and assign locations for d in range(1,12): depMask = df_moored_data['deployment'] == 'dep' + str(d) for stationIDX, station in enumerate(df_moored['depCurr']): stationMask = (df_moored_data['station'] == station) & depMask # Assign geography to station plus deployment df_moored_data.loc[stationMask, 'Northing'] = df_moored.loc[stationIDX, 'Northing'] df_moored_data.loc[stationMask, 'Easting'] = df_moored.loc[stationIDX, 'Easting'] # Create geodataframe and convert to USSP gdf_moored = gp.GeoDataFrame( df_moored_data, geometry=gp.points_from_xy(df_moored_data.Easting, df_moored_data.Northing), crs="EPSG:2263") #convert data to CRS of D3D gdf_moored.geometry = gdf_moored.geometry.to_crs("EPSG:32118") gdf_moored['date'] = pd.to_datetime(gdf_moored['year'].astype(str) + '-' + gdf_moored['month'].astype(str).str.zfill(2) + '-' + gdf_moored['day'].astype(str).str.zfill(2) + ' ' + gdf_moored['hour'].astype(str).str.zfill(2) + ':' + gdf_moored['minute'].astype(str).str.zfill(2)) # %% Plot moored data # Column names to distances data for plotting depths1 = gdf_moored.columns[8:51] # ADCP sensor depth reading column names to list depths_moored = np.array([float(d.replace('m', '')) for d in depths1]) colormap = 'jet' vmin = 0 vmax = 0.2 fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(16, 8)) fig.tight_layout(pad=2) # Loop through Station IDs for deployments 1-9 and assign locations for d in range(1,12): depMask = gdf_moored['deployment'] == ('dep' + str(d)) for stationIDX, station in enumerate(df_moored['depCurr']): stationMask = (gdf_moored['station'] == station) & depMask V = gdf_moored.values[stationMask, 8:51].astype(float) if len(gdf_moored.loc[stationMask, 'date']) != 0: pltDat = axes[stationIDX].pcolormesh(gdf_moored.loc[stationMask, 'date'], depths_moored, np.transpose(V), shading='auto', vmin=vmin, vmax=vmax, cmap=colormap) axes[stationIDX].set_ylim(0, 7) fmt_half_year = mdates.MonthLocator(interval=1) axes[stationIDX].xaxis.set_major_locator(fmt_half_year) axes[stationIDX].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) if d == 1: axes[stationIDX].set_ylabel('Elevation [m]') axes[stationIDX].set_title(station) cbar = fig.colorbar(pltDat, ax=axes[stationIDX], shrink=1.1, aspect=3) cbar.set_label('Magnitude [m/s]') pltDat.set_clim([0, 0.2]) plt.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/ADCP_2012.png', bbox_inches='tight', dpi=300) # %% Load in ADCP2 data # Read in locations df_adcp2_locs = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NCP2_ADCP_D1-D3/AQ_LocationsSO_ADCP_ADV_overview_Coords_20150126_AJMR.xlsx', sheet_name='ADCP') gdf_adcp2_locs = gp.GeoDataFrame(df_adcp2_locs, geometry=gp.points_from_xy(df_adcp2_locs.X_NYSPLI, df_adcp2_locs.Y_NYSPLI), crs="EPSG:2263") gdf_adcp2_locs = gdf_adcp2_locs.to_crs("EPSG:32118") adcp2_data_path = '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/' adcp2_paths = ['NCP2_ADCP_D1-D3/ADCP_D1_070314_090814', 'NCP2_ADCP_D1-D3/ADCP_D2_090914_110214', 'NCP2_ADCP_D1-D3/ADCP_D3_110414_010615', 'NCP2_ADCP_D4-D5/D4_010715_030315', 'NCP2_ADCP_D4-D5/D5_030415_050515'] adcp2_gdfs = dict() for depIDX, adcp2_path in enumerate(adcp2_paths): adcp2_files = os.listdir(adcp2_data_path + adcp2_path) # returns list of files in adv folder for adcp2_file in adcp2_files: if '.xls' in adcp2_file or '.csv' in adcp2_file: if '.xls' in adcp2_file: df_in = pd.read_excel(adcp2_data_path + adcp2_path + '/' + adcp2_file) else: df_in = pd.read_csv(adcp2_data_path + adcp2_path + '/' + adcp2_file) for stationIDX, station in enumerate(df_adcp2_locs['parent_loc_code']): advStrIDX = adcp2_file.find('_')+1 if adcp2_file[advStrIDX:advStrIDX+3] in station: adcp2_geo_x = np.ones([len(df_in), 1]) * df_adcp2_locs.X_NYSPLI[stationIDX] adcp2_geo_y = np.ones([len(df_in), 1]) * df_adcp2_locs.Y_NYSPLI[stationIDX] df_in['date'] = pd.to_datetime(df_in['Year'].astype(str) + '-' + df_in['Month'].astype( str).str.zfill(2) + '-' + df_in['Day'].astype( str).str.zfill(2) + ' ' + df_in['Hour'].astype( str).str.zfill(2) + ':' + df_in['Minute'].astype( str).str.zfill(2) + ':' + df_in['Second'].astype(str).str.zfill(2)) gdf_in = gp.GeoDataFrame( df_in, geometry=gp.points_from_xy(adcp2_geo_x, adcp2_geo_y), crs="EPSG:2263") if adcp2_file[advStrIDX:advStrIDX + 3] not in adcp2_gdfs: adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]] = dict() if adcp2_file[0:advStrIDX-1] not in adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]]: adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]] = dict() if depIDX+1 not in adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]]: adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]][depIDX+1] = gdf_in # %% Plot ADCP2 Data fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(9, 8)) fig.tight_layout(pad=2) colormap = 'jet' vmin = 0 vmax = 0.2 for stationIDX, stat in enumerate(adcp2_gdfs): for depIDX, depdat in enumerate(adcp2_gdfs[stat]['mag']): depths1 = adcp2_gdfs[stat]['mag'][depdat].columns[6:-2] # ADCP sensor depth reading column names to list depths_moored = np.array([float(d.replace('m', '')) for d in depths1]) V = adcp2_gdfs[stat]['mag'][depdat].values[:, 6:-2].astype(float) pltDat = axes[stationIDX].pcolormesh(adcp2_gdfs[stat]['mag'][depdat].loc[:, 'date'], depths_moored, np.transpose(V), shading='auto', vmin=vmin, vmax=vmax, cmap=colormap) axes[stationIDX].set_ylim(0, 7) axes[stationIDX].set_xlim(pd.to_datetime("2014-07-01"), pd.to_datetime('2015-05-15')) #axes[stationIDX, depIDX].format_xdata = mdates.DateFormatter('%Y-%m') fmt_half_year = mdates.MonthLocator(interval=1) axes[stationIDX].xaxis.set_major_locator(fmt_half_year) axes[stationIDX].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) axes[stationIDX].set_title(str(stat)) axes[stationIDX].set_ylabel('Elevation [m]') cbar = fig.colorbar(pltDat, ax=axes[stationIDX], shrink=1.1, aspect=3) cbar.set_label('Magnitude [m/s]') pltDat.set_clim([0, 0.2]) plt.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/ADCP_2014.png', bbox_inches='tight', dpi=300) # %% Load in Water Level Data # Gauge Locations from map gdf_gaugeLoc = gp.read_file('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/NTC6.kml') gdf_gaugeLocUSSP = gdf_gaugeLoc.to_crs("EPSG:32118") # All in Feet NAVD88 df_wl_FFG = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP1_Gauge_Elev_Data_20130521/NC_Gauge_Elev_Data_20130521.xlsx', sheet_name='Field_Facility_Gauge') gdf_wl_FFG = gp.GeoDataFrame(df_wl_FFG, geometry=gp.points_from_xy(np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].x[2], np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].y[2]), crs="EPSG:4326") gdf_wl_FFG.geometry = gdf_wl_FFG.geometry.to_crs("EPSG:32118") df_wl_NGG1 = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP1_Gauge_Elev_Data_20130521/NC_Gauge_Elev_Data_20130521.xlsx', sheet_name='National_Grid_Gauge') gdf_wl_NGG1 = gp.GeoDataFrame(df_wl_NGG1, geometry=gp.points_from_xy(np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].x[3], np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs="EPSG:4326") gdf_wl_NGG1.geometry = gdf_wl_NGG1.geometry.to_crs("EPSG:32118") # Also includes temperature df_wl_NGG2 = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP2_NG_TideGauge_20160113/NCP2_NG_TideGauge_20160113.xlsx', sheet_name='Sheet1') gdf_wl_NGG2 = gp.GeoDataFrame(df_wl_NGG2, geometry=gp.points_from_xy(np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].x[3], np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs="EPSG:4326") gdf_wl_NGG2.geometry = gdf_wl_NGG2.geometry.to_crs("EPSG:32118") # %% Plot Water Level Data fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4)) axes.plot(gdf_wl_FFG.Date_time, gdf_wl_FFG['Water Surface Elevation (ft)']*0.3048, label='Field Facility Gauge') axes.plot(gdf_wl_NGG1.Date_time, gdf_wl_NGG1['Water Surface Elevation (ft)']*0.3048, label='National Grid Gauge Deployment 1') axes.plot(gdf_wl_NGG2.datetime, gdf_wl_NGG2['water_surface_elevation']*0.3048, label='National Grid Gauge Deployment 2') fmt_half_year = mdates.MonthLocator(interval=6) axes.xaxis.set_major_locator(fmt_half_year) axes.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) axes.set_ylabel('Water Surface Elevation [mNAVD88]') axes.set_title('Water Surface Elevation') axes.legend() fig.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/WaterLevel.png', bbox_inches='tight', dpi=300) # %% Load temperature data df_tdat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/10_Salinity/tData.xlsx', skiprows=[1]) # Day Month order is reversed df_tdat['date'] = pd.to_datetime(df_tdat['CollectionDate']) df_tsample = pd.read_csv('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/10_Salinity/tSampleLocation.csv') # Create geodataframe and convert to USSP gdf_tsample = gp.GeoDataFrame( df_tsample, geometry=gp.points_from_xy(df_tsample.EastCoordinate, df_tsample.NorthCoordinate), crs="EPSG:2263") #convert data to CRS of D3D gdf_tsample.geometry = gdf_tsample.geometry.to_crs("EPSG:32118") tdat_geo = list() tdat_facility = list() for i in range(0, len(df_tdat)): #sampleID = df_tdat['LocationID'][i][0:df_tdat['LocationID'][i].find('_')] #geoFind = df_tsample.loc[df_tsample['ParentLocationID'] == sampleID] geoFind = df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].geometry.values if len(geoFind) !=0: tdat_geo.append(df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].geometry.values[0]) tdat_facility.append(df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].SourceArea.values[0]) else: tdat_geo.append(Point()) tdat_facility.append('') # Create geodataframe gdf_tdat = gp.GeoDataFrame(df_tdat, geometry=tdat_geo, crs="EPSG:32118") gdf_tdat.loc[:, 'SourceArea'] = tdat_facility gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='ft', 'DepthM'] = gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='ft', 'BeginDepth']*0.3048 gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='cm', 'DepthM'] = gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='cm', 'BeginDepth']*0.01 # Import spring observations df_springDat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/11_Temperature/NC_Spr2012_Benthic_Water_Quality_FieldData&Observations_20121022.xlsx') # Create geodataframe and convert to USSP gdf_springDat = gp.GeoDataFrame( df_springDat, geometry=gp.points_from_xy(df_springDat.x_coord_as_numeric, df_springDat.y_coord_as_numeric), crs="EPSG:2263") gdf_springDat.geometry = gdf_springDat.geometry.to_crs("EPSG:32118") # Import Summer observations df_summerDat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/11_Temperature/NC_Sum2012_Benthic_Water_Quality_FieldData&Observations_20130205.xlsx') # Create geodataframe and convert to USSP gdf_summerDat = gp.GeoDataFrame( df_summerDat, geometry=gp.points_from_xy(df_summerDat.x_coord_as_numeric, df_summerDat.y_coord_as_numeric), crs="EPSG:2263") gdf_summerDat.geometry = gdf_summerDat.geometry.to_crs("EPSG:32118") # Merged df_SpringSummerDat = pd.concat([df_springDat, gdf_summerDat], ignore_index=True) gdf_SpringSummerDat = gp.GeoDataFrame( df_SpringSummerDat, geometry=gp.points_from_xy(df_SpringSummerDat.x_coord_as_numeric, df_SpringSummerDat.y_coord_as_numeric), crs="EPSG:2263") gdf_SpringSummerDat.geometry = gdf_SpringSummerDat.geometry.to_crs("EPSG:32118") # %% Plot Salinity Time series fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 8)) plotMaskStation = (gdf_tdat.SourceArea == 'Newtown Creek') | (gdf_tdat.SourceArea == 'East Branch of Newtown Creek') plotMaskWater = (gdf_tdat.SampleMedium == 'Surface Water') plotMask = plotMaskStation & plotMaskWater pltDat = axes.scatter(gdf_tdat.loc[plotMask, 'date'], gdf_tdat.loc[plotMask, 'DepthM'], 10, gdf_tdat.loc[plotMask, 'NumericResult']) axes.set_ylim(10, 0) axes.set_ylabel('Depth below water surface [m]') axes.set_title('Surface Water Salinity Samples') axes.set_xlabel('Date') fmt_half_year = mdates.MonthLocator(interval=12) axes.xaxis.set_major_locator(fmt_half_year) axes.xaxis.set_major_formatter(mdates.DateFormatter('%Y')) cbar = fig.colorbar(pltDat, ax=axes, shrink=0.95) cbar.set_label('Salinity [PSU]') pltDat.set_clim([5, 40]) fig.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Salinity.png', bbox_inches='tight', dpi=300) # %% Plot Temperature Time series ## Additional salinity obs are here! fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(6, 8)) fig.tight_layout(pad=2) for i, miles in enumerate(np.arange(0, 3, 0.5)): plotMaskDistance = (gdf_SpringSummerDat.miles_from_NC_mouth > miles) & ( gdf_SpringSummerDat.miles_from_NC_mouth < miles + 0.5) plotMaskStation = (gdf_SpringSummerDat.subfacility_code == 'Newtown Creek') plotMaskVar = (gdf_SpringSummerDat.chemical_name == 'Temperature (field)')# Temperature (field)' plotMask = plotMaskDistance & plotMaskStation & plotMaskVar pltDat = axes[i].scatter(gdf_SpringSummerDat.loc[plotMask, 'location_start_date'], gdf_SpringSummerDat.loc[plotMask, 'elev']*0.3048, 10, gdf_SpringSummerDat.loc[plotMask, 'result_value']) axes[i].set_ylim(-10, 0) axes[i].set_ylabel('Elevation [m]') cbar = fig.colorbar(pltDat, ax=axes[i], shrink=0.95) cbar.set_label('Temp [degC]') pltDat.set_clim([5, 30]) fmt_half_year = mdates.MonthLocator(interval=1) axes[i].xaxis.set_major_locator(fmt_half_year) axes[i].xaxis.set_major_formatter(mdates.DateFormatter('%m-%Y')) axes[0].set_title('Surface Water Temperature Samples') axes[i].set_xlabel('Date') fig.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Temperature.png', bbox_inches='tight', dpi=300) # %% Plot Map 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=(8, 8)) axes.set_xlim(303500, 306500) axes.set_ylim(61000, 63750) gdf.plot(ax=axes, markersize=10, color='blue', label='Mobile ADCP') gdf_SpringSummerDat.plot(ax=axes, markersize=12, color='magenta', label='Temperature & Salinity') gdf_moored_loc.plot(ax=axes, markersize=20, color='red', label='Moored ADCP 2012') gdf_adcp2_locs.plot(ax=axes, markersize=20, color='orange', label='Moored ADCP 2014') gdf_gaugeLocUSSP.loc[2:3, 'geometry'].plot(ax=axes, markersize=20, color='green', label='Water Level Gauge') ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118') axes.set_xlabel('New York State Plane Easting [m]') axes.set_ylabel('New York State Plane Northing [m]') axes.legend() # axes[1].set_xlim(303500, 306500) # axes[1].set_ylim(61000, 63750) # gdf_SpringSummerDat.plot(ax=axes[1], markersize=12, color='magenta', label='Temperature & Salinity') # axes[1].set_xlabel('New York State Plane Easting [m]') # axes[1].legend() # ctx.add_basemap(axes[1], source=mapbox, crs='EPSG:32118') fig.show() fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/DataMap.png', bbox_inches='tight', dpi=300)