search

metno / titanlib

Automatic data quality control software library
GNU Lesser General Public License v3.0
32 stars 9 forks source link

sct expected at most 15 arguments, got 16 #25

Open gpbalsamo opened 1 year ago

gpbalsamo commented 1 year ago

Hi, In the example posted here https://tnipen.github.io/2020/06/15/titanlib-gridpp.html the call to sct generates an error and I report here a log:

TypeError Traceback (most recent call last) Cell In[142], line 14 11 t2neg = np.full(len(obs), 4) 12 eps2 = np.full(len(obs), 0.5) ---> 14 flags, sct, rep = titanlib.sct(obs_lats, obs_lons1, obs_elevs, obs, num_min, num_max, inner_radius, outer_radius, num_iterations, num_min_prof, dzmin, dhmin , dz, t2pos, t2neg, eps2)

File ~/.local/lib/python3.10/site-packages/titanlib.py:883, in sct(args) 882 def sct(args): --> 883 return _titanlib.sct(*args)

TypeError: sct expected at most 15 arguments, got 16 Additional information: Wrong number or type of arguments for overloaded function 'sct'. Possible C/C++ prototypes are: titanlib::sct(titanlib::Points const &,titanlib::vec const &,int,int,float,float,int,int,float,float,float,titanlib::vec const &,titanlib::vec const &,titanlib::vec const &,titanlib::vec &,titanlib::vec &,titanlib::ivec const &) titanlib::sct(titanlib::Points const &,titanlib::vec const &,int,int,float,float,int,int,float,float,float,titanlib::vec const &,titanlib::vec const &,titanlib::vec const &,titanlib::vec &,titanlib::vec &)

cristianlussana commented 1 year ago

Hi Gianpaolo,

I believe we changed the function signature for titanlib.sct() since that example has been posted. Thomas Nipen (@tnipen) can confirm this. I suggest you look more into the examples in the two repositories:

https://github.com/metno/gridpp (wiki https://github.com/metno/gridpp/wiki) tests: https://github.com/metno/gridpp/tree/master/tests

https://github.com/metno/titanlib (wiki https://github.com/metno/titanlib/wiki) queste due pagine sono sicuro che sono aggiornate. tests: https://github.com/metno/titanlib/tree/master/tests

For instance, in the file https://github.com/metno/titanlib/blob/master/tests/sct_test.py you may find an example of how to use titanlib.sct(), which is updated with the most recent version of the function. See here: https://github.com/metno/titanlib/blob/master/include/titanlib.h (from line 67 to line 102)

Hope this helps, Cristian

gpbalsamo commented 1 year ago

[like] Gianpaolo Balsamo reacted to your message:

From: Cristian Lussana @.> Sent: Tuesday, June 13, 2023 8:53:24 AM To: metno/titanlib @.> Cc: Gianpaolo Balsamo @.>; Author @.> Subject: Re: [metno/titanlib] sct expected at most 15 arguments, got 16 (Issue #25)

Hi Gianpaolo,

I believe we changed the function signature for titanlib.sct() since that example has been posted. Thomas Nipen @.***https://github.com/tnipen) can confirm this. I suggest you look more into the examples in the two repositories:

https://github.com/metno/gridpp (wiki https://github.com/metno/gridpp/wiki) tests: https://github.com/metno/gridpp/tree/master/tests

https://github.com/metno/titanlib (wiki https://github.com/metno/titanlib/wiki) queste due pagine sono sicuro che sono aggiornate. tests: https://github.com/metno/titanlib/tree/master/tests

For instance, in the file https://github.com/metno/titanlib/blob/master/tests/sct_test.py you may find an example of how to use titanlib.sct(), which is updated with the most recent version of the function. See here: https://github.com/metno/titanlib/blob/master/include/titanlib.h (from line 67 to line 102)

Hope this helps, Cristian

— Reply to this email directly, view it on GitHubhttps://github.com/metno/titanlib/issues/25#issuecomment-1588846587, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AA6FPMDEXBBXOJWNCEKUP73XLATAJANCNFSM6AAAAAAZDW3ZRQ. You are receiving this because you authored the thread.Message ID: @.***>

gpbalsamo commented 1 year ago

Hi Cristian and Thomas,

Many thanks for your prompt reply.

First of all, congrats on the gridpp titanlib developments, which are really useful and easy to use. I add Tiago who is leading development of several libraries here at ECMWF and we both got interested.

I was trying to update the example on this page as it was closer to what I was trying to go, a T2m OI analysis on the IFS grid fields and the new Earthkit library Tiago has produced.

I have succeeded reproducing the OI analysis for T2m.

Attached a python example (as jupyter notebook, but I also reproduce it below in ascii).

What we would like to know is how the optimal interpolation are done and if possible to access scripts example codes that allow to reproduce for instance what is done in this study https://www.diva-portal.org/smash/get/diva2:1613464/FULLTEXT01.pdf

We are also generally interested in the QC aspects of observations. Personally I am also very interested in the ensemble version documented in Cristian paper.

All suggestions welcome. All the best, Gianpaolo

[optimum-interpolation-of-2m-air-temperature]] == Optimum Interpolation of 2m Air Temperature

