lib/numru/gphys/interpolate.rb in gphys-1.1.1 vs lib/numru/gphys/interpolate.rb in gphys-1.2.2
- old
+ new
@@ -1,7 +1,8 @@
require "numru/gphys/gphys"
-require "numru/dcl" # math1/gt2dlib is used for pure-2D interpolation
+require "numru/dcl" # math1/gt2dlib is used for pure-2D interpolation.
+ # Also for dcl_fig_cut
require "numru/dcl_mouse" # for mouse_cut, mouse_cut_repeat
require "narray_miss"
if $0 == __FILE__
require "numru/gphys" # for test
@@ -11,13 +12,25 @@
class GPhys
@@interpo_previous_cutter = nil
@@interpo_previous_modifier = nil
- @@interpo_missval = -1e30
+ @@interpo_missval = 9.9692099683868690e+36 # NC_FILL_DOUBLE/FLOAT ~15*2^119
+ @@interpo_extrapolation = false
- # = Makes a subset interactively by specifying a (poly-)line on the DCL viewport
+ # Change the behavior of the interpolation methods to extrapolate
+ # outside the grid coverage.
+ #
+ # ARGUMENTS
+ # * extrapo : true or false --- the default behaviour is false
+ # (not to extrapolate), so use this method if you
+ # want to set it to true.
+ def self.extrapolation=(extrapo)
+ @@interpo_extrapolation = extrapo
+ end
+
+ # Makes a subset interactively by specifying a (poly-)line on the DCL viewport
#
# ARGUMENTS
# * dimx {String] : name of number (0,1,..) of the dimension
# corresponding to the X coordinate in the current window of DCL
# * dimy {String] : name of number (0,1,..) of the dimension
@@ -26,11 +39,11 @@
# (2 or greater -- if 2, a single line segment; if 3 or more, a
# poly-line)
#
# RETURN VALUE
# * a GPhys
- def mouse_cut(dimx, dimy, num=2)
+ def mouse_cut(dimx, dimy, num=2, line_type=1, line_index=1)
# < preparation >
dimx = dim_index(dimx)
dimy = dim_index(dimy)
@@ -44,16 +57,44 @@
puts "** The points specified include one(s) outside the U window. Do it again."
else
break
end
end
- line.draw
+ line.draw(line_type, line_index)
vx = line.vx
vy = line.vy
ux = line.ux
uy = line.uy
+ gpnew = dcl_fig_cut(dimx,dimy,ux,uy)
+ [gpnew, line]
+ end
+
+ # Interpolation on the DCL window (automatic iso-interval interpolation
+ # along a poly line that can be drawn in the current viewport of the
+ # DCL window). Used in mouse_cut.
+ #
+ # ARGUMENTS
+ # * dimx [Integer or String] : specifies the dimension corresponding
+ # to the UX coordinate.
+ # (Here, the UX coordinate is the X coordinate of the DCL's USER
+ # coordinate. For exapmle, longitude if map projection.)
+ # * dimy [Integer or String] : specifies the dimension corresponding
+ # to the UY coordinate.
+ # (Here, the UY coordinate is the Y coordinate of the DCL's USER
+ # coordinate. For exapmle, latitude if map projection.)
+ # * ux [Array] : x values in terms of the UX coordinate
+ # * uy [Array] : y values in terms of the UY coordinate
+ # Lengths of ux and uy must be the same and greter or equal to 2.
+ #
+ def dcl_fig_cut(dimx,dimy,ux,uy)
len = ux.length
+ raise("ux and uy must be arrays with the (same) length >= 2") if len<=1
+ raise("ux's len (#{len}) != uy's len (#{uy.length})") if uy.length != len
+ vx=Array.new; vy=Array.new
+ for i in 0...len
+ vx[i],vy[i] = NumRu::DCL.stftrf(ux[i],uy[i])
+ end
kx = Array.new
ky = Array.new
cut = [true]*rank
for i in 0...len
cut[dimx] = ux[i]
@@ -64,17 +105,22 @@
end
ndiv = Array.new
ndsum = [0]
for i in 0...len-1
ndiv[i] = Math.sqrt( (kx[i+1]-kx[i])**2 + (ky[i+1]-ky[i])**2).to_i
+ ndiv[i] += 1 if i==len-2
ndsum.push ndsum[-1] + ndiv[i] # 0, ndiv[0], ndiv[0]+ndiv[1], ...
end
ndtot = ndsum[-1]
vxdiv = NArray.float(ndtot)
vydiv = NArray.float(ndtot)
for i in 0...len-1
- a = NArray.float(ndiv[i]).indgen / ndiv[i]
+ if i!=len-2
+ a = NArray.float(ndiv[i]).indgen / ndiv[i]
+ else
+ a = NArray.float(ndiv[i]).indgen / (ndiv[i]-1)
+ end
vxdiv[ndsum[i]...ndsum[i+1]] = (1.0-a)*vx[i] + a*vx[i+1]
vydiv[ndsum[i]...ndsum[i+1]] = (1.0-a)*vy[i] + a*vy[i+1]
end
uxdiv = NArray.float(ndtot)
uydiv = NArray.float(ndtot)
@@ -113,14 +159,14 @@
# < do the job >
gpnew = interpolate(cutter)
gpnew = modifier[gpnew] if modifier
- [gpnew, line]
+ gpnew
end
- # = Interpolation onto grid points specified by the previous call of GPhys#mouse_cut
+ # Interpolation onto grid points specified by the previous call of GPhys#mouse_cut
def mouse_cut_repeat
if @@interpo_previous_cutter.nil?
raise("You must first use GPhys#mouse_cut. This method repeats it")
end
gpnew = interpolate(@@interpo_previous_cutter)
@@ -148,11 +194,11 @@
{"long_name"=>"distance along great circle","units"=>"degrees"},
"dist")
end
private :__sp_dist
- # = Interpolate to conform the grid to a target GPhys object
+ # Interpolate to conform the grid to a target GPhys object
#
# ARGUMENTS
# * to [GPhys] : the target gphys
#
# RETURN VALUE
@@ -161,12 +207,41 @@
def regrid(to)
coords = to.axnames.collect{|nm| to.coord(nm)}
interpolate(*coords)
end
- # = Wide-purpose multi-dimensional linear interpolation
+ # Reverse the main data (i.e., the dependent variable) and one of the
+ # coordinates (an independent variable) through interpolation.
#
+ # Returns a GPhys in which the main data is the specfied coordinate
+ # (argument: axname) sampled at specified locations (argument: pos)
+ # in terms of the main data of self. The main data of self is expected
+ # to be quai-monotonic with respect to the specfied coordinate.
+ #
+ # ARGUMENTS
+ # * axname [String] : one of the names of the axes (i.e. main
+ # coordinates. Auxiliary coordinates are not supported as the target.)
+ # * pos [NArray] : grid locations. For example, if the current data is
+ # potential temperature theta, pos consists of the theta levels to
+ # make sampling.
+ #
+ # RETURN VALUE
+ # * a GPhys
+ #
+ def coord_data_reverse(axname,pos)
+ gp = self.axis(axname).to_gphys
+ gp = self.shape_coerce_full(gp)[0] # conform the shape to that of self
+ gp = GPhys.new( gp.grid.copy, gp.data ) # copy grid to avoid side effect
+ # on the grid of self
+ gp.set_assoc_coords([self])
+ pos = NArray[*pos].to_type(NArray::FLOAT) if pos.is_a?(Array)
+ newcrd = VArray.new(pos,self.data,self.name) # succeeds the attributes
+ gp.interpolate(axname=>newcrd)
+ end
+
+ # Wide-purpose multi-dimensional linear interpolation
+ #
# This method supports interpolation regarding combinations of
# 1D and 2D coordinate variables. For instance, suppose self is
# 4D with coordinates named ["x", "y", "z", "t"] and associated
# coordinates "sigma"["z"] ("sigma" is 1D and its axis is "z"),
# "p"["x","y"], "q"["x","y"] ("p" and "q" are 2D having the
@@ -318,20 +393,29 @@
z = val
if z.is_a?(NArrayMiss)
missval = ( (a=get_att('_FillValue')) ? a[0] : nil ) ||
( (a=get_att('missing_value')) ? a[0] : nil ) ||
@@interpo_missval
- z = z.to_na(missval)
+ z = z.to_na(missval)
+ input_nomiss = false
else
- missval = nil
+ input_nomiss = true
+ if @@interpo_extrapolation
+ missval = nil
+ else
+ missval = @@interpo_missval
+ end
end
- na = c_interpo_do(newgrid.shape, idxmap, z, missval) # [C-extension]
+ na = c_interpo_do(newgrid.shape, idxmap, z, missval,
+ @@interpo_extrapolation) # [C-extension]
- if missval
+ if !input_nomiss || !@@interpo_extrapolation
mask = na.ne(missval)
- na = NArrayMiss.to_nam_no_dup(na,mask)
+ if !input_nomiss || mask.min == 0
+ na = NArrayMiss.to_nam_no_dup(na,mask)
+ end
end
va = VArray.new(na, data, name)
ret = GPhys.new(newgrid, va)
@@ -351,11 +435,12 @@
if mp.length==1
idxmap.push( [m[0], cd] ) # simple copying
elsif cd.is_a?(Numeric) && mp[1].is_a?(NArray)
xto = mp[1] # 1-D new coordinate var
xfrom = mp[2] # 1-D original coordinate var
- ids, f = c_interpo_find_loc_1D(xto,xfrom) # [C-extension]
+ ids, f = c_interpo_find_loc_1D(xto,xfrom,@@interpo_missval,
+ @@interpo_extrapolation) # [C-extension]
idxmap.push( [ m[0], cd, nil, ids, f ] ) # mapping from 1D
else
# partially 2D case
cdims = mp[2]
xto = mp[3] # 1-D new coordinate var
@@ -365,11 +450,12 @@
if cdims[i] == od
dimc = i
break
end
end
- ids, f = c_interpo_find_loc_1D_MD(xto,xfrom,dimc) # [C-extension]
+ ids, f = c_interpo_find_loc_1D_MD(xto,xfrom,dimc,@@interpo_missval,
+ @@interpo_extrapolation) # [C-extension]
dims_covd = mp[1] #dimensions covered by the coordinate variable(orig)
idxmap.push( [ od, cd, dims_covd, ids, f ] ) # mapping from 2D
#^^
#will be removed : see (***) below
@@ -460,14 +546,18 @@
crdmap1D.push( [ org_dims[ic], ic ] )
else
for j in 0...coords[ic].length
xto = coords[ic][j].val
xfrom = org_coords[ic][j].val
- xto = xto.to_na if !xto.is_a?(NArray)
+ xto = xto.to_na if !xto.is_a?(NArray)
# missing in the coordinate, if any, is ignored
- xfrom = xfrom.to_na if !xfrom.is_a?(NArray)
- # missing in the coordinate, if any, is ignored
+ if xfrom.is_a?(NArrayMiss)
+ xfrom = xfrom.to_type(NArray::FLOAT).to_na(@@interpo_missval)
+ # if xfrom is NArrayMiss, this fixed missing value is set.
+ # The conversion into double is just in case (not needed for
+ # the default @@interpo_missval).
+ end
if org_dims[ic][j].length == 1
crdmap1D.push( [ org_dims[ic][j][0], ic, xto,
xfrom] )
else
crdmap2D.push( [ org_dims[ic][j], ic, xto,
@@ -729,31 +819,33 @@
ga = gd + gx
ga.name = "a"
gd.set_assoc_coords([gp,gq,gr,ga])
- print "GPhys with associated coordinates:\n"
- p gd
+ #print "GPhys with associated coordinates:\n"
+ #p gd
DCL.swpset('iwidth',700)
DCL.swpset('iheight',700)
#DCL.sgscmn(4) # set colomap
DCL.gropn(1)
+ DCL.sgpset('isub', 96) # control character of subscription: '_' --> '`'
DCL.glpset("lmiss",true)
DCL.sldiv("y",2,2)
GGraph::set_fig "viewport"=>[0.15,0.85,0.15,0.85]
GGraph::tone gd
GGraph::color_bar
GGraph::tone gd[true,ny/2,true]
GGraph::color_bar
+ gdd = gd.copy
#< prepare coordinates to interpolate >
xi = NArray[1.0, 2.0, 3.0, 4.0, 5.0]
# yi = NArray[1.0, 4.0, 5.0]
- yi = NArray[-0.1, 2.5, 4.0, 5.5, 6.8] # test of extrapolation
+ yi = NArray[-0.1,2.5, 4.0, 5.5, 6.8] # test of extrapolation
vxi = VArray.new( xi, {"units"=>"m"}, "x") # "0.5m" to test unit conversion
vyi = VArray.new( yi, {"units"=>"m"}, "y") # "0.5m" to test unit conversion
# pi = NArray[10.0, 13.0, 15.0, 17.0, 20.0]
# qi = NArray[0.0, 3.0, 5.0, 7.0, 10.0]
@@ -762,30 +854,31 @@
pi = NArray.float(6).indgen!*2+10
qi = NArray.float(6).indgen!*2
vpi = VArray.new( pi, {"units"=>"mm"}, "p")
vqi = VArray.new( qi, {"units"=>"mm"}, "q")
- ai = NArray[2.0, 4.0]
+ ai = NArray[2.0, 2.5]
vai = VArray.new( ai ).rename("a")
#< test of interpolate >
gxi = vxi.to_g1D
gyi = vyi.to_g1D
gp = GPhys.new(xygrid,vp)
gq = GPhys.new(xygrid,vq)
+ #GPhys::extrapolation = true
gi = gd.interpolate(vxi,vyi,{"z"=>0.5})
GGraph::tone gi,true,"color_bar"=>true
###gd.interpolate(vxi,vyi,vr,vz) # nust fail by over-determination
gi = gd.interpolate([vxi,vyi])
#p gi.max, gi.min
- GGraph::tone gd,true,"min"=>-1.2,"max"=>1.2,"int"=>0.1
+ GGraph::tone gd[false,0],true,"min"=>-1.2,"max"=>1.2,"int"=>0.1
GGraph::scatter gxi, gyi, false,"type"=>4,"size"=>0.027,"index"=>3
- GGraph::color_scatter gxi, gyi, gi, false,"min"=>-1.2,"max"=>1.2,"int"=>0.1,"type"=>10,"size"=>0.029
+ GGraph::color_scatter gxi, gyi, gi[true,0], false,"min"=>-1.2,"max"=>1.2,"int"=>0.1,"type"=>10,"size"=>0.029
GGraph::color_bar
gi = gd.interpolate(vyi,vxi)
GGraph::tone gi,true,"color_bar"=>true
@@ -831,24 +924,55 @@
gi = gd.interpolate(vqi,vpi,{"z"=>0.5})
GGraph::tone gi,true,"color_bar"=>true
###gd.interpolate(vpi) # must fail by insufficient specification
-=begin
-=end
mask=d.lt(0.7)
missv = -999.0
d[mask.not] = missv
- p d[false,0]
+ #p d[false,0]
dm = NArrayMiss.to_nam(d, mask )
vdm = VArray.new( dm, {"missing_value"=>NArray[missv]}, "d")
gdm = GPhys.new(xyzgrid, vdm)
gi = gdm.interpolate(vpi,vqi)
# gi = gdm.interpolate(vxi,vyi)
GGraph::tone gi,true,"color_bar"=>true
+
+ print "start the test of coord_data_reverse\n"
+ GGraph::tone gp
+ GGraph::color_bar
+ xp = gp.coord_data_reverse("x",NArray.sfloat(30).indgen! )
+ GGraph::tone xp,true,"title"=>"#{xp.name} <-coord_data_reverse"
+ GGraph::color_bar
+ #< test of mouse_cut / dcl_fig_cut >
+ gdc = gdd[1..-1,1..-1,true]
+
+=begin
+ GGraph.set_fig "itr"=>4
+ GGraph.tone gdc,true,"int"=>0.15,"min"=>-1.05,"max"=>1.05
+ gpcut, line = gdc.mouse_cut(0,1,3)
+ x = gpcut.coord("x")
+ y = gpcut.coord("y")
+ DCL.sgpmzu(x.val,y.val,2,1,0.015)
+ GGraph.set_fig "itr"=>1
+ GGraph.tone gpcut,true,"int"=>0.15,"min"=>-1.05,"max"=>1.05
+ GGraph.color_bar
+=end
+
+ GGraph.set_fig "itr"=>4
+ GGraph.tone gdc,true,"int"=>0.15,"min"=>-1.05,"max"=>1.05
+ gpcut = gdc.dcl_fig_cut(0,1,[1.0,3.0,6.0],[1.1,2.2,5.5])
+ x = gpcut.coord("x")
+ y = gpcut.coord("y")
+ DCL.sgplu(x.val,y.val)
+ DCL.sgpmzu(x.val,y.val,2,1,0.015)
+ GGraph.set_fig "itr"=>1
+ GGraph.tone gpcut,true,"int"=>0.15,"min"=>-1.05,"max"=>1.05
+ GGraph.color_bar
+
#< finish >
DCL.grcls
end