[ncl-talk] R: Wrong behaviour of lat/lon coordinates in an interpolated field
Dennis Shea
shea at ucar.edu
Mon Jul 5 20:50:09 MDT 2021
Currently, I don't have the time to look at the coding in detail.
I am going to suggest that you work with a scalar field like pressure to
debug this code.
Then add the wind components.
I did look at the text files. In many of the files the last value is not
present. For example:
vento_direzione_201911.TXT
The 1st 3 lines
12/11/19;00:00;39;42;53;60;35;47;33;272;359
12/11/19;00:05;37;40;51;60;32;48;31;111; <=== last value is
missing [MANY warning messages]
12/11/19;00:10;35;43;48;64;36;48;30;356;23
I added some code to eliminate Warning messages
The attached script illustrates some of the issues.
%> ncl debug_3.ncl | less
On Mon, Jul 5, 2021 at 4:22 PM Federico Buscemi - federico.buscemi--- via
ncl-talk <ncl-talk at mailman.ucar.edu> wrote:
> Hello,
> I hope you can see it now. Thank you.
> ------------------------------
> *Da:* Federico Buscemi - federico.buscemi at studio.unibo.it <
> federico.buscemi at studio.unibo.it>
> *Inviato:* lunedì 5 luglio 2021 23:01
> *A:* ncl-talk at mailman.ucar.edu <ncl-talk at mailman.ucar.edu>
> *Oggetto:* Wrong behaviour of lat/lon coordinates in an interpolated field
>
> Dear all,
>
> I am having some trouble with my script, copied and pasted below. I
> interpolated sea level pressure values (z), read from two of the attached
> text files, recorded in six locations around Venezia, Italy. I wanted to
> plot the interpolated field in the Venezia Lagoon area. Even if values are
> well interpolated by means of "natgrid" function, it seems that they are
> not associated with latitude and longitude coordinates, since contours are
> always drawn in the same way, although I vary lat/lon coordinates of the
> map. Since I don't get errors in my output, I don't know how to solve this
> issue. I tried in several ways, but still can't figure it out. Any
> suggestion would be appreciated. Thank you.
>
> Federico Buscemi
>
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> begin
>
> fname1 = "filetxt/vento_direzione_201911.TXT"
> fname2 = "filetxt/vento_velocita_201911.TXT"
> fname3 = "filetxt/Vento_Direzione_ISPRA.txt"
> fname4 = "filetxt/Vento_Velocità_ISPRA.txt"
> fname5 = "filetxt/pressione_201911.TXT"
> fname6 = "filetxt/Pressioni_ISPRA.txt"
> lines1 = asciiread(fname1,-1,"string")
> lines2 = asciiread(fname2,-1,"string")
> lines3 = asciiread(fname3,-1,"string")
> lines4 = asciiread(fname4,-1,"string")
> lines5 = asciiread(fname5,-1,"string")
> lines6 = asciiread(fname6,-1,"string")
>
> delim = ";"
> ndata = str_fields_count(lines1, delim)
>
> date = str_get_field(lines1, 1, delim)
> time = str_get_field(lines1, 2, delim)
>
> dpiatt = stringtoint(str_get_field(lines1, 3, delim))
> dmalpor = stringtoint(str_get_field(lines1, 6, delim))
> dchipor = stringtoint(str_get_field(lines1, 7, delim))
> dsang = stringtoint(str_get_field(lines1, 9, delim))
> dmeab = stringtoint(str_get_field(lines1, 11, delim))
> dpame = stringtoint(str_get_field(lines3, 3, delim))
> dfopo = stringtoint(str_get_field(lines3, 5, delim))
>
> vpiatt = stringtofloat(str_get_field(lines2, 3, delim))
> vmalpor = stringtofloat(str_get_field(lines2, 6, delim))
> vchipor = stringtofloat(str_get_field(lines2, 7, delim))
> vsang = stringtofloat(str_get_field(lines2, 9, delim))
> vmeab = stringtofloat(str_get_field(lines2, 11, delim))
> vpame = stringtofloat(str_get_field(lines4, 3, delim))
> vfopo = stringtofloat(str_get_field(lines4, 5, delim))
>
> ppiatt = stringtofloat(str_get_field(lines5, 3, delim))
> ppaca = stringtofloat(str_get_field(lines5, 4, delim))
> pmeab = stringtofloat(str_get_field(lines5, 5, delim))
> ppame = stringtofloat(str_get_field(lines6, 3, delim))
> pfopo = stringtofloat(str_get_field(lines6, 4, delim))
> plime = stringtofloat(str_get_field(lines6, 5, delim))
>
> wpiatt = wind_component(vpiatt,dpiatt,0)
> wmalpor = wind_component(vmalpor,dmalpor,0)
> wchipor = wind_component(vchipor,dchipor,0)
> wsang = wind_component(vsang,dsang,0)
> wmeab = wind_component(vmeab,dmeab,0)
> wpame = wind_component(vpame,dpame,0)
> wfopo = wind_component(vfopo,dfopo,0)
>
> latpiatt = 45.3142417
> lonpiatt = 12.5083139
> latmalpor = 45.3398
> lonmalpor = 12.2919667
> latchipor = 45.2325389
> lonchipor = 12.2805972
> latsang = 45.4284083
> lonsang = 12.3462639
> latmeab = 45.235
> lonmeab = 12.7716667
> latpame = 45.4027694
> lonpame = 11.8585417
> latfopo = 44.8438333
> lonfopo = 12.4624028
> latpaca = 45.4364324
> lonpaca = 12.3335352
> latlime = 45.43015
> lonlime = 12.38295
>
> x = (/lonpame, lonfopo, lonpiatt, lonmeab, lonpaca, lonlime/)
> y = (/latpame, latfopo, latpiatt, latmeab, latpaca, latlime/)
> z = (/ppame(263), pfopo(263), ppiatt(263), pmeab(263), ppaca(263),
> plime(263)/)
>
> xo = fspan(11.8000, 12.8000, 4000)
> yo = fspan(44.8000, 45.4000, 2400)
>
> yo at units = "degrees_north"
> xo at units = "degrees_east"
>
> yo!0 = "lat"
> xo!0 = "lon"
>
> p_int = natgrid(x, y, z, xo, yo)
> p_int_real = transpose(p_int)
>
> p_int_real!0 = "lat"
> p_int_real!1 = "lon"
>
> p_int_real&lat = yo
> p_int_real&lon = xo
>
> p_int_real&lat at units = "degrees_north"
> p_int_real&lon at units = "degrees_east"
>
>
> wks = gsn_open_wks("pdf","station") ; send graphics to PNG file
>
> mpres = True
> mpres at gsnFrame = False
> mpres at mpLimitMode = "LatLon"
> mpres at mpMinLatF = 40.5
> mpres at mpMinLonF = 10
> mpres at mpMaxLatF = 48.8
> mpres at mpMaxLonF = 15.8
> mpres at mpDataBaseVersion = "HighRes"
> mpres at mpOutlineBoundarySets = "AllBoundaries"
> mpres at gsnMaximize = True
>
> map = gsn_csm_map(wks,mpres)
>
> cnres = True
> cnres at gsnAddCyclic = False
> cnres at gsnFrame = False
> cnres at cnFillOn = False
> cnres at cnLinesOn = True
> cnres at cnInfoLabelOn = False
> cnres at cnLevelSelectionMode = "ManualLevels"
> cnres at cnMinLevelValF = 980.0
> cnres at cnMaxLevelValF = 1010.0
> cnres at cnLevelSpacing = .5
> cnres at cnFillDrawOrder = "PreDraw"
> cnres at cnLineLabelsOn = True
> cnres at cnLineLabelPlacementMode = "Constant"
> ; cnres at sfXArray = xo
> ; cnres at sfYArray = yo
>
> cnres at tmXBBorderOn = False
> cnres at tmXBOn = False
> cnres at tmXTBorderOn = False
> cnres at tmXTOn = False
> cnres at tmYLBorderOn = False
> cnres at tmYLOn = False
> cnres at tmYRBorderOn = False
> cnres at tmYROn = False
>
>
> wmsetp("wbs", .045)
>
> wmbarbmap(wks, latpiatt, lonpiatt, wpiatt(0,263), wpiatt(1,263))
> wmbarbmap(wks, latsang, lonsang, wsang(0,263), wsang(1,263))
> wmbarbmap(wks, latmalpor, lonmalpor, wmalpor(0,263), wmalpor(1,263))
> wmbarbmap(wks, latchipor, lonchipor, wchipor(0,263), wchipor(1,263))
> wmbarbmap(wks, latmeab, lonmeab, wmeab(0,263), wmeab(1,263))
> wmbarbmap(wks, latpame, lonpame, wpame(0,263), wpame(1,263))
> wmbarbmap(wks, latfopo, lonfopo, wfopo(0,263), wfopo(1,263))
>
> cnp = gsn_csm_contour(wks,p_int_real,cnres)
>
> overlay(map,cnp)
>
> end
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210705/258e6365/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: debug_3.ncl
Type: application/octet-stream
Size: 2414 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210705/258e6365/attachment.obj>
More information about the ncl-talk
mailing list