+In[1]:+ [source, ipython3]

Test a T2m OI

import netCDF4

import matplotlib.pyplot as mpl

import gridpp import titanlib import numpy as np import metview as mv import earthkit.data

+In[2]:+ [source, ipython3]

!pip3 install titanlib gridpp #if the above do not work install gridpp and titanlib

+In[3]:+ [source, ipython3]

date_ext="2023-06-09" time_ext="12" step_ext="0" obs_name="airTemperatureAt2M" fg_name=["2t","z"] nlat=1804 nlon=3604 grid_ext=[180./nlat,360./nlon] print (obs_name,"analysis increment on latlon grid :",grid_ext)

obs_name="windSpeedAt10M"

fg_name="10u"

obs_name="dewpointTemperatureAt2M"

fg_name="2d"

obs_name="totalSnowDepth"

fg_name="sd"

+Out[3]:+

airTemperatureAt2M analysis increment on latlon grid : [0.25, 0.25]

+In[4]:+ [source, ipython3]

fc = earthkit.data.from_source ("mars", param=fg_name,levtype="sfc",grid=grid_ext, date=date_ext, time=time_ext, step=step_ext)

+In[5]:+ [source, ipython3]

obs = earthkit.data.from_source("mars",type="ob",repres="bu",time=time_ext, step="0",date=date_ext)

+In[6]:+ [source, ipython3]

an_inc = mv.retrieve({'expver':"1",'stream':"oper",'levtype':"sfc",'date':date_ext,'time':time_ext,'step':step_ext,'type': "fc",'param':["2t"], 'grid' : grid_ext})

+In[7]:+ [source, ipython3]

from ecmwfapi import ECMWFService

server = ECMWFService("mars")

+In[8]:+ [source, ipython3]

server.execute({"class": "od","date": date_ext,"time" : "6","type": "ob", "repres":"bu"},"synop_obs.bufr")

+In[9]:+ [source, ipython3]

filename = "synop_obs.bufr"

synop = mv.read(filename)

+In[10]:+ [source, ipython3]

synop_t2m = obs = mv.obsfilter(output="ncols",parameter=["heightOfStation","airTemperatureAt2M"],missing_data="exclude",data=synop)

+In[11]:+ [source, ipython3]

synop_t2m.save("synop_obs.bufr")

+In[12]:+ [source, ipython3]

to examine the bufr

obs.save("synop.bufr")

codes_ui -b synop.bufr

+In[13]:+ [source, ipython3]

synop_t2m = obs.to_pandas(columns=["latitude", "longitude", "heightOfStation",obs_name], filters={ obs_name:lambda x: x==x, "heightOfStation":lambda x: x==x })

+In[14]:+ [source, ipython3]

synop_t2m

+Out[14]:+

[cols=",,,,",options="header",] |========================================================== | |latitude |longitude |heightOfStation |airTemperatureAt2M |0 |-4.26 |15.24 |316 |302.8 |1 |-3.44 |114.75 |20 |299.0 |2 |18.41 |103.52 |298 |298.4 |3 |-15.57 |-175.62 |61 |298.0 |4 |-15.95 |-173.77 |3 |297.7 |... |... |... |... |... |4445 |49.63 |6.23 |369 |298.3 |4446 |32.37 |36.25 |683 |303.6 |4447 |31.98 |35.98 |779 |302.2 |4448 |31.72 |35.98 |722 |303.8 |4449 |-62.23 |-58.79 |11 |267.3 |==========================================================

4450 rows × 4 columns

+In[15]:+ [source, ipython3]

obs_lats = synop_t2m["latitude"] obs_lons = synop_t2m["longitude"] obs_elevs = synop_t2m["heightOfStation"]*1.0 obs = synop_t2m[obs_name]

+In[16]:+ [source, ipython3]

obs_lons1 =np.mod(obs_lons,360)

