[ncl-talk] Winds overlaid with shipping track data

Mary Haley haley at ucar.edu
Wed Aug 5 09:18:05 MDT 2015


Hi Melissa,

Sorry about the delay. I've got a lot of emails to deal with.

I think that you can get rid of many of your do loops, and that really, you
should only need to loop across time, and then add all the markers at one
point to each time step.

Right now, it looks like you are creating one plot at a time, and only
adding one marker to each plot. Is this what you want?

I'm thinking your code should look something more like this:

data_indexes = ind(year  .ge. 1881 .and. year  .le. 1885 .and.\
                   month .ge.    8 .and. month .le.   19 .and.\
                   day   .ge.    0 .and. day   .le.   30 .and.\
                   hour  .ge.    0 .and. hour  .le.   29)

lons = slon(data_indexes)
lats = slat(data_indexes)

do gg = 0,dimsizes(time)-1
   plot(gg) = gsn_csm_map_ce(wks,mapres)
   plot1(gg)= gsn_csm_vector(wks,u(gg,:,:),v(gg,:,:),res)
   overlay(plot(gg),plot1(gg))
   plot2(gg) = gsn_add_polymarker(wks,plot(gg),lons,lats,polyres)
end do

Note that I don't think the "data_indexes" line is quite correct, because I
don't understand why your months would go from 8 to 19.  I would have to
see your various year, month, etc arrays to better understand what time
range you are trying to select.

But, hopefully you understand what the above code is doing, and you can fix
it as necessary.

--Mary


On Mon, Aug 3, 2015 at 9:34 AM, Melissa Lazenby <M.Lazenby at sussex.ac.uk>
wrote:

> Hi Mary
>
> Sorry I was not specific enough.
> The main problem is with the loops to create the conditions when reading
> in the text file.
> Below is a snippet of the code that isn't running how I would like it too.
>
> yrs=4
> mons=12
> days=31
> hrs=23
>
> do j=0, yrs-1
>    jj=j+1881
>   do i=0, mons-1
>      ii=i+8
>     do m=0, days-1
>       do n = 0, hrs+6
>
> data_index = ind((year .eq. jj) .and. (month .eq. ii) .and. (day .eq. m)
> .and. (hour .eq. n))
>
>   if(.not. any(ismissing(data_index))) then
>   x=slon(data_index)
>   y=slat(data_index)
>
>     end if
>
> do gg = 0,dimsizes(time)-1
>
>  plot(gg) = gsn_csm_map_ce(wks,mapres)
>  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,y,polyres)
>
> end do
>  end do
>    end do
>       end do
>         end do
>
> When I put actual values in for jj, ii, m and n, the code runs perfectly
> and plots just those shipping tracks I want according to the conditional
> statement. But because there is a lot of data, I want to create loops that
> go through the years, months, days and hours (only every 6 hours) and the
> output is not looping through the different winds values and different
> shipping track data. It either gets the wind to change and the ship tracks
> to stay the same or vise versa.
>
> I do not have much experience with how to write do loops and I think my
> syntax is not correct so if there are any examples to help me or your
> advice would be great.
>
> Many thanks!
>
> Kindest Regards
> Melissa
>
> ------------------------------
> *From:* Mary Haley [haley at ucar.edu]
> *Sent:* 03 August 2015 16:09
> *To:* Melissa Lazenby
> *Cc:* ncl-talk at ucar.edu
> *Subject:* Re: [ncl-talk] Winds overlaid with shipping track data
>
> Melissa,
>
> Could you be more specific than "it's not working correctly?"  It's hard
> for us to debug a problem when we don't know where to look.
>
> It's also helpful to include a sample image with an explanation of what's
> wrong with it, or an error message.
>
> --Mary
>
>
> On Fri, Jul 31, 2015 at 10:28 AM, Melissa Lazenby <M.Lazenby at sussex.ac.uk>
> wrote:
>
>> 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
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150805/5cd2216d/attachment.html 


More information about the ncl-talk mailing list