;---------------------------------------------------------------------- ; This example shows how to read geographic data ; from Natural Earth shapefiles ; and plot them as polylines and polygons. ;---------------------------------------------------------------------- ; This particular example plots data for Switzerland. ;---------------------------------------------------------------------- ; Download the shapefiles from http://www.naturalearthdata.com/ ; Unzip to a directory ;---------------------------------------------------------------------- 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" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl" begin ;***************************************** ; Open workstation and define colormap * ;***************************************** wks_type = "png" ; use "newpdf" instead of "pdf" for smaller files wks_type@wkWidth = 1200 wks_type@wkHeight = 1200 wks = gsn_open_wks(wks_type,"temp") gsn_define_colormap(wks,"rainbow+white+gray") ; choose colormap ;****************** ; read gfs data * ;****************** systemdate = systemfunc("date +%H%M") if ((systemdate.eq."0645"))then url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" + systemfunc("date +%Y%m%d") + "/" filename = url + "gfs_hd_00z" end if if ((systemdate.eq."1245")) then url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" + systemfunc("date +%Y%m%d") + "/" filename = url + "gfs_hd_06z" end if if ((systemdate.eq."1845")) then url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" + systemfunc("date +%Y%m%d") + "/" filename = url + "gfs_hd_12z" end if if ((systemdate.eq."0045")) then url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd" + systemfunc("date -d '-1 day' +'%Y%m%d'") + "/" filename = url + "gfs_hd_18z" end if url = "http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20111001/" filename = url + "gfs_hd_00z" exists = isfilepresent(filename) if(.not.exists) then print("OPeNDAP isfilepresent test unsuccessful.") print("Either file doesn't exist, or NCL does not have OPeNDAP capabilities on this system") else print("OPeNDAP isfilepresent test successful.") gfs = addfile(filename,"r") vars = getfilevarnames(gfs) end if ;********************************* ; risorse per la mappa di base * ;********************************* mappa = True ; plot mods desired mappa@mpProjection = "LambertConformal" ; choose projection mappa = True ; plot mods desired mappa@gsnMaximize = True mappa@gsnDraw = False ; do not draw the plot mappa@gsnFrame = False ; do not flip the page mappa@mpProjection = "LambertConformal" ; choose projection mappa@mpLimitMode = "LatLon" mappa@mpMaxLonF = 47.0 mappa@mpMaxLatF = 88.5 mappa@mpMinLatF = 27.0 mappa@mpMinLonF = -43.0 mappa@mpGridAndLimbOn = True ; turn on grid lines mappa@mpGridLineDashPattern = 10 ; lat/lon lines dashed ;****************************** ; risorse generali dei plot * ;****************************** plot_generale= True plot_generale@gsnMaximize = True plot_generale@gsnDraw = False ; do not draw the plot plot_generale@gsnFrame = False ; do not flip the page plot_generale@gsnAddCyclic = True plot_generale@gsnSpreadColors = False ; turn off the use entire color map plot_generale@gsnRightString = " " ; turn off right string plot_generale@gsnLeftString = " " ;******************************* ; risorse generali dei testi * ;******************************* testo_generale = True testo_generale@txFont = "times-roman" ; font testo_generale@txBackgroundFillColor = "Transparent" drawNDCGrid(wks) ;****************************** ; generazione mappa di base * ;****************************** map = gsn_csm_map(wks,mappa) fnames = "/mnt/internetserver/map/shapefile/europa/aaa_full/" + (/"europa"/) + ".shp" ; Files to be open linecolors = (/"Black"/) ; color definition for each shapefile fillcolors = (/"Transparent"/) ; Fill color definition for each shapefile thicks = (/2/) ; Thickness for each shapefile lnres = True ; resources for polylines plres = True prims = True lines = True do n=0,dimsizes(fnames)-1 ; Loop through files that we want to read geographic information from. f = addfile(fnames(n),"r") ; Open the shapefile number n. segments = f->segments ; Read data off the shapefile geometry = f->geometry segsDims = dimsizes(segments) geomDims = dimsizes(geometry) geom_segIndex = f@geom_segIndex ; Read global attributes geom_numSegs = f@geom_numSegs segs_xyzIndex = f@segs_xyzIndex segs_numPnts = f@segs_numPnts geometry_type = f@geometry_type numFeatures = geomDims(0) lon = f->x lat = f->y if (geometry_type.eq."polygon") then ; Put if statement outside the loop plres@gsFillColor = fillcolors(n) plres@gsEdgesOn = True ; Draw border around polygons plres@gsEdgeColor = linecolors(n) plres@gsEdgeThicknessF = thicks(n) do i=0, numFeatures-1 ; Section to draw polygons on map. startSegment = geometry(i, geom_segIndex) numSegments = geometry(i, geom_numSegs) do seg=startSegment, startSegment+numSegments-1 startPT = segments(seg, segs_xyzIndex) endPT = startPT + segments(seg, segs_numPnts) - 1 dumstr = unique_string("lines") ; This call adds the polygon. map@$dumstr$ = gsn_add_polygon(wks, map , lon(startPT:endPT), lat(startPT:endPT), plres) end do end do else lnres@gsLineThicknessF = thicks(n) lnres@gsLineColor = linecolors(n) do i=0, numFeatures-1 ; Section to draw polylines on map. startSegment = geometry(i, geom_segIndex) numSegments = geometry(i, geom_numSegs) do seg=startSegment, startSegment+numSegments-1 startPT = segments(seg, segs_xyzIndex) endPT = startPT + segments(seg, segs_numPnts) - 1 dumstr = unique_string("primitive") ; This call adds the line segment. map@$dumstr$ = gsn_add_polyline(wks, map, lon(startPT:endPT), lat(startPT:endPT), lnres) end do end do end if delete(lat) ; Clean up before we read in same variables again. delete(lon) delete(segments) delete(geometry) delete(segsDims) delete(geomDims) end do ;******************************* ; 2 metres temperature plot * ;******************************* temperatura_a_due_metri = plot_generale testo1_temperatura_a_due_metri = testo_generale testo2_temperatura_a_due_metri = testo_generale testo3_temperatura_a_due_metri = testo_generale testo4_temperatura_a_due_metri = testo_generale ;testo1_temperatura_a_due_metri@txFontHeightF = 0.08 ; set the label size ;testo1_temperatura_a_due_metri@txJust = "CenterLeft" ;testo1_temperatura_a_due_metri@txPosXF = 0.3 ; Rough approximation ;testo1_temperatura_a_due_metri@txPosYF = 0.9 testo2_temperatura_a_due_metri@txFontHeightF = 0.015 testo2_temperatura_a_due_metri@txJust = "BottomRight" testo2_temperatura_a_due_metri@txPosXF = 0.98 ; Rough approximation testo2_temperatura_a_due_metri@txPosYF = 0.85 temperatura_a_due_metri@cnLineDashSegLenF = 0.04 temperatura_a_due_metri@cnLinesOn = True ; turn on contour lines temperatura_a_due_metri@cnLineLabelsOn = True ; turn on contour labels temperatura_a_due_metri@cnLineLabelFont = "times-roman" ; set font on contour labels temperatura_a_due_metri@cnLineLabelFontHeightF = 0.006 ; set font height on contour labels temperatura_a_due_metri@cnLineLabelPlacementMode = "Constant" ; set placement mode of labels on contour temperatura_a_due_metri@cnLineLabelBackgroundColor = "Transparent" ; set background color of the contour label temperatura_a_due_metri@cnLevelSelectionMode = "ManualLevels" ; set how manage contour line levels temperatura_a_due_metri@cnMinLevelValF = -50 ; max level do draw contour lines temperatura_a_due_metri@cnMaxLevelValF = 50 ; min level do draw contour lines temperatura_a_due_metri@cnLevelSpacingF = 1 ; set spacing between contour lines temperatura_a_due_metri@cnLevelFlags = (/"LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly","LineAndLabel","LineOnly","LineOnly","LineOnly","LineOnly", \ "LineAndLabel"/) ; set how draw contour lines temperatura_a_due_metri@cnFillOn = True ; turn on color temperatura_a_due_metri@cnFillColors = (/12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78,80,82,84,86,88,90,92,94,96,98,100,102,104,106,108,110,112,114,116,118,120,122,124,128,132,136,140,144,148,152,156,160,164,166,168,170,172,174,176,178,180,182,184,186,188,190,192,194,196,198,200,202,204,206,208,210,212,214,216,218,220,222,224,226,228,230,232/) temperatura_a_due_metri@cnMonoLineColor = False ; turn off auto contour line colors temperatura_a_due_metri@cnLineColors = (/"Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","White","White","Black","White","White","White","White", \ "Black","White","White","white","White","Black","White","White","White","White", \ "Black"/) ; set contour lines colors temperatura_a_due_metri@cnMonoLineThickness = False ; turn off auto contour line thickness temperatura_a_due_metri@cnLineThicknesses = (/1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0,0.5,0.5,0.5,0.5,1.0,0.5,0.5,0.5,0.5, \ 1.0/) ; set contour lines thickness ; Create a labelbar object. temperatura_a_due_metri@lbLabelsOn = True ; turn on the labels temperatura_a_due_metri@lbLabelAutoStride = True ; let NCL determine label spacing temperatura_a_due_metri@lbLeftMarginF = 0.02 temperatura_a_due_metri@lbOrientation = "horizontal" ; label orientation temperatura_a_due_metri@lbBoxLinesOn = True ; turn on lines between labelbar colors temperatura_a_due_metri@lbLabelFont= "times-roman" ; font temperatura_a_due_metri@lbLabelFontHeightF= 0.0045 ; set the label size temperatura_a_due_metri@lbLabelAlignment = "BoxCenters" temperatura_a_due_metri@pmLabelBarHeightF = 0.04 temperatura_a_due_metri@lbLabelStride = 2 ; label every other box ;temperatura_a_due_metri@lbMonoFillColor = False temperatura_a_due_metri@lbBoxCount = 101 temperatura_a_due_metri@lbBoxLinesOn = True ;temperatura_a_due_metri@lbFillColors = (/12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78,80,82,84,86,88,90,92,94,96,98,100,102,104,106,108,110,112,114,116,118,120,122,124,128,132,136,140,144,148,152,156,160,164,166,168,170,172,174,176,178,180,182,184,186,188,190,192,194,196,198,200,202,204,206,208,210,212,214,216,218,220,222,224,226,228,230,232/) ;**************************** ; risorse per la label bar * ;**************************** ;---data for 2 metres temperature plot TEMP2M = gfs->tmp2m(:,:,:) ; temperature NTIMES = dimsizes(gfs->time) ; number of times in the file TEMP2M = TEMP2M - 273.15 ;TEMP2M@units = "(C)" ;TEMP2M@long_name = "Temperatura a 2 metri dal suolo" ;testo1_temperatura_a_due_metri_plot = gsn_create_text(wks,"Temperatura a 2 metri dal suolo",testo1_temperatura_a_due_metri) do it = 0,10 ;NTIMES-1 temperatura_a_due_metri_plot = gsn_csm_contour(wks,TEMP2M(0,:,:),temperatura_a_due_metri) overlay(map,temperatura_a_due_metri_plot) draw(map) format ="%Y%N%D %H" date = ut_string(gfs->time(it),format) print(date) date1 = systemfunc("date -d '20111002 00' +'%A %H'") date2 = systemfunc("date -d '(date)' +'%A %H'") print(date1) print(date2) testo2_temperatura_a_due_metri_plot = gsn_create_text(wks,systemfunc("date -d '(ut_string(gfs->time(it),format))' +'%A %H'"),testo2_temperatura_a_due_metri) draw(testo2_temperatura_a_due_metri_plot) ; Draw text box. frame(wks) end do do it = 0,10 ;NTIMES-1 if (it.le.8) then system("mv temp.00000" + (it+1) + ".png" + " temperatura0" + (it) + "_" + systemfunc("date +%Y%m%d%H") + ".png") end if if (it.eq.9) then system("mv temp.0000" + (it+1) + ".png" + " temperatura0" + (it) + "_" + systemfunc("date +%Y%m%d%H") + ".png") end if if (it.ge.10) then system("mv temp.0000" + (it+1) + ".png" + " temperatura" + (it) + "_" + systemfunc("date +%Y%m%d%H") + ".png") end if end do end