+In[17]:+ [source, ipython3]

inner_radius = 50000 outer_radius = 100000 num_min = 5 num_max = 100 num_iterations = 1 num_min_prof = 20 dzmin = 100 dhmin = 10000 dz = 200 t2pos = np.full(len(obs), 4) t2neg = np.full(len(obs), 4) eps2 = np.full(len(obs), 0.5)

flags = 0

+In[18]:+ [source, ipython3]

obs_points=titanlib.Points(obs_lats, obs_lons1, obs_elevs)

+In[19]:+ [source, ipython3]

flags, sct, rep = titanlib.sct(obs_points, obs, num_min, num_max, inner_radius, outer_radius, num_iterations, num_min_prof, dzmin, dhmin , dz, t2pos, t2neg, eps2)

+In[20]:+ [source, ipython3]

index_valid_obs = np.where(flags == 0)[0] index_invalid_obs = np.where(flags != 0)[0]

+In[21]:+ [source, ipython3]

print ("Number of invalid observations from sct test:",np.sum(flags))

+Out[21]:+

Number of invalid observations from sct test: 36

+In[22]:+ [source, ipython3]

valid_obs_lats=np.copy(obs_lats[index_valid_obs]) valid_obs_lons=np.copy(obs_lons1[index_valid_obs]) valid_obs_elevs=np.copy(obs_elevs[index_valid_obs]) valid_obs=np.copy(obs[index_valid_obs])

+In[23]:+ [source, ipython3]

print (obs_elevs)

print (obs)

+In[24]:+ [source, ipython3]

Load First Guess - Metview version

index_control = 0

blats = fc.select(shortName="z").latitudes().reshape((721, 1440))

blons = fc.select(shortName="z").longitudes().reshape((721, 1440))

belevs = fc.select(shortName="z").values().reshape((721, 1440))/ 9.81

bgrid = gridpp.Grid(blats, blons, belevs)

background = fc.select(shortName=fg_name).values().reshape((721, 1440))

+In[25]:+ [source, ipython3]

Load First Guess - Earthkit version

index_control = 0 latlon = fc[0].to_points(flatten=False) blats = latlon["lat"].reshape((nlat+1, nlon)) blons = latlon["lon"].reshape((nlat+1, nlon)) belevs= fc[1].to_numpy(flatten=False).reshape((nlat+1, nlon)) bgrid = gridpp.Grid(blats, blons, belevs) background =fc[0].to_numpy(flatten=False).reshape((nlat+1, nlon))

+In[26]:+ [source, ipython3]

blons1 =np.mod(blons,360)

+In[27]:+ [source, ipython3]

Load Analysis template

bgrid = gridpp.Grid(blats, blons1, belevs)

+In[28]:+ [source, ipython3]

Downscaling FG on Analysis template

gradient = -0.0065 background = gridpp.simple_gradient(bgrid, bgrid, background, gradient) points = gridpp.Points(valid_obs_lats, valid_obs_lons, valid_obs)

points = gridpp.Points(obs_lats, obs_lons1, obs_elevs)

pbackground = gridpp.bilinear(bgrid, points, background)

variance_ratios = np.full(points.size(), 0.5)

h = 300000 v = 800 structure = gridpp.BarnesStructure(h, v)

max_points = 50

+In[29]:+ [source, ipython3]

print (pbackground)

print (obs)

print (pbackground-obs)

+In[30]:+ [source, ipython3]

Perform the OI Analysis

analysis = gridpp.optimal_interpolation(bgrid, background, points, valid_obs, variance_ratios, pbackground, structure, max_points)

+In[31]:+ [source, ipython3]

background

+In[32]:+ [source, ipython3]

Calculate the increments

diff = analysis - background

+In[33]:+ [source, ipython3]

Plotting the increments with Matplotlib

diff = analysis - background

mpl.pcolormesh(blons, blats, diff, cmap="RdBu_r", vmin=-15, vmax=15, shading="auto")

mpl.xlim(0,360)

mpl.ylim(-90, 90)

mpl.gca().set_aspect(2)

cb = mpl.colorbar()

cb.set_label(r"Increment ($\degree C$)")

+In[34]:+ [source, ipython3]

diff.max()

+Out[34]:+ ----3.740265----

+In[35]:+ [source, ipython3]

diff.min()

+Out[35]:+ -----3.257904----

+In[36]:+ [source, ipython3]

lons=blons.flatten() lats=blats.flatten() vals=diff.flatten()

+In[37]:+ [source, ipython3]

