[ncl-talk] Time series plot
Toni Klemm
toni-klemm at tamu.edu
Tue Jan 29 11:54:42 MST 2019
Hi!
Welcome to NCL! This task might be quite a steep learning curve, but it is possible. I have not worked with your exact kind of data, but others might chime in with more specific steps and examples.
To separate your data into severe and very severe, you could load the original data twice and use an if-condition and do-loops (page 2 here <https://www.ncl.ucar.edu/Document/Reference_Cards/NCL_scripting_language_reference_card_legal.pdf>) based on the wind speed number to check every data point and delete data points that do not meet the criteria. Then save those two variables with your original time/latitude/longitude coordinates in a new file. To plot them both in one map <https://www.ncl.ucar.edu/Applications/vector.shtml> you can use an overlay plot <https://www.ncl.ucar.edu/Applications/overlay.shtml>, which means you plot one map first and then another map on top of it with the lines/arrows in a different color. Finally, add a legend to define the two colors.
Below and attached are a very basic NCL script which I used to create a spatial subset of global daily precipitation data. It might give you a general idea of the structure of an NCL script. First you load the data, then you manipulate it in the way you need to, then you save it in a new file, and in this case it finally plots the data (however, not using an overlay plot).
I hope this will help to get you started.
Best of luck,
Toni
Toni Klemm, Ph.D.
Postdoctoral Research Associate
Department of Ecosystem Science and Management
College of Agriculture and Life Sciences
Texas A&M University, College Station, TX
Contributor to the Early Career Climate Forum <http://www.eccforum.org/>
www.toni-klemm.de <http://www.toni-klemm.de/> | @toniklemm <http://twitter.com/toniklemm>
; ---------------------------------------------------
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"
; ---------------------------------------------------
; ***** GFDL FLOR-B01 daily precipitation data ******
; ---------------------------------------------------
begin
do year = 1980,2013 ; do-loop through the record years, 1980-2013
do m = 1,12 ; do-loop through the 12 calendar months
do r = 1,12 ; do-loop through the 12 model runs
if (m.lt.10) then
m_zero = "0" + m ; turn "1"..."9" into "01"..."09" to match filename structure (months have 2 digits)
end if
if(m.gt.9) then
m_zero = m ; leave "10..."12" unchanged
end if
print("**** START: " + year + m_zero + ", run " + r + " ****")
; ***********************************************
; READ IN MULTIPLE FILES
; ***********************************************
nmme_files = systemfunc("ls pr_day_GFDL*-v3.1-" + m_zero + year + "_r" + r + "i1p1" + "*.nc")
nmme_f = addfiles(nmme_files,"r")
ListSetType(nmme_f,"join")
print("Finished: 1/8 Read in multiple files")
; ***********************************************
; READ IN A SINGLE FILE
; ***********************************************
; year = 1980
; m = 3
; r = 1
; nmme_f = addfile("pr_day_GFDL-FLORB01_FLORB01-P1-ECDA-v3.1-031980_r1i1p1_19800301-19810228.nc","r")
lat = nmme_f[0]->lat
lon = nmme_f[0]->lon
; print("Finished: 1/8 Read in a single file")
; ***********************************************
; PRINT THE DIMENSIONS OF THE READ-IN FILE
; ***********************************************
; print("the dimensions of latitude are:")
; print(dimsizes(lat))
; print("the dimensions of longitude are:")
; print(dimsizes(lon))
; printMinMax(lon,True) ; prints the minimum and maximum longitude (0.3125 ... 359.6875)
; print(lat) ; prints all 360 latitude entries
print("Finished: 2/8 Print the dimensions of the read-in file")
; ***********************************************
; READ IN THE VARIABLE (precipitation)
; ***********************************************
precip = nmme_f[:]->pr*3600*24 ; extract precipitation into new variable, output is mm/day
; print("variable dimensions")
; print(dimsizes(precip))
; num_precip = num(precip) ; finds out the number of total values in all 12 model runs
; max_precip = max(precip) ; finds the smallest value in any of the 12 model runs
; print("The total number of values in the 12 1980 GFDL model runs are:")
; print(num_precip)
; print("Largest number in the GFDL 1980 files:")
; print(max_precip)
; since we used join command above, the dimensions of precip should be (number files,ndays,nlat,nlon)
print("Finished: 3/8 Read in the variable (precipitation)")
; ***********************************************
; CREATE US-SUBSET OF THE GLOBAL DATA
; ***********************************************
; let's choose latitudes of
latind = (/20,51/) ; latind is a variable
; let's choose longitudes of
lonind = (/230,306/) ; lonind is a variable
lti = ind_nearest_coord(latind,lat,0) ; lti, lni = variables, ind_nearest_coord = pre-defined functions to determine the indices of locations closest to the coordinate array, needed later to define precip_subdomain_4D
lni = ind_nearest_coord(lonind,lon,0)
lat_subdomain = lat(lti(0):lti(1))
lon_subdomain = lon(lni(0):lni(1))
precip_subdomain_3D = precip(:,lti(0):lti(1),lni(0):lni(1)) ; define the subdomain based on latind and lonind
; print(lat_subdomain_4D) ; displays the latitudes of the chosen area
; print(lon_subdomain_4D) ; displays the longitudes of the chosen area
; print(dimsizes(precip_subdomain_3D)) ; prints all lats and lons of the subset
day = nmme_f[:]->time
ensemble_member = ispan(1,12,1) ; there are 12 ensemble members
print("Finished: 4/8 Create US-subset of the global data")
; ***********************************************
; WRITE THE SUBSET INTO NEW netCDF FILE
; ***********************************************
x = (/year,100/)
y = product(x)
z = (/y,m/)
init = sum(z) ; filename element in YYYYMM format to indicate initiation month
netCDF1 = True ; Output format is NetCDF
if (netCDF1) then
diro = "./subset/nc/" ; output in subdirectory "subset" within input directory
filename = "GFDL_US_pr_daily_init_" + init + "run_" + r +".nc" ; output filename for subset using initialization year and month to differentiate
filo = filename
end if
setfileoption("nc","Format","LargeFile")
if (netCDF1) then
system("/bin/rm -f "+diro+filo) ; removes old files with the same name if they are present
ncdf = addfile(diro+filo,"c")
setfileoption(ncdf,"DefineMode",True)
; create attributes
fAtt = True
fAtt at title = "Data: Daily GFDL FLOR-B01, run: " + r + " i1p1, for the Contiguous U.S., initialized " + init
fAtt at source = "Vecchi et al. 2014"
fAtt at Conventions = "None"
fAtt at creation_date = systemfunc("date")
fileattdef(ncdf,fAtt)
dimNames = (/"year","ensemble_member","days","lat","lon"/)
dimSizes = (/-1,12,dimsizes(day),dimsizes(lat_subdomain),dimsizes(lon_subdomain)/) ; r = model run (1-12), was "12" before
dimUnlim = (/True,False,False,False,False/)
filedimdef(ncdf,dimNames,dimSizes,dimUnlim)
filevardef(ncdf,"year",typeof(year),(/"year"/))
filevardef(ncdf,"lat",typeof(lat_subdomain),(/"lat"/))
filevardef(ncdf,"lon",typeof(lon_subdomain),(/"lon"/))
; print("Variable day: ")
; Print(day)
filevardef(ncdf,"days",typeof(day),(/"days"/))
filevardef(ncdf,"ensemble_member",typeof(ensemble_member),(/"ensemble_member"/))
; get attributes of our variable
filevardef(ncdf,"precipitation",typeof(precip_subdomain_3D),(/"days","lat","lon"/))
; now write all the variables to the file
ncdf->year = (/year/)
ncdf->ensemble_member = (/ensemble_member/)
ncdf->days = (/day/)
ncdf->lat = (/lat_subdomain/)
ncdf->lon = (/lon_subdomain/)
ncdf->precipitation = (/precip_subdomain_3D/)
end if
print("Finished: 5/8 Write the subset into new netCDF file")
; ***********************************************
; REDUCE 3D ARRAY TO PLOTABLE 2D ARRAY
; ***********************************************
; print(dimsizes(precip_subdomain_3D)) ; prints all lats and lons of the subset
; printVarSummary(precip_subdomain_3D) ; print summary of subdomain precip: var name, type, byte size, # of values, dimensions (x,y,z,temp)
precip_subdomain_2D = dim_avg_n(precip_subdomain_3D, 0) ; averaging all daily precipitation values into one annual average for every lat and lon, creating 2D array, needed to plot single map
; printVarSummary(precip_subdomain_3D) ; print size (byte, values) and dimensions (x, y, z)
; print(precip_subdomain_3D) ; prints *ALL* precip values for the subset, about 5 million lines total!!!
; print(precip_subdomain_3D(:,:,:)) ; array subscripting for all subset lats and lons for 1980, all 12 model runs
print("Finished: 6/8 Reduce 3D array to plottable 2D array")
; ***********************************************
; CREATE MAP
; ***********************************************
; Define boundaries for map
precip_subdomain_2D!0 = "lat_subdomain" ; latitude information come from variable lat_subdomain
precip_subdomain_2D&lat_subdomain = lat_subdomain
precip_subdomain_2D!1 = "lon_subdomain" ; longitude information come from variable lon_subdomain
precip_subdomain_2D&lon_subdomain = lon_subdomain
precip_subdomain_2D at long_name = "Precipitation" ; map subtitle, printed on the top left
precip_subdomain_2D at units = "mm/day" ; map units, printed on the top right
diro = "./subset/png/" ; rename output directory to save map in png folder
wks = gsn_open_wks("png", diro + "GFDL_US_pr_daily_init_" + init + "run_" + r) ; map name and file type
cmap = read_colormap_file("BlueDarkRed18")
res = True ; plot mods desired
res at tiMainString = "Run " + r + ", init. " + m + "/" + year + " - GFDL-Hindcast, daily" ; main title
res at cnFillOn = True ; turn on color fill, works with gsn_csm_map() function 3 lines down
res at gsnAddCyclic = False
res at cnLinesOn = False ; contour lines of the map fill on or off
res at mpProjection = "Mercator" ; Mercator map projection
res at cnFillPalette = cmap(::-1,:) ; reverse color map
res at mpLimitMode = "LatLon" ; define plotted area by lats and lons, namely those of the selection (20-50N, 230-300E)
res at mpMaxLatF = 50
res at mpMinLatF = 20
; res at mpCenterLatF = 35 ; define center latitude, only use with res at mpPrjection is conical equidistant
res at mpMaxLonF = 305
res at mpMinLonF = 230
; res at mpCenterLonF = 265 ; define center longitude, only use with res at mpPrjection is conical equidistant
res at mpOutlineBoundarySets = "GeophysicalAndUSStates" ; lines and boundaries of USA and US States
res at mpGeophysicalLineThicknessF = 2 ; thickness of the USA lines
res at mpUSStateLineThicknessF = 2 ; thickness of the state lines
res at mpSpecifiedFillColors = (/0,100/) ; fill with background color, works with gsn_map() 2 lines down
print("Finished: 7/8 Create map")
; ***********************************************
; PLOT MAP
; ***********************************************
; plot = gsn_csm_contour_map_ce(wks,precip_subdomain_2D, res) ; conical equidistant projection, if used take out res at mpProjection line
plot = gsn_csm_contour_map(wks,precip_subdomain_2D, res) ; universal map projection, define in res at mpProjection
print("Finished: 8/8 Plot map")
print("**** ALL DONE: " + init + ", run " + r + " ****")
delete(precip) ; delete before the next loop or program will fail
delete(precip_subdomain_3D) ; delete before the next loop or program will fail
; delete(precip_subdomain_2D) ; delete before the next loop or program will fail
delete(day) ; delete before the next loop or program will fail
end do ; end the loop through the 12 model runs (1-12)
end do ; end the loop through the 12 calendar months (1-12)
end do ; end the loop through the years of 1980-2013
end
> On Jan 29, 2019, at 11:50 AM, prithvi raj <rajjallu at gmail.com <mailto:rajjallu at gmail.com>> wrote:
>
> Hello all,
> I am rather new to NCL. I have cyclone track data from 1945-2018 as ascii files. In those files the date, lat, lon, wind speed at that particular lat and lon is mentioned. I want to segregate the them based on the intensity of cyclones into severe and very severe and plot them as bars in a time series. Can it be done using NCL?
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu <mailto: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/20190129/b4eef052/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: GFDL_subset_save_and_plot_DAILY.ncl
Type: application/octet-stream
Size: 10524 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190129/b4eef052/attachment.obj>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190129/b4eef052/attachment-0001.html>
More information about the ncl-talk
mailing list