sct expected at most 15 arguments, got 16 #25
Comments
|
Hi Gianpaolo, I believe we changed the function signature for titanlib.sct() since that example has been posted. Thomas Nipen (@tnipen) can confirm this. https://github.com/metno/gridpp https://github.com/metno/titanlib For instance, in the file Hope this helps, |
|
[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 GitHub<#25 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AA6FPMDEXBBXOJWNCEKUP73XLATAJANCNFSM6AAAAAAZDW3ZRQ>.
You are receiving this because you authored the thread.Message ID: ***@***.***>
|
|
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=180*4
nlon=360*4
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 GitHub<#25 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AA6FPMDEXBBXOJWNCEKUP73XLATAJANCNFSM6AAAAAAZDW3ZRQ>.
You are receiving this because you authored the thread.Message ID: ***@***.***>
|
|
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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 &)
The text was updated successfully, but these errors were encountered: