codar/getcodar_drifter.py at master · xhx509/codar · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
# -*- coding: utf-8 -*-
"""
Created on Fri Feb 01 14:15:21 2013

@author: huanxin
"""
########################################################################
#plot codar and drifter,
#It gets input value from control file "getcodar_bydrifter_ctl.txt" or
# a python function "getcodar_ctl_file_py()"
#From lat lon of drifter data, it gets a range area first,
# then it gets codar and sst data in this range.
# At last, it plots a picture or pictures
#Input values:datetime_wanted,filename,driftnumber,url,model_option,num,interval_dtime,interval,step_size
#output values:gbox,drifter_data,id,lat_vel,lon_vel,u,v
#function uses:getcodar_ctl_file,getdrift_raw_range_latlon,getcodar_ctl_lalo,getcodar_ctl_id,getdrift_raw,getcodar_edge
########################################################################

"""
Notice: make sure to select appropriate data of codar and drifter, especially same time period.
"""
import pytz
from matplotlib.dates import date2num, num2date
import datetime
import pylab
import sys
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import numpy as np

pydir='../'

sys.path.append(pydir)
from hx import getcodar_ctl_file,getdrift_raw_range_latlon,getcodar_ctl_lalo,getcodar_ctl_id,getdrift_raw,getcodar_edge
###############################################
#### input #########################################################

utc = pytz.timezone('UTC')
png_num=0 #for save picture
inputfilename='./getcodar_bydrifter_ctl.txt'
#id=str(int(id))
(datetime_wanted,filename,driftnumber,url,model_option,num,interval_dtime,interval,step_size)=getcodar_ctl_file(inputfilename)

id3=int(driftnumber) #change format for id
datetime_wanted_1=datetime_wanted
(maxlon,minlon,maxlat,minlat)=getdrift_raw_range_latlon(filename,id3,interval,datetime_wanted_1,num,step_size)
for i in range(5): #make sure the picture can show lat and lon clearly
    if maxlat-minlat<=0.1:
        maxlat=maxlat+0.02
        minlat=minlat-0.02
    if maxlon-minlon<=0.1:
        maxlon=maxlon+0.02
        minlon=minlon-0.02

(lat_max_i,lon_max_i,lat_min_i,lon_min_i)=getcodar_ctl_lalo(model_option,maxlat,maxlon,minlat,minlon)# get index of max min lat lon
gbox=[minlon-0.03,maxlon+0.03, minlat-0.03, maxlat+0.03] # get edge for get sst
for x in range(num):

  id=getcodar_ctl_id(model_option,url,datetime_wanted) #get index of codar
  ask_input=num2date(datetime_wanted) #get time for getsst
  (drifter_data)=getdrift_raw(filename,id3,interval,datetime_wanted) #get drifter data
  lon=drifter_data['lon']
  lat=drifter_data['lat']
  (lat_vel,lon_vel,u,v)=getcodar_edge(url,id,lat_max_i,lon_max_i,lat_min_i,lon_min_i) #get codar data
  id=str(id)
  idg1=list(ml.find(np.array(u)<>-999.0/100.))
  idg2=list(ml.find(np.array(lat_vel)>=minlat))
  idg12=list(set(idg1).intersection(set(idg2)))
  idg3=list(ml.find(np.array(lon_vel)>=minlon))
  idg=list(set(idg12).intersection(set(idg3))) #get index for codar data based on edge and id(time)

  png_num=png_num+1 #for save movie picture

  if len(idg)<>0:

    fig = plt.figure()

    ax = fig.add_subplot(111)

    plt.title(str(num2date(datetime_wanted).strftime("%d-%b-%Y %H"))+'h') #plot title
    pylab.ylim([minlat-0.03,maxlat+0.03]) #limit edge
    pylab.xlim([minlon-0.03,maxlon+0.03])
    #plot codar
    codar_scale_index=np.average(np.average(np.reshape(u,np.size(u))[idg]),np.average(np.reshape(v,np.size(v))[idg]))
    q=plt.quiver(np.reshape(lon_vel,np.size(lon_vel))[idg],np.reshape(lat_vel,np.size(lat_vel))[idg],np.reshape(u,np.size(u))[idg],np.reshape(v,np.size(v))[idg],angles='xy',scale=2,color='b',label='codar')
    p = plt.quiverkey(q,minlon-0.01,maxlat+0.06,0.2,str(round(0.2,2))+"m/s",coordinates='data',color='b') #plot key
    #codar_scale_index/10 get the suit size of arrow
    lat_wanted=lat[-1]
    lon_wanted=lon[-1]

    plt.plot(lon_wanted,lat_wanted,'.',markersize=20,color='r',label='end')#end time
    plt.plot(np.reshape(lon,np.size(lon)),np.reshape(lat,np.size(lat)),color='black')
    plt.plot(lon[0],lat[0],'.',markersize=20,color='g',label='start') # start time
    plt.legend( numpoints=1)

    plt.title(str(num2date(datetime_wanted).strftime("%d-%b-%Y %H"))+'h')

  else:
    print "Sorry. No good codar data at this time"

    plt.title(str(num2date(datetime_wanted).strftime("%d-%b-%Y %H"))+'h')
    lat_wanted=lat[-1]
    lon_wanted=lon[-1]
    #find wanted point lat,lon

    plt.plot(lon_wanted,lat_wanted,'.',markersize=30,color='r',label='end')
    plt.plot(np.reshape(lon,np.size(lon)),np.reshape(lat,np.size(lat)),color='black')
    pylab.ylim([minlat-0.02,maxlat+0.02])
    pylab.xlim([minlon-0.02,maxlon+0.02])
    plt.plot(lon[0],lat[0],'.',markersize=20,color='g',label='start') # start time
    plt.legend()
  #for another forloop
  id=int(id)+interval
  plt.savefig('./'+str('%03d' % png_num) + '.png')
  datetime_wanted=date2num(num2date(datetime_wanted)+datetime.timedelta( 0,step_size*60*60 ))
  plt.show()