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/wrf/WRFUserARW.ncl" load "./map_tickmarks.ncl" ;=============================================================; ; Example of script to plot ZS field from MESONH netcdf4 file ; Could be use with any files generated by MESONH if LCARTESIAN=False ;=============================================================; begin ;=============================================================; ; Open file ;=============================================================; mnh_file="PEARL0915_PGD_10km.nc" a = addfile(mnh_file, "r") ;=================================================; ; Get information on variable sizes ;=================================================; jphext = a->JPHEXT mdims = getfilevardimsizes(a,"ZS") ; get some dimension sizes for the file nd = dimsizes(mdims) imax=mdims(nd-1)-2*jphext ; -2*jphext : to remove non-physical values jmax=mdims(nd-2)-2*jphext ; -2*jphext : to remove non-physical values ;=================================================; ; Read the variables ;=================================================; ; Note: do not read first and last value which are non physical values ; ------------------------------------- lat=a->LAT(jphext:jmax+jphext-1,jphext:imax+jphext-1) lon=a->LON(jphext:jmax+jphext-1,jphext:imax+jphext-1) ZS = a->ZS(jphext:jmax+jphext-1,jphext:imax+jphext-1) ZS@long_name="Orography" ZS@units="Unit: metre" ; Associate coordinates to variable ; Necessary to plot the field at the right geographical position ; ------------------------------------- ZS@lat2d = lat ZS@lon2d = lon ; Read projection parameters ; ------------------------------------- RPK = a->RPK BETA= a->BETA LON0= a->LON0 LAT0= a->LAT0 ;=================================================; ; CREATE PLOT : Orography (ZS, filled contour) ;=================================================; ; Open workstation and define colormap ; ------------------------------------- type = "png" ; open a x11 window ; change type to png, ps or pdf to get the plot into a file wks = gsn_open_wks(type,"plt_BasicMap") ; gsn_define_colormap(wks,"topo_15lev") ; Choose colormap ;========================================== ; Map ressources ;========================================== resmap = True resmap@gsnDraw = False resmap@gsnFrame = False ; If there is an error on HighRes, it means that you don't have the HighRes data ; You need to download them or the change HighRes by MediumRes ; See https://www.ncl.ucar.edu/Document/Graphics/rangs.shtml for info ; ------------------------------------- resmap@mpDataBaseVersion = "HighRes" ; choose highres map data version (must be donwloaded) resmap@mpGridAndLimbOn = True ; turn on lat/lon lines ; resmap@mpGridLatSpacingF = 2. ; spacing for lat lines ; resmap@mpGridLonSpacingF = 10. ; spacing for lon lines resmap@pmTickMarkDisplayMode = "Always" ; turn on tickmarks ;=================================================; ; Set map projection ressources using projection parameters ;=================================================; if (RPK.gt.0) ; --------------------------- ; Lambert projection from north pole ; --------------------------- resmap@mpProjection = "LambertConformal" ; projection pole = 1 ; projection for north hemisphere resmap@mpLambertParallel1F = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere resmap@mpLambertParallel2F = resmap@mpLambertParallel1F ; ncl adds from grib file resmap@mpLambertMeridianF = LON0 ; ncl adds from grib file end if if (RPK.lt.0) ; --------------------------- ; Lambert projection from south pole ; --------------------------- resmap@mpProjection = "LambertConformal" ; projection pole = -1 ; projection for south hemisphere resmap@mpLambertParallel1F = pole*asin(RPK)*180/3.14 ; with pole=1 for north hemisphere and -1 for south hemisphere resmap@mpLambertParallel2F = resmap@mpLambertParallel1F ; ncl adds from grib file resmap@mpLambertMeridianF = LON0 ; ncl adds from grib file end if if (RPK.eq.1) ; --------------------------- ; Stereographic projection north ; --------------------------- resmap@mpProjection = "Stereographic" resmap@mpRelativeCenterLon = True resmap@mpRelativeCenterLat = True resmap@mpCenterLonF = LON0 resmap@mpCenterRotF = BETA ; resmap@mpCenterLatF = LAT0 resmap@mpCenterLatF = 90. end if if (RPK.eq.-1) ; --------------------------- ; Stereographic projection south ; --------------------------- resmap@mpProjection = "Stereographic" resmap@mpCenterLonF = LON0 resmap@mpCenterRotF = BETA resmap@mpCenterLatF = -90 end if if (RPK.eq.0) then ; --------------------------- ; Mercator projection ; --------------------------- resmap@mpProjection = "Mercator" end if print("Map projection="+resmap@mpProjection) ;========================================== ; Defining the corners for projection otherwise the plot will be global ;========================================== resmap@mpLimitMode = "Corners" resmap@mpLeftCornerLatF = lat(0,0) resmap@mpLeftCornerLonF = lon(0,0) resmap@mpRightCornerLatF = lat(jmax-1,imax-1) resmap@mpRightCornerLonF = lon(jmax-1,imax-1) print ("Corner (0,0); Lat="+resmap@mpLeftCornerLatF+ \ ", Lon="+resmap@mpLeftCornerLonF) print ("Oppos corner; Lat="+resmap@mpRightCornerLatF+ \ ", Lon= "+resmap@mpRightCornerLonF) ;========================================== ; Create ZS plot (contour) ;========================================== ; --------------------------- ; General ressources ; --------------------------- resmap@gsnMaximize = True ; maximize size resmap@gsnSpreadColors = True ; use full range of colormap ; --------------------------- ; Contour ressources ; --------------------------- resmap@cnFillOn = True ; turn on color fill resmap@cnLinesOn = True ; turn off contour lines resmap@cnLevelSelectionMode = "ManualLevels"; Manually set contour levels. resmap@cnMinLevelValF = 5. ; resmap@cnMaxLevelValF = 3200. ; resmap@cnLevelSpacingF = 200. ;-------------------------------; ; X-axis title (tiY) ;-------------------------------; resmap@tiXAxisFontHeightF = 0.018 ; font height resmap@tiXAxisFont = 21 ; font index resmap@tiXAxisString = "longitude" ; string to use as the X-Axis title ;-------------------------------; ; Y-axis title (tiY) ;-------------------------------; resmap@tiYAxisFontHeightF = 0.018 ; font height resmap@tiYAxisFont = 21 ; font index resmap@tiYAxisString = "latitude" ; string to use as the Y-Axis title resmap@lbOrientation = "Vertical" resmap@pmLabelBarOrthogonalPosF= 0.09 ; move labelbar away from tickmarks resmap@tiDeltaF = 5.0 ; move X/Y axis strings away from tickmarks resmap@gsnRightStringOrthogonalPosF = 0.08 ; Move the subtitles up a bit resmap@gsnLeftStringOrthogonalPosF = 0.08 ;=================================================; ; PLOT ZS ;=================================================; ; resmap@tmYLOn = False ; resmap@tmYROn = False ; resmap@tmYUseLeft = False ; test = lat(:,0) ; resmap@tmYLMode = "Explicit" ; resmap@tmYRValues = test plot_zs=gsn_csm_contour_map(wks,ZS,resmap) ;---We're putting tickmarks on all four axes. tmres = True tmres@tmYLValues = ispan(70,82,1) tmres@tmXBValues = ispan(-101,-52,3) tmres@tmYRValues = ispan(70,82,1) tmres@tmXTValues = ispan(-150,-50,5) tmres@tmXBLabelFontHeightF = 0.01 tmres@tmYLLabelFontHeightF = 0.01 ;---Attach the new map tickmarks. plot_zs = add_map_tickmarks2(wks,plot_zs,tmres) coordlateu = 79.5920 coordloneu =-85.5627 polyres = True polyres@gsMarkerIndex = 16 ; polymarker style polyres@gsMarkerColor = "red" ; Color of marker polyres@gsMarkerSizeF = 8. ; polymarker size polyres@gsMarkerOpacityF = 1. plottest = gsn_add_polymarker(wks, plot_zs, coordloneu, coordlateu, polyres) txres = True txres@txJust = "TopLeft" txres@txFontHeightF = 0.012 txid = gsn_add_text(wks,plot_zs," "+"Station Eureka", coordloneu, coordlateu, txres) draw(plot_zs) frame(wks) end