an_inc=an_inc.set_values(vals)

+In[38]:+ [source, ipython3]

mv.write("aninc_2t.grib",an_inc)

+Out[38]:+ ----0.0----

[[plotting]] == Plotting

+In[39]:+ [source, ipython3]

set a single contour set-up for differences comparison

cont_diff = mv.mcont(contour_shade="on", contour_shade_colour_method="palette", contour_shade_method="area_fill", contour_level_count=20, contour_shade_palette_name="eccharts_black_red_21", contour_shade_technique = "grid_shading", contour_max_level=10.0, contour_min_level=-10.0, legend="on", contour="off", contour_label="off",)

+In[40]:+ [source, ipython3]

fieldc = mv.mcont( contour_level_selection_type="INTERVAL", contour_shade="ON", contour_interval=5, contour_shade_max_level=300., contour_shade_min_level=200, contour_shade_method="AREA_FILL", contour_shade_colour_method="CALCULATE", contour_shade_colour_direction="CLOCKWISE", contour_shade_max_level_colour="RED", contour_shade_min_level_colour="BLUE", contour_line_colour="GREY", contour_line_thickness=2, contour_highlight="OFF", contour_label="OFF", legend="ON" )

+In[41]:+ [source, ipython3]

define coastlines

coastlines = mv.mcoast( map_coastline_resolution="high", map_coastline_sea_shade="on", map_coastline_sea_shade_colour="RGB(0.4845,0.6572,0.9351)", )

define coastlines

coastlines = mv.mcoast( map_coastline_resolution="high", map_coastline_sea_shade="on", map_coastline_sea_shade_colour="RGB(0.4845,0.6572,0.9351)", map_grid_frame_colour="grey", map_grid_frame_line_style="dash", map_grid_frame="off", map_coastline="on", map_label="off",

map_grid="off",

)

define views

view_coast = mv.geoview(coastlines=coastlines, page_frame="off",subpage_frame="off") view_nh = mv.make_geoview(area="north_pole") view_ce = mv.make_geoview(area="central_europe") view_na = mv.make_geoview(area="north_america") view_sa = mv.make_geoview(area="southern_asia") view_eu = mv.make_geoview(area="europe") view_eu_coast = mv.geoview(area_mode="name", area_name="europe", coastlines=coastlines, page_frame="off",subpage_frame="off") view_alps = mv.geoview(map_area_definition="corners", area=[40, 0, 55, 20]) view_himalaya=mv.geoview(map_area_definition="corners", area=[22, 68, 46, 110]) view_alps_details = mv.geoview(map_area_definition="corners", area=[43, 4, 49, 17])

+In[42]:+ [source, ipython3]

mv.setoutput("jupyter", output_width=1200, output_font_scale=3, plot_widget=False)

+In[43]:+ [source, ipython3]

mv.plot(view_nh,an_inc,cont_diff)

From: Cristian Lussana @.> Date: Tuesday, 13 June 2023 at 09:53 To: metno/titanlib @.> Cc: Gianpaolo Balsamo @.>, Author @.> Subject: Re: [metno/titanlib] sct expected at most 15 arguments, got 16 (Issue #25)

Hi Gianpaolo,

I believe we changed the function signature for titanlib.sct() since that example has been posted. Thomas Nipen @.***https://github.com/tnipen) can confirm this. I suggest you look more into the examples in the two repositories:

https://github.com/metno/gridpp (wiki https://github.com/metno/gridpp/wiki) tests: https://github.com/metno/gridpp/tree/master/tests

https://github.com/metno/titanlib (wiki https://github.com/metno/titanlib/wiki) queste due pagine sono sicuro che sono aggiornate. tests: https://github.com/metno/titanlib/tree/master/tests

For instance, in the file https://github.com/metno/titanlib/blob/master/tests/sct_test.py you may find an example of how to use titanlib.sct(), which is updated with the most recent version of the function. See here: https://github.com/metno/titanlib/blob/master/include/titanlib.h (from line 67 to line 102)

Hope this helps, Cristian

— Reply to this email directly, view it on GitHubhttps://github.com/metno/titanlib/issues/25#issuecomment-1588846587, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AA6FPMDEXBBXOJWNCEKUP73XLATAJANCNFSM6AAAAAAZDW3ZRQ. You are receiving this because you authored the thread.Message ID: @.***>

tnipen commented 1 year ago

Dear Gianpaolo.

Thanks for letting us know about the incorrect code in the blog post. I have updated it to reflect the changes that were made to titanlib.

Thomas