[ncl-talk] R: Wrong behaviour of lat/lon coordinates in an interpolated field

Mon Jul 5 16:22:19 MDT 2021

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210705/02c78aa6/attachment.html>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Vento_Velocit?_ISPRA.txt
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210705/02c78aa6/attachment.txt>