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/contrib/cd_string.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "/al11/andrea/research/Elliptical/EllipticalFit.ncl" ;For the dates set below, this code finds the 30km hght contour at 10-hPa and feeds the lat/lon points of the contour ;to a python scripts to calculate the best fit ellipse of that 30 km contour line. (Needed python for imaginary number work) ;This script then reads in the output from that python script (which is in a polar sterographic coord system) ;and converts it into a plotable form to get the lat&lon center, the axes lengths and the angle of rotation of the best fit ellipse ;This data is then written to the text file EMetricsYYYY1YYYY2 ;For 1996: May 1: 122; Jun 1: 153; July 1: 183; Aug 1: 214; Sept 1: 245; Oct 1: 275; Nov 1: 306; Dec 1:336 begin lev2plot = 10 lowerlev = 1000 year1 = 1996 year1@calendar = "366_day" ;change to "365_day" for non leap years mm1 = 11 dd1 = 01 ;nens = 11 ;nff = 32 print("y: "+year1+" "+mm1+" "+dd1) ;forecast intiialization date datestr = "19961101" ; "1996111100" ;for GFSre ;; systemfunc("date -u '+%y%m%d'") ; for realtime date modellist = (/"CMA"/);"HMCR"/) ; the list of models to use modelcolors = (/"forestgreen","orangered","deepskyblue3"/) ; colors for each models plots = new(3,graphic) nmodels = dimsizes(modellist) ;filename = "ellipseS2S_"+lev2plot+"_65_BOMCMAMetFr"+year1+mm1+dd1 maxff = 1488 ; max forecast hr of all models included modelnff = new(nmodels,integer,0) modelnens = new(nmodels,integer,0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; loop through modellist ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do j=0,nmodels-1 modelname = modellist(j) datestr = "19961101" ; "1996111100" ;for GFSre ;; systemfunc("date -u '+%y%m%d'") ; for realtime date fstart = 0 if (modelname.eq."BOM") then fstart = 1 end if ; if first forecast is F00, set fstart to 0, if the first forecast in the data is F24, set fstart to 1 ... ; CMA -> F0 start ; BOM -> F24 start ; ECMWF -> F0 start ; MetFr -> F0 start ; HMCR -> F0 start ;;;open from S2S files ; I've saved my S2S data in the format yyyymmddhh_modelname.nc g_file := addfile("/langlab_rit/andrea/Daledata/S2Sdata/"+datestr+"_"+modelname+".nc","r") ; Dale data is in langlab_rit ;;grab the data ; print(getfilevarnames(g_file)) ;;; shows content of g_file ginfo := g_file->gh printVarSummary(ginfo) garray := g_file->gh(:,:,{lev2plot},{0:90},:) ; var names checked against the print statments above uarray := g_file->u(:,:,{lev2plot},{65},:) ; (forecast/valid time, ens member, level, lat, lon) tarray := g_file->t(:,:,{lev2plot},{75:90},:) ens := g_file->number(:) ; array of ens member numbers ftimes := g_file->time(:) ; array of forecast times ; units : hours since 1900-01-01 00:00:0.0 ; calendar : gregorian ; printVarSummary(garray) lat := tarray&latitude lon := garray&longitude ylat := dimsizes(lat) xlon := dimsizes(lon) nens := dimsizes(ens) nff := dimsizes(ftimes) ;nff := 63 print("N: "+nens+" F: "+nff) modelnff(j) = nff modelnens(j) = nens ; nens = 11 ;--> 11 ens members in GFS Reforecast ; nff = 32 ;--> 32 forecast times (192hr forecast) if (modelname.eq."BOM") then nff = 58 end if print("y: "+year1+" "+mm1+" "+dd1) ;day1 = day_of_year(year1,mm1,dd1) date1 = cd_inv_calendar( tointeger(year1), tointeger(mm1), tointeger(dd1), 00, 00, 00, "hours since 1800-01-01 00:00", 0 ) ;change date format --> S2S data in hrs since 1900, cfsr in hrs since 1800 ;date2 = cd_inv_calendar( year1, mm2, dd2, 18, 00, 00, "hours since 1900-01-01 00:00", 0 ) ;calc zonal means and polar cap avgs ; gfcast := dim_avg_wgt_Wrap(dim_avg_Wrap(garray),cos((3.1415/180)*lat),0) ; calc zonal mean then polar cap lat-wgt'd avg ufcast := dim_avg_Wrap(uarray) ; calc zonal mean (time, ens,lat lon) tfcast := dim_avg_wgt_Wrap(dim_avg_Wrap(tarray),cos((3.1415/180)*lat),0) ; calc zonal mean then polar cap lat-wgt'd avg ;calc ens means ensmean_g := dim_avg_n_Wrap(garray,1) ;;; pretty up the forecast time strings figstart = 0 nmtrc = (nff) * (nens+1) metrics = new(nmtrc,string) ;;;;;;;;;;;Convert times to readable dates;;;;;;;;;;; ; Array to hold month abbreviations. Don't store anything in index ; '0' (i.e. let index 1=Jan, 2=Feb, ..., index 12=Dec). ; month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \ "Oct","Nov","Dec"/) ; Convert to UTC time. utc_date = cd_calendar(ftimes, 0) ; Store return information into more meaningful variables. ; yy = tointeger(utc_date(:,0)) ; Convert to integer for mm = tointeger(utc_date(:,1)) ; use sprinti dd = tointeger(utc_date(:,2)) hh = tointeger(utc_date(:,3)) ;dayn = day_of_year(yy,mm,dd) yystr = sprinti("%0.4i", yy) mmstr = month_abbr(mm) ddstr = sprinti("%0.2i", dd) hhstr = sprinti("%0.2i", hh) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;Begin looping over times;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do n=0,nff-1 ;forecast and analysis times fignum = figstart + n print(n) fhr = (n+fstart)*24 ; forecast times * 24 gives us forecast hour for S2S data ;BOM data starts with f24 fff = sprinti("%0.4i",fhr) ; string of forecast hour for titles/filenames if (fhr.lt.1000) then ; check forecast hr to set string correctly for ff time fff = sprinti("%0.3i",fhr) ; gfs reforecast ens data has forecast hours names e.g., 06 09 12 60 108 if (fhr.lt.100) then ; check forecast hr to set string correctly for ff time fff = sprinti("%0.2i",fhr) ; gfs reforecast ens data has forecast hours names e.g., 06 09 12 60 108 end if end if datestr := "F"+fff+" - "+hhstr(n)+" UTC "+ddstr(n)+" "+mmstr(n)+" "+yystr(n) ensmean := ensmean_g(n,:,:) printVarSummary(ensmean) do mbr=0,nens ; extra for ensmean cntr = new(nens+1,graphic) print("F: "+fff+" Mbr: "+mbr) if(mbr.eq.0) then ;--------------start with ens mean run-------------- ;************************************************ ; create plot ;************************************************ wks_type = "png" wks_type@wkWidth = 800 wks_type@wkHeight = 600 filename = "ellipse"+modelname+"_"+n ;filename = "Test" wks1 = gsn_open_wks("png","tmp") ; open a png file wks = gsn_open_wks(wks_type ,filename) gsn_define_colormap(wks, "GMT_polar") ; choose colormap res = True res@gsnMaximize = True ;make that map bigger res@gsnDraw = False ; Don't draw plot res@gsnFrame = False ; Don't advance frame res@vpHeightF = 0.9 ; change aspect ratio of plot res@vpWidthF = 0.9 res@gsnRightString = "" ; turn off right string res@gsnLeftString = "" ; turn off left string res@cnFillOn = False ; turns on the color res@mpFillOn = False ; turns off continent gray res@cnLinesOn = True ; turn off contour lines cnres = res ;make contour specific resource list plyres = res ;make polyline specific resource list fres = res res@mpGeophysicalLineColor = "gray60" res@mpGeophysicalLineThicknessF = 2. res@mpProjection = "Orthographic" ; choose projection res@mpPerimOn = False ; turn off box around plot res@mpFillOn = False ; fres@mpMonoFillColor = True ; fres@mpFillColor = "gray40" res@mpOutlineOn = True res@mpOutlineDrawOrder = "PostDraw" res@mpLimitMode = "LatLon" res@mpMinLatF = 10. res@mpMaxLatF = 90. res@mpCenterLonF = 180. ; choose center lon res@mpCenterLatF = 90. ; choose center lat res@mpGridAndLimbOn = True ; turn on lat/lon lines res@mpGridSpacingF = 10. res@mpGridLineColor = "gray75" res@mpFillOn = True res@mpLandFillColor = "gray75" fres = res fres@tiMainString = datestr;s" ";lev2plot+"-hPa 30km contour and Best-Fit Ellipse" ; plot title fres@tiMainFontHeightF = 0.021 ; fres@gsnCenterString = datestr ; center title ; fres@gsnCenterStringFontHeightF = 0.019 ; fres@cnFillOn = False ; color plot desired ; fres@cnLineLabelsOn = False ; turn off contour line labels ; fres@cnLineLabelFontHeightF = 0.005 ; fres@cnLinesOn = True ; turn off contour lines ; fres@gsnSpreadColors = True ; use full range of color map ; fres@cnLevelSelectionMode = "ExplicitLevels" ; fres@cnLevelSelectionMode = "ManualLevels" ; fres@cnInfoLabelOn = False ; fres@cnLevelSpacingF = .25 ; spacing ; fres@cnMinLevelValF = -3.5 ; min value ; fres@cnMaxLevelValF = 3.5 ; max value ; fres@cnLevels = (/-2,-1.5,-1.,-0.5,0.5,1,1.5,2/) ; fres@cnMonoLineColor = False ; fres@cnLineThicknessF = 1.5 ; fres@cnLineColors = (/"blue","blue","blue","blue","red","red","red","red"/) ; fres@gsnLeftString = "" ; fres@gsnRightString = "" ; fres@gsnContourNegLineDashPattern = 11 ; fres@gsnContourZeroLineThicknessF = 0.0 ; no zero line print("fills") fills = gsn_csm_map(wks,fres) ; create the contour plot cnres@cnLevelSelectionMode = "ExplicitLevels" ; set manual/explicit contour levels? cnres@cnLevels = (/28500., 29000.,29500.,30000.,30500.,31000./) ; cnres@cnLevels = (/21500.,22000.,22500.,23000.,23500.,24000./) ; cnres@cnLevels = (/18500.,19000.,19500.,20000.,20500.,21000./) cnres@cnMonoLineThickness = True cnres@gsLineThicknessesF = (/2.,2.,2., 4., 2., 2./) cnres@gsLineThicknessF = 2. cnres@cnLineThicknessF = 3 cnres@gsLineColor = "darkblue" cnres@cnLineLabelFontHeightF = 0.00 cnres@cnFillDrawOrder = "Predraw" cnres@cnInfoLabelString = "" cnres@cnLineLabelBackgroundColor = "transparent" ;the line break at the contour label cnres@cnInfoLabelOn = False cnres@cnLineLabelBackgroundColor = "transparent" cnres2=cnres cntr(mbr) = gsn_csm_contour(wks1,ensmean,cnres2) ; create the plot print("plotted cnres") ; info about contourline isoline := get_isolines(cntr(mbr),30000.) ;30000 for 10hpa, 23000 for 30hpa, 20000 move for 50 cntr(mbr) = gsn_csm_contour(wks,ensmean,cnres) ; create the plot overlay(fills,cntr(mbr)) plyres@cnLevelSelectionMode = "ManualLevels" plyres@gsLineThicknessF = 6. count = 0 segments := isoline@segment_count if (segments .gt. 2) then ;if there are more than 2 its probably an error segments := 2 ;so only do 2 end if ;IF USED FOR LEVELS .GT. 50 THIS PROBABLY NEEDS TO BE CHANGED error = 0 ;reset error to 0 every time plyres@gsLineColor = "mediumblue" ; change to blue every time ylat := isoline(0,:) ; b:e) <--------------- let try reading all points in xlon := isoline(1,:) ; b:e) <--------------- and fitting ellipse to all points if (dimsizes(ylat) .le. 50) then ;it its an error set error = 1 error = 1 end if ;;;;;;;remove comments of lines with i for doing segemnt things plyplot1 = gsn_add_polyline(wks,fills,xlon,ylat,plyres) ; count = count + isoline@n_points(i) ; print(isoline@level + " has " + count + " total points in " + isoline@segment_count + " segments" ) max_npts = max(isoline@n_points) lats := ylat*(3.1415/180) lons := xlon*(3.1415/180) ;convert lat and lon to new coords with NP at origin and 0 lon as y ; Following Waugh (1997) ;;;;;;;Write contour lat/lon to file;;;;;;;;;;;;;;; x := (cos(lons) * cos(lats))/(1 + sin(lats)) asciiwrite("x.txt",x) y := (sin(lons) * cos(lats))/(1 + sin(lats)) asciiwrite("y.txt",y) ;;;;;;;;Run Python code to find Ellipse;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;uses x.txt and y.txt to find ellipse then makes xx.txt and yy.txt;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; cmd = "python fitEllipse2.py" system(cmd) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print("Ran Python") xx := asciiread("xx.txt",-1,"float") yy := asciiread("yy.txt",-1,"float") ;convert polar sterographic cartesian coords to lat/lon xxlon := atan(yy/xx) xxlon := where(xx.lt.0, where(yy.gt.0, atan(yy/xx) +3.1415 ,atan(yy/xx)-3.14159), atan(yy/xx)) yysinxxlon := yy / sin(xxlon) pi4 = 3.1415 / 4 yylat := -2 * ( atan(yysinxxlon) - pi4) lats := yylat*(180/3.1415) lons := xxlon*(180/3.1415) do l=1,dimsizes(lats)-1 if(abs(lats(l)-lats(l-1)).gt.1.5)then lats(l) = lats(l-1) end if end do ;print(lats) ;print(lons) ;read in major and minor axis lengths in cart. coords, center of ellipse, and rotation of major axis in radians vals := asciiread("axescenterphi.txt", -1,"float") print(vals) aax := vals(0) bax := vals(1) centerx := vals(2) centery := vals(3) phi := vals(4) ;convert cartisian to lat/lon coords ; if x = ( cos(lon)cos(lat) ) / (1 + sin(lat)) and y = ( sin(lon)cos(lat) ) / (1 + sin(lat)) ; then lon = atan(y/x) and lat = 2[atan( y/sin(lon) ) - pi/4] ; via tangent half angle formula (see http://en.wikipedia.org/wiki/Tangent_half-angle_formula ) cenlon := atan(centery / centerx) cenlon := where(centerx.lt.0, where(centery.gt.0, atan(centery/centerx) +3.1415 ,atan(centery/centerx)-3.14159), atan(centery/centerx)) cenlondeg := cenlon*180/3.1415 ysinlon := centery / sin(cenlon) cenlat := -2 * ( atan(ysinlon) - pi4) cenlatdeg := cenlat*(180/3.1415) print("center of ellipse = "+cenlatdeg+" N, "+cenlondeg+" E.") ;Calc enpoints of the axes of the vortex xa = (aax) * cos(phi) ya = (aax) * sin(phi) xb = (bax) * sin(phi) yb = (bax) * cos(phi) endx = (/ centerx+xa , centerx-xa , centerx+xb , centerx-xb/) endy = (/ centery+ya , centery-ya , centery-yb , centery+yb/) endlon = atan(endy / endx) endlon = where(endx.lt.0, where(endy.gt.0, atan(endy/endx) +3.1415 ,atan(endy/endx)-3.14159), atan(endy/endx)) endlondeg = endlon*180/3.1415 ysinelon = endy / sin(endlon) endlat = -2 * ( atan(ysinelon) - pi4) endlatdeg = endlat*(180/3.1415) ;print(endlatdeg+" "+endlondeg) ;add lines for the a and b ellipse axes resp = True ; polyline mods desired resp@gsLineColor = "mediumblue" ; color of lines resp@gsLineThicknessF = 6.0 ; thickness of lines resp@gsLineDashPattern = 0 ; change to blue every time ; create array of dummy graphic variables. This is required, b/c each line ; must be associated with a unique dummy variable. ; find distance and points on the great circles that represent the major and minor ellipse axes a1gcircle := gc_latlon(endlatdeg(0),endlondeg(0),cenlatdeg,cenlondeg,100,4) a2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(1),endlondeg(1),100,4) b1gcircle := gc_latlon(endlatdeg(2),endlondeg(2),cenlatdeg,cenlondeg,100,4) b2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(3),endlondeg(3),100,4) ; dum := new(4,graphic) ; draw each line separately. Each line must contain two points. resp@gsLineLabelString= a1gcircle+" km" ; adds a line label string ; dum(0)=gsn_add_polyline(wks,fills,a1gcircle@gclon,a1gcircle@gclat,resp) ;a1 axis ; dum(1)=gsn_add_polyline(wks,fills,a2gcircle@gclon,a2gcircle@gclat,resp) ;a1 axis resp@gsLineLabelString= b1gcircle+" km" ; adds a line label string ; dum(2)=gsn_add_polyline(wks,fills,b1gcircle@gclon,b1gcircle@gclat,resp) ;b axis ; dum(3)=gsn_add_polyline(wks,fills,b2gcircle@gclon,b2gcircle@gclat,resp) ;b axis plyres@gsLineColor = "deepskyblue" ; change to blue ellipseplot = gsn_add_polyline(wks,fills,lons,lats,plyres) overlay(fills,ellipseplot) plyres@gsMarkerColor = "deepskyblue" plyres@gsMarkerIndex = 1 plyres@gsMarkerSizeF = 0.09 ellipseplot2 = gsn_add_polymarker(wks,fills,cenlondeg,cenlatdeg,plyres) overlay(fills,ellipseplot2) ;; make sure we write major axis first and that phi corresponds to angle of major axis phideg = phi * 180./3.1415 ; if a1 is greater than b1 if (a1gcircle .lt. b1gcircle) then ; if a1 is less than b1, adjust phi by 90 deg phideg = phi * 180./3.1415 - 90 end if print("Emetrics phi: "+phideg) newphideg = phideg if (phideg.lt.-45) then newphideg = phideg+180 end if ratio = a1gcircle/b1gcircle if (ratio.lt.1) then ratio = 1/ratio end if size=3.14159*a1gcircle*b1gcircle Ulabel = sprintf("%5.2f",dim_avg(ufcast(n,mbr))) Tlabel = sprintf("%5.2f",dim_avg(tfcast(n,mbr))) slabel = sprintf("%5.2e",size) philabel = sprintf("%5.2f",newphideg) rlabel = sprintf("%5.2f",ratio) lonlabel = sprintf("%5.2f",cenlondeg) ;;;;Create legend with Elliptical Diagnostics lgres = True lgres@lgLineColors = (/"white","white","white","white", "white", "white","white","white"/) lgres@lgLineThicknessF = 1. lgres@lgLabelFontHeightF = .38 ; set the legend label font thickness lgres@vpWidthF = 0.25 ; width of legend (NDC) lgres@vpHeightF = 0.25 ; height of legend (NDC) lgres@lgMonoDashIndex = True lgres@lgPerimColor = "gray60" ; draw the box perimeter in orange ;lgres@lgPerimThicknessF = 2.0 ; thicken the box perimeter lgres@lgPerimFill = 0 ; solid fill lgres@lgPerimFillColor = "White" labels = (/ "U65: "+Ulabel, "Polar Cap T: "+Tlabel, "Size: "+slabel, "Rotation: "+philabel+"~S~o~N~","Ratio: "+rlabel, "Center Lon: "+lonlabel+"~S~o~N~E ", "Center Lat: "+sprintf("%5.2f",cenlatdeg)+"~S~o~N~N ","~F22~ Ens Mean" /) ; Create the legend. lbid = gsn_create_legend(wks,8,labels,lgres) ; create legend ; Set up resources to attach legend to map. amres = True amres@amParallelPosF = 0.33 ; positive move legend to the right amres@amOrthogonalPosF = 0.33 ; positive move the legend up annoid1 = gsn_add_annotation(fills,lbid,amres) ; attach legend to plot enssize := new(nens+1,float,0) ensphi := new(nens+1,float,0) enslat := new(nens+1,float,0) enslon := new(nens+1,float,0) ensratio := new(nens+1,float,0) ensu := new(nens+1,float,0) enst := new(nens+1,float,0) ensseg := new(nens+1,float,0) ;if (isoline@n_points(i).eq.max_npts) then <-------add back in for segment looping ensseg(0) = segments enssize(0) = size ensphi(0) = newphideg enslat(0) = cenlatdeg enslon(0) = cenlondeg ensratio(0) = ratio aaa = a1gcircle bbb = b1gcircle ensu(0) = dim_avg(ufcast(n,mbr)) enst(0) = dim_avg(tfcast(n,mbr)) ;end if <-------add back in for segment looping ; end do; end do for segments in member = 0 <-------------------add back in for segments loops end if ;when mbr=0--------------------------------------------------- ;--------------------------------------------------------------------------------------------- ;------------------------Plotting for ENS members--------------------------- ;--------------------------------------------------------------------------------------------- if (mbr.ne.0) then g := garray(n,mbr-1,:,:) cnres@cnLevels = (/2850., 2900.,2950.,3000.,3050.,3100./) cntr(mbr) = gsn_csm_contour(wks1,g,cnres2) ; create the plot print("plotted cnres") ; info about contourline isoline := get_isolines(cntr(mbr),30000.) ;30000 for 10hpa, 23000 for 30hpa, 20000 move for 50 cntr(mbr) = gsn_csm_contour(wks,g,cnres) ; create the plot overlay(fills,cntr(mbr)) plyres@cnLevelSelectionMode = "ManualLevels" plyres@gsLineColor = -1 ; change to transparent plyres@gsLineThicknessF = 8. count = 0 segments := 0 segments := isoline@segment_count if (segments .gt. 2) then ;if there are more than 2 its probably an error segments := 2 ;so only do 2 end if ;IF USED FOR LEVELS .GT. 50 THIS PROBABLY NEEDS TO BE CHANGED do i = 0,segments -1 error = 0 ;reset error to 0 every time plyres@gsLineColor = -1;"mediumblue" ; change to blue every time b := isoline@start_point(i) e := b + isoline@n_points(i) - 1 ylat := isoline(0,:) xlon := isoline(1,:) if (dimsizes(ylat) .le. 50) then ;it its an error set error = 1 error = 1 end if if (i .eq. 0 .and. error .eq. 0) then ;there's only 1 vortex or its the first one plyplot1 = gsn_add_polyline(wks,fills,xlon,ylat,plyres) end if if (i .eq. 1 .and. error .eq. 0) then ;it really is a split vortex plyplot2 = gsn_add_polyline(wks,fills,xlon,ylat,plyres) end if count = count + isoline@n_points(i) print(isoline@level + " has " + count + " total points in " + isoline@segment_count + " segments" ) ; printVarSummary(xlon) ; print("Avg lon: "+dim_avg(xlon)) ; printVarSummary(ylat) ; print("Avg lat: "+dim_avg(ylat)) ; print(x) ; print(y) lats := ylat*(3.1415/180) lons := xlon*(3.1415/180) ;convert lat and lon to new coords with NP at origin and 0 lon as y ; Following Waugh (1997) ;;;;;;;Write contour lat/lon to file;;;;;;;;;;;;;;; x := (cos(lons) * cos(lats))/(1 + sin(lats)) asciiwrite("x.txt",x) y := (sin(lons) * cos(lats))/(1 + sin(lats)) asciiwrite("y.txt",y) ;;;;;;;;Run Python code to find Ellipse;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;uses x.txt and y.txt to find ellipse then makes xx.txt and yy.txt;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; cmd = "python fitEllipse2.py" system(cmd) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print("Ran Python") xx := asciiread("xx.txt",-1,"float") yy := asciiread("yy.txt",-1,"float") ;convert polar sterographic cartesian coords to lat/lon xxlon := atan(yy/xx) xxlon := where(xx.lt.0, where(yy.gt.0, atan(yy/xx) +3.1415 ,atan(yy/xx)-3.14159), atan(yy/xx)) sinxxlon = sin(xxlon) sinxxlon = where(sinxxlon.eq.0, sinxxlon@_FillValue,sinxxlon) ; don't divide by 0 yysinxxlon := yy / sinxxlon pi4 = 3.1415 / 4 yylat := -2 * ( atan(yysinxxlon) - pi4) lats := yylat*(180/3.1415) lons := xxlon*(180/3.1415) do l=1,dimsizes(lats)-1 if(abs(lats(l)-lats(l-1)).gt.1.5 .and. .not.ismissing(lats(l)) .and. .not.ismissing(lats(l-1)))then lats(l) = lats(l-1) end if end do ;read in major and minor axis lengths in cart. coords, center of ellipse, and rotation of major axis in radians vals := asciiread("axescenterphi.txt", -1,"float") print(vals) aax := vals(0) bax := vals(1) centerx := vals(2) centery := vals(3) phi := vals(4) ;convert cartisian to lat/lon coords ; if x = ( cos(lon)cos(lat) ) / (1 + sin(lat)) and y = ( sin(lon)cos(lat) ) / (1 + sin(lat)) ; then lon = atan(y/x) and lat = 2[atan( y/sin(lon) ) - pi/4] ; via tangent half angle formula (see http://en.wikipedia.org/wiki/Tangent_half-angle_formula ) cenlon := atan(centery / centerx) cenlon := where(centerx.lt.0, where(centery.gt.0, atan(centery/centerx) +3.1415 ,atan(centery/centerx)-3.14159), atan(centery/centerx)) cenlondeg := cenlon*180/3.1415 ysinlon := centery / sin(cenlon) cenlat := -2 * ( atan(ysinlon) - pi4) cenlatdeg := cenlat*(180/3.1415) print("center of ellipse = "+cenlatdeg+" N, "+cenlondeg+" E.") ;Calc enpoints of the axes of the vortex xa = (aax) * cos(phi) ya = (aax) * sin(phi) xb = (bax) * sin(phi) yb = (bax) * cos(phi) endx = (/ centerx+xa , centerx-xa , centerx+xb , centerx-xb/) endy = (/ centery+ya , centery-ya , centery-yb , centery+yb/) endlon = atan(endy / endx) endlon = where(endx.lt.0, where(endy.gt.0, atan(endy/endx) +3.1415 ,atan(endy/endx)-3.14159), atan(endy/endx)) endlondeg = endlon*180/3.1415 ysinelon = endy / sin(endlon) endlat = -2 * ( atan(ysinelon) - pi4) endlatdeg = endlat*(180/3.1415) ;print(endlatdeg+" "+endlondeg) ;add lines for the a and b ellipse axes resp = True ; polyline mods desired resp@gsLineColor = "firebrick" ; color of lines resp@gsLineThicknessF = 2.0 ; thickness of lines ; create array of dummy graphic variables. This is required, b/c each line ; must be associated with a unique dummy variable. ; find distance and points on the great circles that represent the major and minor ellipse axes a1gcircle := gc_latlon(endlatdeg(0),endlondeg(0),cenlatdeg,cenlondeg,100,4) a2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(1),endlondeg(1),100,4) b1gcircle := gc_latlon(endlatdeg(2),endlondeg(2),cenlatdeg,cenlondeg,100,4) b2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(3),endlondeg(3),100,4) if (i .eq. 0 .and. error .eq. 0) then plyres@gsLineThicknessF = 2.0 ; thickness of lines plyres@gsLineColor = modelcolors(j) ; change to blue ellipseplot := gsn_add_polyline(wks,cntr(mbr),lons,lats,plyres) overlay(fills,ellipseplot) end if if (i .eq. 1 .and. error .eq. 0) then plyres@gsLineThicknessF = 2.0 ; thickness of lines plyres@gsLineColor = modelcolors(j) ; change to blue ellipseplot2 := gsn_add_polyline(wks,cntr(mbr),lons,lats,plyres) overlay(fills,ellipseplot2) end if ;; make sure we write major axis first and that phi corresponds to angle of major axis phideg = phi * 180./3.1415 ; if a1 is greater than b1 if (a1gcircle .lt. b1gcircle) then ; if a1 is less than b1, adjust phi by 90 deg phideg = phi * 180./3.1415 - 90 end if print("Emetrics phi: "+phideg) newphideg = phideg if (phideg.lt.-45) then newphideg = phideg+180 end if ratio = a1gcircle/b1gcircle if (ratio.lt.1) then ratio = 1/ratio end if size=3.14159*a1gcircle*b1gcircle max_npts = max(isoline@n_points) print("n: "+n+", mbr: "+mbr) if (isoline@n_points(i).eq.max_npts) then ensseg(mbr) = segments enssize(mbr) = size ensphi(mbr) = newphideg enslat(mbr) = cenlatdeg enslon(mbr) = cenlondeg aaa = a1gcircle bbb = b1gcircle ensratio(mbr) = ratio ensu(mbr) = ufcast(n,mbr-1) enst(mbr) = tfcast(n,mbr-1) end if if (mbr.eq.nens) then do w=0,dimsizes(enslat)-1 print("In loop at L743, w = "+w) plyres@gsMarkerColor = modelcolors(j) ; if (ensseg(mbr).gt.1) then ; plyres@gsMarkerColor = "lavenderblush3" ; if split change line color <-----add for segment ; end if plyres@gsMarkerIndex = 4 plyres@gsMarkerSizeF = 0.015 ellipseplot2 := gsn_add_polymarker(wks,cntr(mbr),enslon(w),enslat(w),plyres) if(.not.ismissing(enslon(w)) .and. .not.ismissing(enslat(w)) ) then overlay(fills,ellipseplot2) end if end do ;;;;Create legend with Elliptical Diagnostics lgres = True lgres@lgLineColors = (/"white","white","white","white", "white", "white","white","white"/) lgres@lgLineThicknessF = 1. lgres@lgLabelFontHeightF = .38 ; set the legend label font thickness lgres@vpWidthF = 0.25 ; width of legend (NDC) lgres@vpHeightF = 0.28 ; height of legend (NDC) lgres@lgMonoDashIndex = True lgres@lgPerimColor = "slategray" ; draw the box perimeter in orange ;lgres@lgPerimThicknessF = 3.0 ; thicken the box perimeter lgres@lgPerimFill = 0 ; solid fill lgres@lgPerimFillColor = "White" labels = (/"U65: "+sprintf("%5.2f",dim_avg(ensu(1:)) ),"Polar Cap T: "+sprintf("%5.2f",dim_avg(enst(1:)) ),"Size: "+sprintf("%5.2e",dim_avg(enssize(1:))),"Rotation: "+sprintf("%5.2f",dim_avg(ensphi(1:)))+"~S~o~N~", "Ratio: "+sprintf("%5.2f",dim_avg(ensratio(1:))), "Center Lon: "+sprintf("%5.2f",dim_avg(enslon(1:)))+"~S~o~N~E ", "Center Lat: "+sprintf("%5.2f",dim_avg(enslat(1:)))+"~S~o~N~N ","~F22~ Ens. Mean" /) ; Create the legend. ; lbid = gsn_create_legend(wks,8,labels,lgres) ; create legend ; Set up resources to attach legend to map. ; amres = True ; amres@amParallelPosF = 0.4 ; positive move legend to the right ; amres@amOrthogonalPosF = 0.5 ; positive move the legend up ; annoid1 = gsn_add_annotation(fills,lbid,amres) ; attach legend to plot end if ;mbr=10 end do ;ending segment do end if ;------------------------------------------------------------- nmbr = n*nens + mbr ;;;;write the Elliptical diagnostics to a file;;;;;;; temp := hhstr(n)+" "+ddstr(n)+" "+mmstr(n)+" "+yystr(n)+" "+fff+" "+modellist(j)+" "+mbr+" "+enslon(mbr)+" "+enslat(mbr)+" "+aaa+" "+bbb+" "+ensphi(mbr)+" "+ensu(mbr)+" "+enst(mbr)+" "+ensseg(mbr) print(temp) metrics(nmbr) = temp end do; loop through members draw(fills) pres = True maximize_output(wks,pres) delete(wks1) end do ;;;;;;;;;;;;;;;;;;;;;;;;;end the time do loop;;;;;;;;; asciiwrite(modelname+"metrics"+year1+mm1+dd1+".txt", metrics) end do ;;;j end