import os.path import numpy as np import math import csv from datetime import datetime, timedelta from pytz import timezone import pytz import time #strptime for e-mail alarms import json #for loading json file manually import requests # http://docs.python-requests.org/en/master/ import pandas as pd #for reading master spreadsheet xls file from zeep import Client import binascii import struct import gsw from scipy.interpolate import interp1d import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec from pathlib import Path wds_token = "eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJ1bmlxdWVfbmFtZSI6IjEzMDY1Iiwib3JnSWQiOiIyNCIsImp0aSI6IjY2ODZhMzYyMDc0NDQ2ZGFhNmYzM2QwZjk3NjEzYzg4IiwibmJmIjoxNzY3MjA3NDIyLCJleHAiOjE3NzQ5ODM0MjIsImlhdCI6MTc2NzIwNzQyMn0.NE36UJp42Td0XFPclEmAJBAPbVK5k0i0_BWGBATECrc" # alarm thresholds alarmEquilPressThresh = 87 #kPa #----------------------------------------------------# # Define functions def convertToDict(zeepOutput, result): jsonDict = dict() if zeepOutput.__values__.get("Error") == False: jsonString = zeepOutput.__values__.get(result) jsonDict = json.loads(jsonString) return jsonDict #define fn to convert datetime object to unix time stamp def dt2unix(dt): dt = dt.replace(tzinfo=pytz.utc) #make aware #ep = datetime.utcfromtimestamp(0) ep = datetime.fromtimestamp(0, timezone('UTC')) ep = ep.replace(tzinfo=pytz.utc) #make aware return (dt - ep).total_seconds() * 1000.0 #define fn to create and/or append to csv file def writeFile(myDirWeb, tBeg,myTitle, myParamNames, myTime, myParams, rnds): # global alarmCTDpressMsg, alarmEquilPressMsg, alarmTagIDmsg, tBeg, VR2CtagReading myPath = myDirWeb + myTitle +'.csv' myTimeOld = list() #time stamps for existing file if os.path.isfile(myPath): #if exists, read-in existing time stamps first f = open(myPath, 'r+') reader = csv.reader(f) f.readline() #skip header line for row in reader: myTimeOld.append(int(row[0])) else: #if doesn't exist, create new file & write header line f = open(myPath,'w') #write header line f.write('Time') for j in range(0,len(myParams)): f.write(',' +myParamNames[j]) #write empty line with tBeg (work-around for Highstock issue) f.write('\n') f.write(str(int(dt2unix(tBeg)) - 60*60*1000)) #**subtract an hour per CO2 issue** for j in range(0,len(myParams)): f.write(',') #write data for i in range(0, len(myTime)): if myTime[i] not in myTimeOld: #write only unique time stamps f.write('\n') f.write(str(myTime[i])) for j in range(0, len(myParams)): if np.isnan(myParams[j][i]): f.write( ',') else: f.write( ',' +str(round(myParams[j][i], rnds[j])) ) #check for alarms # Equil # if myTitle=='CO2equil': # if myParams[2][i] < alarmEquilPressThresh: #Equil-ON (3=OFF) # alarmEquilPressMsg = ( alarmEquilPressMsg + # '\n' +time.strftime("%m/%d/%Y %H:%M",time.gmtime(myTime[i]/1000)) + # ', ' +str(round(myParams[2][i], rnds[2])) +' kPa' ) # # tagID # elif myTitle=='tagID': # if len(myParams[0]) > 0: # alarmTagIDmsg = ( alarmTagIDmsg + # '\n' +time.strftime("%m/%d/%Y %H:%M",time.gmtime(myTime[i]/1000)) + # ', ' +str(myParams[0][i]) ) # alarmTagIDmsg = ( alarmTagIDmsg +'\n' +VR2CtagReading[i] ) f.close() #Wave Energy write function (mix of strings and numbers) def writeFileWave(myDirWeb, tBeg, myTitle, myParamNames, myTime, myParams, rnds): #global alarmCTDpressMsg, alarmEquilPressMsg, alarmTagIDmsg, tBeg, VR2CtagReading myPath = myDirWeb +'csv/' +myTitle +'.csv' myTimeOld = list() #time stamps for existing file if os.path.isfile(myPath): #if exists, read-in existing time stamps first f = open(myPath, 'r+') reader = csv.reader(f) f.readline() #skip header line for row in reader: myTimeOld.append(int(row[0])) else: #if doesn't exist, create new file & write header line f = open(myPath,'w') #write header line f.write('Time') for j in range(0,len(myParams)): f.write(',' +myParamNames[j]) #write empty line with tBeg (work-around for Highstock issue) f.write('\n') f.write(str(int(dt2unix(tBeg)) - 60*60*1000)) #**subtract an hour per CO2 issue** for j in range(0,len(myParams)): f.write(',') #write data for i in range(0, len(myTime)): if myTime[i] not in myTimeOld: #write only unique time stamps f.write('\n') f.write(str(myTime[i])) for j in range(0, len(myParams)): if isinstance(myParams[j][i],str): if myParams[j][i] =='SWELL': f.write( ',' + str(1)) elif myParams[j][i]=='AVERAGE' : f.write( ',' + str(2)) elif myParams[j][i] =='STEEP': f.write( ',' + str(3)) else: f.write (',' + str(4)) elif np.isnan(myParams[j][i]): f.write( ',') else: f.write( ',' +str(round(myParams[j][i], rnds[j])) ) #check for alarms # Equil # if myTitle=='CO2equil': # if myParams[2][i] < alarmEquilPressThresh: #Equil-ON (3=OFF) # alarmEquilPressMsg = ( alarmEquilPressMsg + # '\n' +time.strftime("%m/%d/%Y %H:%M",time.gmtime(myTime[i]/1000)) + # ', ' +str(round(myParams[2][i], rnds[2])) +' kPa' ) # # tagID # elif myTitle=='tagID': # if len(myParams[0]) > 0: # alarmTagIDmsg = ( alarmTagIDmsg + # '\n' +time.strftime("%m/%d/%Y %H:%M",time.gmtime(myTime[i]/1000)) + # ', ' +str(myParams[0][i]) ) # alarmTagIDmsg = ( alarmTagIDmsg +'\n' +VR2CtagReading[i] ) f.close() def writeFileXY(myDirWeb,myTitle, myParamNames, myParams, rnds): myPath = myDirWeb +'csv/' +myTitle +'.csv' #create new file & write header line f = open(myPath,'w') #write header line f.write(myParamNames[0] +',' +myParamNames[1]) #write data for i in range(0, len(myParams[0])): if (not np.isnan(myParams[0][i])) and (not np.isnan(myParams[1][i])): f.write('\n') f.write( str(round(myParams[0][i], rnds[0])) +',' +str(round(myParams[1][i], rnds[1])) ) f.close() #define fn to convert raw SBE o2 frequency to umol/kg #Soc, Foffset, Tau20, A, B, C, E def convO2(o2Hz, o2sol, O2cal, p, T, sigt): if o2Hz == np.nan or o2sol == np.nan or O2cal == np.nan or p == np.nan or T == np.nan or sigt == np.nan: return np.nan o2mll = ( O2cal[0] * (o2Hz + O2cal[1]) * (1.0 + O2cal[3]*T + O2cal[4]*T**2 + O2cal[5]*T**3) * o2sol * math.exp(O2cal[6] * p / (T+273.15)) ) o2umolkg = (o2mll * 1.4276) * (1e6/((sigt+1000)*32)) return o2umolkg #define fn to convert Eco-Puck counts to 'beta' (or chl) def convEP(counts, EPcalN): EPout = EPcalN[0] * (counts - EPcalN[1]) return EPout #define fn to convert Eco-Puck 'beta' to 'bb' def counts2bb(beta, wavel, S): #get beta and bb from pure sw (Morel, 1974) th = 117*math.pi/180; #all WetLabs sensors measure at 117 degrees d = 0.09; term1 = 1.38 * ( (wavel/500)**-4.32 ) term3 = 1 + ( (math.cos(th)**2) * (1-d)/(1+d) ) beta_w = term1 * ( (1+(0.3*S/37)) * 10**-4 ) * term3 bb_w = 8 * math.pi/3 * beta_w * (2+d)/(1+d) #calc betap, bbp, and estimate bb (ONLY KEEPING BB FOR NOW, ADD BETAP, BBP LATER) X = 1.1 betap = beta - beta_w bbp = 2*math.pi*X*betap bb = bbp + bb_w return bb #define fn to convert pH sensor voltage 'V1' to pH units # input T should be from CTD (for the WG or moorings) def pHV12pH(pHV1, pHcal, T): R = 8.31451 #gas constant (J mol-1 K-1) F = 96487 #converts between Coulombs and Faradays (1 F = 96487 C) ST = ( R * (T+273.15) * math.log(10) ) / F E0T = pHcal[0]-0.001 * ( (T+273.15)-pHcal[1] ) pH = (pHV1-E0T)/ST return pH #define fn to apply CO2 temperature-cal # input T should be from CO2 system, not CTD # standard correction is done now too def convCO2(xCO2meas, xCO2zeroMeas, CO2cal, T, xCO2stndMeas, xCO2stndKn): xCO2 = xCO2meas - xCO2zeroMeas xCO2conv = CO2cal[0] + CO2cal[1]*T + CO2cal[2]*xCO2 + CO2cal[3]*T*xCO2 xCO2convStnd = xCO2conv * xCO2stndKn / xCO2stndMeas return xCO2convStnd # #define fn to apply NEW CO2 offset- and temperature-cal # # new version also does standard correction in real-time (e.g., now) # def convCO2_v2(xCO2meas, xCO2zeroMeas, xCO2stndKn, xCO2stndMeas, CO2cal, T): # yOff = ( xCO2meas - ( CO2cal[0]*xCO2meas**2 + CO2cal[1]*xCO2meas + CO2cal[2] + #cal-derived # ( (xCO2stndMeas - xCO2stndKn - xCO2zeroMeas)/xCO2stndKn ) * xCO2meas + #field-nudged # xCO2zeroMeas ) ) # yTemp = ( (CO2cal[3]*yOff**2 + CO2cal[4]*yOff + CO2cal[5]) * T**2 + # (CO2cal[6]*yOff**2 + CO2cal[7]*yOff + CO2cal[8]) * T + # (CO2cal[9]*yOff**2 + CO2cal[10]*yOff + CO2cal[11]) ) # yFinal = yOff + (yOff - yTemp) # return yFinal #define fn to convert sea water xCO2 to pCO2 # input T, S here should be from CTD def calcpCO2(xCO2conv, T, S): TempK = T + 273.15 VPWP = math.exp(24.4543 - 67.4509 * (100/TempK) - 4.8489 * math.log(TempK/100)) VPCorrWP = math.exp(-0.000544 * S) VPSWWP = VPWP * VPCorrWP VPFac = 1 - VPSWWP pCO2 = xCO2conv * VPFac return pCO2 #define fn to write lonlat.txt in geoJSON fmt def writeLonLat(myDirWeb): myPathPos = myDirWeb +'csv/Position.csv' #Position.csv if os.path.isfile(myPathPos): #position.csv must exist! myPathLL = myDirWeb +'log/lonlat.txt' myPathRecentLL = myDirWeb +'log/recentlonlat.txt' fLL = open(myPathLL, 'w') fRecentLL = open(myPathRecentLL, 'w') fPos = open(myPathPos, 'r') #write json-style header fLL.write('{\n') fLL.write('\t"type": "FeatureCollection",\n') fLL.write('\t"features": [{\n') fLL.write('\t\t"type": "Feature",\n') fLL.write('\t\t"geometry": {\n') fLL.write('\t\t\t"type": "LineString",\n') fLL.write('\t\t\t"coordinates": [\n') #now, read-in rows from position.csv, write to lonlat.txt reader = csv.reader(fPos) #Can't have a trailing comma, get numRows first numRows = 0 for row in reader: numRows = numRows + 1 fPos.seek(0) #back to beginning fPos.readline() #skip header line fPos.readline() #skip empty line w/ tBeg for i in range(0,numRows-3): line = fPos.readline().replace('\n','').split(',') if line[1] != '' and line[2] != '' : fLL.write('\t\t\t\t[' +line[1] +', ' +line[2] +'],\n') line = fPos.readline().replace('\n','').split(',') fLL.write('\t\t\t\t[' +line[1] +', ' +line[2] +']\n') #write closing braces fLL.write('\t\t\t]\n') fLL.write('\t\t}\n') fLL.write('\t}]\n') fLL.write('}') fLL.close() fRecentLL.write('{\n') fRecentLL.write('\t"type": "Point",\n') fRecentLL.write('\t"coordinates": ') fRecentLL.write('[' +line[1] +', ' +line[2] +']\n') #write closing braces fRecentLL.write('}') fRecentLL.close() fPos.close() #define fn to write .json file for table def writeJSONtable(myPath, myList): with open(myPath,'w') as f: f.write('[\n') for i in range(0, len(myList) ): myRow = '\t{' myRow += '"id":' +str(i+1) +', ' myRow += '"yr":"' +myList[i][0] +'", ' myRow += '"durDay":' +str(myList[i][1]) +', ' myRow += '"dist":' +str(myList[i][2]) +'}' if i < len(myList)-1: myRow += ',' f.write(myRow +'\n') f.write(']') f.close() #define fn to burst-average T, S, O2 def burstAvgCTD(myCTDtime, myCTDparam): #convert lists to numpy arrays (eventually will move entirely to np arrays) myCTDtime = np.array(myCTDtime) myCTDparam = np.array(myCTDparam) #index bursts as diff's in time > 10 min myDiff = np.subtract(myCTDtime[1:], myCTDtime[0:-1]) f = np.concatenate( (np.array([-1]), np.where(myDiff >= 600000)[0], np.array([len(myCTDtime)-1])) ) #pre-allocate and fill averages myCTDtimeAvg = np.zeros(len(f)-1, dtype='int64') myCTDparamAvg = np.zeros(len(f)-1) for i in range(0, len(f)-1): ind = range( f[i]+1, f[i+1]+1 ) #have to add 1 to final index (python indexing) if len(ind) == 8 or len(ind) == 4: myCTDparamAvg[i] = np.mean(myCTDparam[ind]) myCTDtimeAvg[i] = np.mean(myCTDtime[ind]) zeroInd = np.where(myCTDtimeAvg==0) myCTDtimeAvg = np.delete(myCTDtimeAvg,zeroInd) myCTDparamAvg = np.delete(myCTDparamAvg,zeroInd) return myCTDtimeAvg.tolist(), myCTDparamAvg.tolist() #define fn to align one time stamp to another def alignCTD( myCTDtime1, myCTDparam1, myCTDtime2 ): myCTDparam1a = list() for i in range(0, len(myCTDtime2)): myDiff = list() for j in range(0, len(myCTDtime1)): myDiff.append( myCTDtime1[j] - myCTDtime2[i] ) if len(myDiff) > 0: ind = np.argmin(np.absolute(myDiff)) if np.absolute(myDiff[ind]) < 10*60*1000: myCTDparam1a.append(myCTDparam1[ind]) else: myCTDparam1a.append(np.nan) else: myCTDparam1a.append(np.nan) return myCTDparam1a def processTimeStamp(time): format = '%Y-%m-%dT%H:%M:%S' if 'AM' in time or 'PM' in time: format = '%m/%d/%Y %I:%M:%S %p' return int((datetime.strptime(time,format) - datetime(1970, 1, 1)).total_seconds() * 1000.0) def getBuoyData(buoy, myDirRT, start_time,end_time): #download buoy wave data into same folder and add to wavedata.txt if buoy =='PointSur': host=('https://www.ndbc.noaa.gov/data/realtime2/46239.spec') #gets last 45 days if buoy =='MontereyCanyon': host=('https://www.ndbc.noaa.gov/data/realtime2/46236.spec') #gets last 45 days #host=('https://www.ndbc.noaa.gov/data/realtime2/46279.spec') #gets last 45 days headers = { 'Accept-Encoding': 'gzip, deflate, sdch', 'Accept-Language': 'en-US,en;q=0.8', 'Upgrade-Insecure-Requests': '1', 'User-Agent': 'Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/56.0.2924.87 Safari/537.36', 'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8', 'Cache-Control': 'max-age=0', 'Connection': 'keep-alive',} wx=requests.get(host,headers).content status_code = requests.get(host,headers).status_code if status_code != 404: sl=wx.splitlines()#separate into lines wavefile=os.path.join(myDirRT, buoy + '_buoydata.txt') if os.path.isfile(wavefile): #get existing data with open(wavefile,'r') as file: #open in read mode lastline=file.readlines()[-1] #extract last line with open(wavefile,'a+') as file:#open in append mode for line in reversed(sl[2:]): #ignores two header lines datenow=datetime.strptime( line.decode('ascii').rstrip()[0:16],'%Y %m %d %H %M').replace(tzinfo=pytz.utc) lastdate=datetime.strptime( lastline[0:16],'%Y %m %d %H %M').replace(tzinfo=pytz.utc) if datenow>lastdate and datenow=start_time - timedelta(minutes=120): #write lines relevant to current dep file.write('\n'+line.decode('ascii').rstrip()) #load data, turning MM into nan data=pd.read_csv(open(wavefile,'r+'),sep='\s+', header=None, skiprows=2,keep_default_na=False,na_values=['MM','99.0','99.00','999.0','999'], names=['year','month','day','hour','minute','wvht','swh','swp','wwh','wwp', 'swd','wwd','steepness','apd','mwd'])#create headers to remove spaces for indexing data['date']=pd.to_datetime(data[['year','month','day','hour','minute']]) #get rid of white space data.columns=data.columns.str.strip() #convert time to hex first unix_time=(data.date - pd.Timestamp("1970-01-01")) // pd.Timedelta('1s')*1000 #add unix time to data frame data.insert(16,'unix_time',unix_time) #last column return data def getCTDDataRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") #CTD1 CTD1time = list() CTD1temp = list() CTD1sal = list() CTD1press = list() CTD1o2Hz = list() CTD1o2sol = list() #CTD2 CTD2time = list() CTD2temp = list() CTD2sal = list() CTD2press = list() CTD2o2Hz = list() CTD2o2sol = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Seabird CTD Records with D.O", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) withDo = convertToDict(zeepOutput, "GetReportDataResult") for record in withDo: if float(record['pressure (dbar)']) < 1.0 or float(record['pressure (dbar)']) > 3.0: #> 3 this accounts for bad pressure sensors CTD1time.append(processTimeStamp(record['timeStamp'])) CTD1temp.append(float(record['temperature (degC)'])) CTD1sal.append(float(record['salinity (PSU)'])) CTD1press.append(float(record['pressure (dbar)'])) CTD1o2Hz.append(float(record['oxygen (freq)'])) sol = gsw.O2sol_SP_pt(record['salinity (PSU)'], record['temperature (degC)'])/43.570; CTD1o2sol.append(sol) else: CTD2time.append(processTimeStamp(record['timeStamp'])) CTD2temp.append(float(record['temperature (degC)'])) CTD2sal.append(float(record['salinity (PSU)'])) CTD2press.append(float(record['pressure (dbar)'])) CTD2o2Hz.append(float(record['oxygen (freq)'])) sol = gsw.O2sol_SP_pt(record['salinity (PSU)'], record['temperature (degC)'])/43.570; CTD2o2sol.append(sol) start_date = end_date return CTD1time, CTD1temp, CTD1sal, CTD1press, CTD1o2Hz, CTD1o2sol, CTD2time, CTD2temp, CTD2sal, CTD2press, CTD2o2Hz, CTD2o2sol def getCTDwoO2DataRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") #CTD1 CTD1time = list() CTD1temp = list() CTD1sal = list() CTD1press = list() CTD1o2Hz = list() CTD1o2sol = list() #CTD2 CTD2time = list() CTD2temp = list() CTD2sal = list() CTD2press = list() CTD2o2Hz = list() CTD2o2sol = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Seabird CTD Records without D.O", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) withoutDo = convertToDict(zeepOutput, "GetReportDataResult") for record in withoutDo: if float(record['pressure (dbar)']) < 1.0 or float(record['pressure (dbar)']) > 3.0: #> 2 this accounts for bad pressure sensors CTD1time.append(processTimeStamp(record['timeStamp'])) CTD1temp.append(float(record['temperature (degC)'])) CTD1sal.append(float(record['salinity (PSU)'])) CTD1press.append(float(record['pressure (dbar)'])) CTD1o2Hz.append(np.nan) CTD1o2sol.append(np.nan) else: CTD2time.append(processTimeStamp(record['timeStamp'])) CTD2temp.append(float(record['temperature (degC)'])) CTD2sal.append(float(record['salinity (PSU)'])) CTD2press.append(float(record['pressure (dbar)'])) CTD1o2Hz.append(np.nan) CTD2o2sol.append(np.nan) start_date = end_date return CTD1time, CTD1temp, CTD1sal, CTD1press, CTD1o2Hz, CTD1o2sol, CTD2time, CTD2temp, CTD2sal, CTD2press, CTD2o2Hz, CTD2o2sol def getWindRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") METtime = list() METwindDir = list() METwindSpd = list() METairPress = list() METairTemp = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Wind Record 2", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: for record in data: METtime.append(processTimeStamp(record['timeStamp'])) METwindDir.append(float(record['avgWindDir(deg)'])) METwindSpd.append(float(record['avgWindSpeed(kt)']) * 0.514444) #knots to m/s METairPress.append(float(record['avgPress(mbar)'])) METairTemp.append(float(record['avgTemp(C)'])) start_date = end_date return METtime, METwindDir, METwindSpd, METairPress, METairTemp def getTelemRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") WGtime = list() WGlon = list() WGlat = list() WGdist = list() WGwaypoint = list() WGHead = list() WGDesHead = list() WGCOG = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Telemetry 6 Report", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: if type(data) is dict and "recordData" in data.keys(): recordData = data["recordData"] # Loop over each record for record in recordData: def bearing(lat1,lon1,lat2,lon2): dLon = lon2 - lon1; y = math.sin(dLon) * math.cos(lat2); x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(dLon); brng = np.rad2deg(math.atan2(y, x)); if brng < 0: brng+= 360 return brng WGtime.append(processTimeStamp(record['gliderTimeStamp'])) WGlon.append(float(record['longitude'])) WGlat.append(float(record['latitude'])) WGdist.append(int(record['distanceOverGround'])) WGwaypoint.append(int(record['targetWayPoint'])) WGHead.append(int(record['headingSubDegrees'])) WGDesHead.append(int(record['desiredBearingDegrees'])) try: WGCOG.append(bearing(WGlat[-2], WGlon[-2], WGlat[-1], WGlon[-1])) except IndexError: WGCOG.append(0) start_date = end_date return WGtime, WGlon, WGlat, WGdist, WGwaypoint, WGHead, WGDesHead, WGCOG def getPowermRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") PwrTime = list() PwrSolarPwrGen = list() PwrBatChargePwr = list() PwrOutPwrGen = list() PwrTotalBatPwr = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Amps Power Summary Report", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: if type(data) is dict and "recordData" in data.keys(): recordData = data["recordData"] # Loop over each record for record in recordData: PwrTime.append(processTimeStamp(record['gliderTimeStamp'])) PwrSolarPwrGen.append(int(record['solarPowerGenerated'])/1000) PwrBatChargePwr.append(int(record['batteryChargingPower'])/1000) PwrOutPwrGen.append(int(record['outputPortPower'])/1000) PwrTotalBatPwr.append(int(record['totalBatteryPower'])/1000) start_date = end_date return PwrTime, PwrSolarPwrGen, PwrBatChargePwr, PwrOutPwrGen, PwrTotalBatPwr def getEcoPuckRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") EPtime = list() EPcounts470 = list() EPcounts650 = list() EPcountsChl = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Eco Puck Data", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: # Loop over each record for record in data: EPtime.append(processTimeStamp(record['timeStamp'])) EPcounts470.append(float(record['column4'])) EPcounts650.append(float(record['column6'])) EPcountsChl.append(float(record['column8'])) start_date = end_date return EPtime, EPcounts470, EPcounts650, EPcountsChl def getVemcoStatRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") VR2Cstatus = '' start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Vemco Vr2c Status", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: # print(vr2cstatus) for record in data: VR2Cstatus = record['status'] start_date = end_date return VR2Cstatus def getVemcoTagRecord(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") VR2Ctime = list() VR2CtagID = list() VR2CtagReading = list() start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Vemco Vr2c Tags", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: for record in data: #print(record) VR2Ctime.append(int(record['time'])) a = record['tag Reading'].split(',') VR2CtagID.append(int(a[4])) #thought it was supposed to be a[5]? VR2CtagReading.append( record['kind'] +',' +record['latitude'] + ',' + record['longitude'] +',' +record['time'] +',"' +record['tagReading'] +'"' ) start_date = end_date return VR2Ctime, VR2CtagID, VR2CtagReading def getCO2Record(CCUnum,tBeg,tEnd): client = Client("https://wds.wgms.com/WGMSData.asmx?wsdl") CO2time = list() CO2equilOffCO2 = list() CO2equilOffTemp = list() CO2zeroOffCO2 = list() CO2airOffCO2 = list() CO2airOffTemp = list() pHtime = list() pHV1 = list() pHtimeInc = (np.arange(0,51,10)*60*1000).tolist() #activate this line for 10 min pH sampling - cwahl #pHtimeInc = (np.arange(-25,1,5)*60*1000).tolist() #JS 8/31/2016 otherwise get 'future' timestamps--DS 4/7/21 updated to 5min sampling, activate this line for 5 min pH sampling - cwahl CO2equil = np.zeros((0,6)) CO2zero = np.zeros((0,6)) CO2air = np.zeros((0,6)) CO2stnd = np.zeros((0,6)) CO2humEquilOff = list() #6 CO2humAirOff = list() #19 CO2humStndOff = list() #24 start_date = tBeg end_date = tBeg while start_date < tEnd: end_date = start_date + timedelta(hours=24) if end_date > tEnd: end_date = tEnd zeepOutput = client.service.GetReportData( token=wds_token, ReportName="Raw Data Report", SerialNo=CCUnum, StartDate=start_date.isoformat(), EndDate=end_date.isoformat(), ResultFormat=1) data = convertToDict(zeepOutput, "GetReportDataResult") if len(data) != 0: for record in data['recordData']: # print(str(record["structureType"]) + " " + str(record["payloadLength"])) if record["structureType"] == 0x91 and record["payloadLength"] == 212: temp = record["payload"].encode('utf-8') local_record = bytearray.fromhex(str(temp, 'utf-8')) local_record = ''.join(f"{n:02X}" for n in local_record) payload = local_record[156:420] #format = '%Y-%m-%dT%H:%M:%S.%f' #int_time = int((datetime.strptime(record["gliderTimeStamp"],format) - datetime(1970, 1, 1)).total_seconds() * 1000.0) int_time = struct.unpack('!i', binascii.unhexlify(payload[0:8]))[0] * 1000.0 equal_pump_on_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[8:16]))[0] equal_pump_on_temp = struct.unpack('!f', binascii.unhexlify(payload[16:24]))[0] equal_pump_on_pres = struct.unpack('!f', binascii.unhexlify(payload[24:32]))[0] equal_pump_off_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[32:40]))[0] equal_pump_off_temp = struct.unpack('!f', binascii.unhexlify(payload[40:48]))[0] equal_pump_off_pres = struct.unpack('!f', binascii.unhexlify(payload[48:56]))[0] equal_pump_off_humid = struct.unpack('!f', binascii.unhexlify(payload[56:64]))[0] zero_pump_on_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[64:72]))[0] zero_pump_on_temp = struct.unpack('!f', binascii.unhexlify(payload[72:80]))[0] zero_pump_on_pres = struct.unpack('!f', binascii.unhexlify(payload[80:88]))[0] zero_pump_off_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[88:96]))[0] zero_pump_off_temp = struct.unpack('!f', binascii.unhexlify(payload[96:104]))[0] zero_pump_off_pres = struct.unpack('!f', binascii.unhexlify(payload[104:112]))[0] air_pump_on_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[112:120]))[0] air_pump_on_temp = struct.unpack('!f', binascii.unhexlify(payload[120:128]))[0] air_pump_on_pres = struct.unpack('!f', binascii.unhexlify(payload[128:136]))[0] air_pump_off_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[136:144]))[0] air_pump_off_temp = struct.unpack('!f', binascii.unhexlify(payload[144:152]))[0] air_pump_off_pres = struct.unpack('!f', binascii.unhexlify(payload[152:160]))[0] air_pump_off_humid = struct.unpack('!f', binascii.unhexlify(payload[160:168]))[0] stand_flow_on_pres = struct.unpack('!f', binascii.unhexlify(payload[168:176]))[0] stand_flow_off_pCO2 = struct.unpack('!f', binascii.unhexlify(payload[176:184]))[0] stand_flow_off_temp = struct.unpack('!f', binascii.unhexlify(payload[184:192]))[0] stand_flow_off_pres = struct.unpack('!f', binascii.unhexlify(payload[192:200]))[0] stand_flow_off_humid = struct.unpack('!f', binascii.unhexlify(payload[200:208]))[0] ph1 = struct.unpack('!f', binascii.unhexlify(payload[208:216]))[0] ph2 = struct.unpack('!f', binascii.unhexlify(payload[216:224]))[0] ph3 = struct.unpack('!f', binascii.unhexlify(payload[224:232]))[0] ph4 = struct.unpack('!f', binascii.unhexlify(payload[232:240]))[0] ph5 = struct.unpack('!f', binascii.unhexlify(payload[240:248]))[0] ph6 = struct.unpack('!f', binascii.unhexlify(payload[248:256]))[0] can_humid = struct.unpack('!f', binascii.unhexlify(payload[256:264]))[0] CO2time.append(int(int_time)) CO2equilOffCO2.append(float(equal_pump_off_pCO2)) CO2equilOffTemp.append(float(equal_pump_off_temp)) CO2zeroOffCO2.append(float(zero_pump_off_pCO2)) CO2airOffCO2.append(float(air_pump_off_pCO2)) CO2airOffTemp.append(float(air_pump_off_temp)) #diagnostics chart CO2equil = np.vstack( [CO2equil, [ float(equal_pump_on_temp), float(equal_pump_off_temp), float(equal_pump_on_pres), float(equal_pump_off_pres), float(equal_pump_on_pCO2), float(equal_pump_off_pCO2)] ]) CO2zero = np.vstack( [CO2zero, [ float(zero_pump_on_temp), float(zero_pump_off_temp), float(zero_pump_on_pres), float(zero_pump_off_pres), float(zero_pump_on_pCO2), float(zero_pump_off_pCO2)] ]) CO2air = np.vstack( [CO2air, [ float(air_pump_on_temp), float(air_pump_off_temp), float(air_pump_on_pres), float(air_pump_off_pres), float(air_pump_on_pCO2), float(air_pump_off_pCO2)] ]) CO2stnd = np.vstack( [CO2stnd, [ np.nan, float(stand_flow_off_temp), float(stand_flow_on_pres), float(stand_flow_off_pres), np.nan, float(stand_flow_off_pCO2)] ]) #CO2 humidity chart CO2humEquilOff.append(float(equal_pump_off_humid)) CO2humAirOff.append(float(air_pump_off_humid)) CO2humStndOff.append(float(stand_flow_off_humid)) #create pH time stamp # for i in range(0,len(pHtimeInc)): # pHtime.append( CO2time[-1] + pHtimeInc[i] ) # pHV1.append( float(dataLine['data'][25 + i]) / 1000 ) for i in range(0,len(pHtimeInc)): pHtime.append( CO2time[-1] + pHtimeInc[i] ) pHV1.append( float(ph1) / 1000 ) pHV1.append( float(ph2) / 1000 ) pHV1.append( float(ph3) / 1000 ) pHV1.append( float(ph4) / 1000 ) pHV1.append( float(ph5) / 1000 ) pHV1.append( float(ph6) / 1000 ) start_date = end_date return CO2time, CO2equilOffCO2, CO2equilOffTemp, CO2zeroOffCO2, CO2airOffCO2, CO2airOffTemp, pHtime, pHV1, pHtimeInc, CO2equil, CO2zero, CO2air, CO2stnd, CO2humEquilOff, CO2humAirOff, CO2humStndOff #Generate JPG Report #----------------------------# def read_and_process_csv(file_path): dn = datetime(1970, 1, 1) # Read CSV skipping first row (header) df = pd.read_csv(file_path, skiprows=1, header=None) # Replace zeros with NaN df.replace(0, np.nan, inplace=True) # Convert first column from ms to days and add dn df[0] = df[0] / 86400 / 1000 # convert ms to days df[0] = df[0].apply(lambda x: dn + timedelta(days=x)) return df def interp_columns(position_times, position_vals, target_times): # Create interpolation function with bounds_error=False to allow extrapolation as NaN f = interp1d(position_times, position_vals, bounds_error=False, fill_value=np.nan) return f(target_times) def qcWG_plotAll(myDirPC): myCol = np.array([ [100, 100, 100], # gray [11, 99, 160], # blue [255, 127, 14], # orange [44, 160, 44], # green [50, 50, 50], # dark gray [134, 159, 192], # lt. blue [20, 20, 20] # black ]) / 255 pos = pd.DataFrame() co2 = pd.DataFrame() ctd = pd.DataFrame() ctd2 = pd.DataFrame() ep = pd.DataFrame() met = pd.DataFrame() metw = pd.DataFrame() co2equil = pd.DataFrame() co2zero = pd.DataFrame() co2air = pd.DataFrame() co2stnd = pd.DataFrame() co2hum = pd.DataFrame() timeAll = np.empty(0) timeAll = timeAll.astype('datetime64[ns]') if Path(myDirPC + 'Position.csv').is_file(): pos = read_and_process_csv(myDirPC + 'Position.csv') pos_times = pos[0].values.astype('datetime64[ns]').astype(np.float64) # convert to float for interp pos_lat = pos[1].values pos_lon = pos[2].values timeAll = np.concatenate([timeAll, pos.iloc[:,0]]) if Path(myDirPC + 'CO2.csv').is_file(): co2 = read_and_process_csv(myDirPC + 'CO2.csv') co2_times = co2[0].values.astype('datetime64[ns]').astype(np.float64) co2[5] = interp_columns(pos_times, pos_lat, co2_times) co2[6] = interp_columns(pos_times, pos_lon, co2_times) timeAll = np.concatenate([timeAll, co2.iloc[:,0]]) if Path(myDirPC + 'CTD.csv').is_file(): ctd = read_and_process_csv(myDirPC + 'CTD.csv') ctd_times = ctd[0].values.astype('datetime64[ns]').astype(np.float64) ctd[5] = interp_columns(pos_times, pos_lat, ctd_times) ctd[6] = interp_columns(pos_times, pos_lon, ctd_times) timeAll = np.concatenate([timeAll, ctd.iloc[:,0]]) if Path(myDirPC + 'CTD2.csv').is_file(): ctd2 = read_and_process_csv(myDirPC + 'CTD2.csv') ctd2_times = ctd2[0].values.astype('datetime64[ns]').astype(np.float64) ctd2[5] = interp_columns(pos_times, pos_lat, ctd2_times) ctd2[6] = interp_columns(pos_times, pos_lon, ctd2_times) timeAll = np.concatenate([timeAll, ctd2.iloc[:,0]]) if Path(myDirPC + 'EcoPuck.csv').is_file(): ep = read_and_process_csv(myDirPC + 'EcoPuck.csv') ep_times = ep[0].values.astype('datetime64[ns]').astype(np.float64) ep[5] = interp_columns(pos_times, pos_lat, ep_times) ep[6] = interp_columns(pos_times, pos_lon, ep_times) timeAll = np.concatenate([timeAll, ep.iloc[:,0]]) if Path(myDirPC + 'MET.csv').is_file(): met = read_and_process_csv(myDirPC + 'MET.csv') met_times = met[0].values.astype('datetime64[ns]').astype(np.float64) met[4] = interp_columns(pos_times, pos_lat, met_times) met[5] = interp_columns(pos_times, pos_lon, met_times) timeAll = np.concatenate([timeAll, met.iloc[:,0]]) if Path(myDirPC + 'METwind.csv').is_file(): metw = read_and_process_csv(myDirPC + 'METwind.csv') metw_times = metw[0].values.astype('datetime64[ns]').astype(np.float64) metw[4] = interp_columns(pos_times, pos_lat, metw_times) metw[5] = interp_columns(pos_times, pos_lon, metw_times) timeAll = np.concatenate([timeAll, metw.iloc[:,0]]) # TagID commented out in Matlab code due to errors if no data # TagID = pd.read_csv(myDirPC + 'tagID.csv', skiprows=1, header=None) if Path(myDirPC + 'CO2equil.csv').is_file(): co2equil = read_and_process_csv(myDirPC + 'CO2equil.csv') if Path(myDirPC + 'CO2zero.csv').is_file(): co2zero = read_and_process_csv(myDirPC + 'CO2zero.csv') if Path(myDirPC + 'CO2air.csv').is_file(): co2air = read_and_process_csv(myDirPC + 'CO2air.csv') if Path(myDirPC + 'CO2stnd.csv').is_file(): co2stnd = read_and_process_csv(myDirPC + 'CO2stnd.csv') if Path(myDirPC + 'CO2hum.csv').is_file(): co2hum = read_and_process_csv(myDirPC + 'CO2hum.csv') #no data found so leave. if timeAll.size == 0: return xmin, xmax = np.min(timeAll), np.max(timeAll) fig = plt.figure(figsize=(22, 15)) gs = GridSpec(7, 1, figure=fig, hspace=0.25, top=0.85, bottom=0.15) # Axes 1 ax1 = fig.add_subplot(gs[0]) ax1.set_ylabel('Lon [\u00B0E]') ax1.set_xlim(xmin, xmax) ax1.tick_params(labelbottom=False) ax1.spines['top'].set_visible(False) ax1.spines['right'].set_visible(False) ax1.yaxis.label.set_color(myCol[0]) ax1a = ax1.twinx() ax1a.set_ylabel('Lat [\u00B0N]') ax1a.tick_params(labelbottom=False) ax1a.spines['top'].set_visible(False) ax1a.spines['left'].set_visible(False) ax1a.yaxis.set_label_position("right") ax1a.yaxis.tick_right() ax1a.yaxis.label.set_color(myCol[1]) ax2 = fig.add_subplot(gs[1]) ax2.set_ylabel('T [\u00B0C]') ax2.set_xlim(xmin, xmax) ax2.tick_params(labelbottom=False) ax2.spines['top'].set_visible(False) ax2.spines['right'].set_visible(False) ax2.yaxis.label.set_color(myCol[1]) ax2a = ax2.twinx() ax2a.set_ylabel('Sal') ax2a.tick_params(labelbottom=False) ax2a.spines['top'].set_visible(False) ax2a.spines['left'].set_visible(False) ax2a.yaxis.set_label_position("right") ax2a.yaxis.tick_right() ax2a.yaxis.label.set_color(myCol[2]) ax2b = ax2.twinx() ax2b.set_ylabel('O2 [\u03BCmol/kg]') ax2b.spines['top'].set_visible(False) ax2b.spines['left'].set_visible(False) ax2b.yaxis.set_label_position("right") ax2b.yaxis.tick_right() ax2b.yaxis.label.set_color(myCol[3]) ax2b.spines['right'].set_position(('outward', 40)) # ax2b1 = ax2.twinx() # ax2b1.set_ylabel('O2 [\u03BCmol/kg]') # ax2b1.spines['top'].set_visible(False) # ax2b1.spines['left'].set_visible(False) # ax2b1.yaxis.set_label_position("right") # ax2b1.yaxis.tick_right() # ax2b1.yaxis.label.set_color(myCol[3]) # # Offset ax2b1 to the right # ax2b1.spines['right'].set_position(('outward', 40)) ax3 = fig.add_subplot(gs[2]) ax3.set_ylabel('T [\u00B0C]') ax3.set_xlim(xmin, xmax) ax3.tick_params(labelbottom=False) ax3.spines['top'].set_visible(False) ax3.spines['right'].set_visible(False) ax3.yaxis.label.set_color(myCol[1]) ax3a = ax3.twinx() ax3a.set_ylabel('Sal') ax3a.tick_params(labelbottom=False) ax3a.spines['top'].set_visible(False) ax3a.spines['left'].set_visible(False) ax3a.yaxis.set_label_position("right") ax3a.yaxis.tick_right() ax3a.yaxis.label.set_color(myCol[2]) ax3b = ax3.twinx() ax3b.set_ylabel('O2 [\u03BCmol/kg]') ax3b.spines['top'].set_visible(False) ax3b.spines['left'].set_visible(False) ax3b.yaxis.set_label_position("right") ax3b.yaxis.tick_right() ax3b.yaxis.label.set_color(myCol[3]) ax3b.spines['right'].set_position(('outward', 50)) ax4 = fig.add_subplot(gs[3]) ax4.set_ylabel('pH') ax4.set_xlim(xmin, xmax) ax4.tick_params(labelbottom=False) ax4.spines['top'].set_visible(False) ax4.spines['right'].set_visible(False) ax4.yaxis.label.set_color(myCol[1]) ax4a = ax4.twinx() ax4a.set_ylabel('pCO2 Water [ppm]') ax4a.yaxis.label.set_color(myCol[1]) ax4a.tick_params(labelbottom=False) ax4a.spines['top'].set_visible(False) ax4a.spines['left'].set_visible(False) ax4a.yaxis.set_label_position("right") ax4a.yaxis.tick_right() ax4a.yaxis.label.set_color(myCol[2]) ax4b = ax4.twinx() ax4b.set_ylabel('pCO2 Air [ppm]') ax4b.yaxis.label.set_color(myCol[1]) ax4b.tick_params(labelbottom=False) ax4b.spines['top'].set_visible(False) ax4b.spines['left'].set_visible(False) ax4b.yaxis.set_label_position("right") ax4b.yaxis.tick_right() ax4b.yaxis.label.set_color(myCol[3]) ax4b.spines['right'].set_position(('outward', 50)) ax5 = fig.add_subplot(gs[4]) ax5.set_ylabel('Fluor [\u03BCg/l]') ax5.set_xlim(xmin, xmax) ax5.yaxis.label.set_color(myCol[3]) ax5.tick_params(labelbottom=False) ax5.spines['top'].set_visible(False) ax5.spines['right'].set_visible(False) ax5a = ax5.twinx() ax5a.set_ylabel('bb470 [1/m]') ax5a.tick_params(labelbottom=False) ax5a.spines['top'].set_visible(False) ax5a.spines['left'].set_visible(False) ax5a.yaxis.set_label_position("right") ax5a.yaxis.tick_right() ax5a.yaxis.label.set_color(myCol[5]) ax5b = ax5.twinx() ax5b.set_ylabel('bb650 [1/m]') ax5b.tick_params(labelbottom=False) ax5b.spines['top'].set_visible(False) ax5b.spines['left'].set_visible(False) ax5b.yaxis.set_label_position("right") ax5b.yaxis.tick_right() ax5b.yaxis.label.set_color(myCol[6]) ax5b.spines['right'].set_position(('outward', 50)) ax6 = fig.add_subplot(gs[5]) ax6.set_ylabel('Wind Dir') ax6.set_xlim(xmin, xmax) ax6.tick_params(labelbottom=False) ax6.spines['top'].set_visible(False) ax6.spines['right'].set_visible(False) ax6a = ax6.twinx() ax6a.set_ylabel('Speed [m/s]') ax6a.yaxis.label.set_color(myCol[1]) ax6a.tick_params(labelbottom=False) ax6a.spines['top'].set_visible(False) ax6a.spines['left'].set_visible(False) ax6a.yaxis.set_label_position("right") ax6a.yaxis.tick_right() ax7 = fig.add_subplot(gs[6]) ax7.set_ylabel('Press [hPa]') ax7.set_xlim(xmin, xmax) if len(pos.index) > 1: ax1.plot(pos.iloc[:,0], pos.iloc[:,1], '.-', color=myCol[0]) ax1.set_ylim(np.min(pos.iloc[:,1]), np.max(pos.iloc[:,1])) ax1a.plot(pos.iloc[:,0], pos.iloc[:,2], '.-', color=myCol[1]) ax1a.set_ylim(np.min(pos.iloc[:,2]), np.max(pos.iloc[:,2])) # Axes 2 if len(ctd.index) > 1: if ctd.iloc[:,1].count() > 1: ax2.plot(ctd.iloc[:,0], ctd.iloc[:,1], '.-', color=myCol[1]) ax2.set_ylim(np.min(ctd.iloc[:,1]), np.max(ctd.iloc[:,1])) if ctd.iloc[:,2].count() > 1: ax2a.plot(ctd.iloc[:,0], ctd.iloc[:,2], '.-', color=myCol[2]) ax2a.set_ylim(np.min(ctd.iloc[:,2]), np.max(ctd.iloc[:,2])) if ctd.iloc[:,3].count() > 1: ax2b.plot(ctd.iloc[:,0], ctd.iloc[:,3], '.-', color=myCol[3]) ax2b.set_ylim(np.min(ctd.iloc[:,3]), np.max(ctd.iloc[:,3])) # Axes 3 if len(ctd2.index) > 1: if ctd2.iloc[:,1].count() > 1: ax3.plot(ctd2.iloc[:,0], ctd2.iloc[:,1], '.-', color=myCol[1]) ax3.set_ylim(np.min(ctd2.iloc[:,1]), np.max(ctd2.iloc[:,1])) if ctd2.iloc[:,2].count() > 1: ax3a.plot(ctd2.iloc[:,0], ctd2.iloc[:,2], '.-', color=myCol[2]) ax3a.set_ylim(np.min(ctd2.iloc[:,2]), np.max(ctd2.iloc[:,2])) if ctd2.iloc[:,3].count() > 1: ax3b.plot(ctd2.iloc[:,0], ctd2.iloc[:,3], '.-', color=myCol[3]) ax3b.set_ylim(np.min(ctd2.iloc[:,3]), np.max(ctd2.iloc[:,3])) # Axes 4 if len(co2.index) > 1: if co2.iloc[:,1].count() > 1: ax4.plot(co2.iloc[:,0], co2.iloc[:,1], '.-', color=myCol[1]) ax4.set_ylim(np.nanmin(co2.iloc[:,1]), np.nanmax(co2.iloc[:,1])) #f = (co2.iloc[:,2]).copy() #f[np.isnan(f)] = 0 f_mask = ~np.isnan(co2.iloc[:,2]) co2_all_y = np.empty(0) if f_mask.count() > 1: ax4a.plot(co2[f_mask][0], co2[f_mask][2], '.-', color=myCol[2]) co2_all_y = np.concatenate([co2_all_y, co2[f_mask][2]]) f_mask = ~np.isnan(co2.iloc[:,3]) if f_mask.count() > 1: ax4b.plot(co2[f_mask][0], co2[f_mask][3], '.-', color=myCol[3]) co2_all_y = np.concatenate([co2_all_y, co2[f_mask][3]]) co2_ymin = np.nanmin(co2_all_y) co2_ymax = np.nanmax(co2_all_y) ax4a.set_ylim(co2_ymin, co2_ymax) ax4b.set_ylim(co2_ymin, co2_ymax) # Axes 5 if ep.iloc[:,1].count() > 1 and ep.iloc[:,2].count() > 1 and ep.iloc[:,3].count() > 1 and ep.iloc[:,4].count() > 1: ax5.plot(ep.iloc[:,0], ep.iloc[:,1], '.-', color=myCol[3]) ax5.set_ylim(np.min(ep.iloc[:,1]), np.max(ep.iloc[:,1])) ax5a.plot(ep.iloc[:,0], ep.iloc[:,2], '.-', color=myCol[5]) ax5a.set_ylim(np.nanmin(ep.iloc[:,2:4]), np.nanmax(ep.iloc[:,2:4])) ax5b.plot(ep.iloc[:,0], ep.iloc[:,3], '.-', color=myCol[6]) ax5b.set_ylim(np.nanmin(ep.iloc[:,2:4]), np.nanmax(ep.iloc[:,2:4])) # Axes 6 if len(metw.index) > 1: ax6.plot(metw.iloc[:,0], metw.iloc[:,1], '.', color=myCol[0]) ax6.yaxis.label.set_color(myCol[0]) ax6.set_ylim(np.min(metw.iloc[:,1]), np.max(metw.iloc[:,1])) ax6a.plot(metw.iloc[:,0], metw.iloc[:,2], '.-', color=myCol[1]) ax6a.set_ylim(np.min(metw.iloc[:,2]), np.max(metw.iloc[:,2])) # Axes 7 if len(met.index) > 1: ax7.plot(met.iloc[:,0], met.iloc[:,1], '.-', color=myCol[0]) ax7.yaxis.label.set_color(myCol[0]) ax7.set_ylim(np.min(met.iloc[:,1]), np.max(met.iloc[:,1])) #plt.show() plt.savefig(myDirPC +'ReportPlots.jpg', dpi=300, bbox_inches='tight')