[ncl-talk] Winds overlaid with shipping track data
Melissa Lazenby
M.Lazenby at sussex.ac.uk
Fri Jul 31 10:28:12 MDT 2015
Hi NCL-Users
I am struggling to get my loops to work in a plot I am trying to create for an animation.
I have u and v wind data that has 7304 timesteps in it and I have a text file of shipping data that I want to use to plot ship tracks for the same wind fields i.e I would like to create 7304 png images to use to create an animation using fmpeg.
The text file of the shipping data has 252391 timesteps so I only want to extract the data that matches the wind data. So I want to plot ship tracks only every 6 hours.
The wind data is 4 times daily at 00,06,12,18
The ship data is hourly 00,01,02,03,04 etc..
I want to create plots to match up with the wind fields i.e. use the winds as the background plot and overlay them with ship tracks using polymarker function in NCL for every 6 hours.
Below is my script but it is not working correctly.
If anyone could please advice how to go about fixing it that would be much appreciated.
Many thanks!
Kindest Regards
Melissa
;----------------------------------------------------------------------
; Krakatoa_1.ncl
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
f=addfile("/mnt/geog/ml382/McGill_Project/ReanalysisV2/uwnd_krakatoa.nc", "r")
g=addfile("/mnt/geog/ml382/McGill_Project/ReanalysisV2/vwnd_krakatoa.nc", "r")
latS = -70.
latN = 50.
lonL = 0.
lonR = 150.
u = f->uwnd(:,0,:,:)
v = g->vwnd(:,0,:,:)
lat = f->lat
lon = f->lon
time = f->time
speed= new(dimsizes(u),float)
do gg = 0,dimsizes(time)-1
speed(gg,:,:) = sqrt(u(gg,:,:)^2+v(gg,:,:)^2)
speed!0 = "time"
speed!1 = "lat"
speed!2 = "lon"
speed&time= time
speed&lat = lat
speed&lon = lon
speed at long_name = "Wind Speed"
speed at units = "m/s"
wks = gsn_open_wks("X11","Krakatoa_Animation_"+time(gg)+"") ; specifies a plot
plot = new (dimsizes(time),"graphic")
plot1= new (dimsizes(time),"graphic")
plot2= new (dimsizes(time),"graphic")
x= new (dimsizes(time),"float")
y= new (dimsizes(time),"float")
;---Create plot
res = True ; plot mods desired
res at cnFillOn = True ; turn on color for contours
res at cnLinesOn = False ; turn off contour lines
res at cnLineLabelsOn = False ; turn off contour line labels
res at gsnScalarContour = True ; contours desired
res at gsnDraw = False ; do not draw the plot
res at gsnFrame = False
res at vcRefMagnitudeF = 2.0 ; define vector ref mag
res at vcRefLengthF = 0.008 ; define length of vec ref
res at vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
res at vcRefAnnoArrowLineColor = "black" ; change ref vector color
res at vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
res at vcGlyphStyle = "CurlyVector" ; turn on curley vectors
res at vcLineArrowColor = "black" ; change vector color
res at vcLineArrowThicknessF = 2.0 ; change vector thickness
res at vcVectorDrawOrder = "PostDraw" ; draw vectors last
res at vcMinDistanceF = 0.02
res at gsnLeftString = ""
res at gsnMainString = ""
res at gsnRightString = ""
res at mpMinLonF = lonL ; select a subregion
res at mpMaxLonF = lonR
res at mpMinLatF = latS
res at mpMaxLatF = latN
pres = True
pres at mpFillOn = False
pres at cnFillOn = True ; color on
pres at cnLinesOn = False ; turn off contour lines
pres at gsnScalarContour = True ; vectors over contours
pres at gsnSpreadColors = True ; use full colormap
pres at gsnAddCyclic = False
pres at gsnDraw = False ; do not draw the plot
pres at gsnFrame = False
pres at lbLabelBarOn = False
pres at gsnMaximize = True
pres at mpMinLonF = lonL ; select a subregion
pres at mpMaxLonF = lonR
pres at mpMinLatF = latS
pres at mpMaxLatF = latN
pres at gsnLeftString = ""
pres at gsnMainString = ""
pres at gsnRightString = ""
pres at cnLevelSelectionMode= "ManualLevels" ; manual levels
pres at cnMinLevelValF = -100 ; min level
pres at cnMaxLevelValF = 120 ; max level
pres at cnLevelSpacingF = 10 ; contour spacing
mres = True
mres at gsMarkerIndex = 3 ; polymarker style
mres at gsMarkerSizeF = 20. ; polymarker size
mres at gsMarkerColor = "red"
mpres = True ; plot mods desired
mpres at xyMarkLineModes = "Markers" ; choose which have markers
mpres at xyMarkers = 16 ; choose type of marker
mpres at xyMarkerColor = "red" ; Marker color
mpres at xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
;Ship track data
ship = "/mnt/geog/ml382/McGill_Project/ReanalysisV2/Krakatoa_ID.txt"
data = asciiread(ship,-1,"string")
year = stringtointeger(str_get_cols(data,1,5))
month = stringtointeger(str_get_cols(data,8,9))
day = stringtointeger(str_get_cols(data,11,13))
hour = stringtointeger(str_get_cols(data,15,20))
slat = stringtointeger(str_get_cols(data,22,28))
slon = stringtointeger(str_get_cols(data,31,37))
k=0
do j=0, 4
jj=j+1881
do i=0, 11
do m=0, 30
n = 0
do nn=0, 7304
data_index = ind((year(k) .eq. jj) .and. (month(k) .eq. i) .and. (day(k) .eq. m) .and. (hour(k) .eq. n))
if(.not. any(ismissing(data_index))) then
x(k)=slon(data_index)
y(k)=slat(data_index)
print(x)
print(y)
end if
polyres = True ; poly marker mods desired
polyres at gsMarkerIndex = 16 ; choose circle as polymarker
polyres at gsMarkerSizeF = 7.0 ; select size to avoid streaking
polyres at gsMarkerColor = "red" ; choose color
pres at gsnLeftString = "Krakatoa: Winds and Shipping Tracks"
plot(gg) = gsn_csm_contour_map(wks,speed(gg,:,:),pres)
plot1(gg)= gsn_csm_vector(wks,u(gg,:,:),v(gg,:,:),res)
overlay(plot(gg),plot1(gg))
plot2(gg) = gsn_add_polymarker(wks,plot(gg),x(k),y(k),polyres) ; draw polymarkers
draw(plot(gg))
frame(wks)
k=k+1
n=nn+6
end do
end do
end do
end do
end do
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150731/4c5d7ae2/attachment.html
More information about the ncl-talk
mailing list