gphys_netcdf_io.rb in gphys-1.4.3

- old
+ new

@@ -523,11 +523,11 @@
       def __files2varray( files, varname, dim=nil )
 	if files.is_a?(NetCDF)
 	  # Single file. Returns a VArrayNetCDF. dim is ignored.
 	  file = files
 	  var = file.var( varname )
-	  raise "variable '#{varname}' not found in #{file}" if !var
+	  raise "variable '#{varname}' not found in #{file.inspect}" if !var
 	  VArrayNetCDF.new( var )
 	elsif files.is_a?(NArray)
 	  # Suppose that files is a NArray of NetCDF. Returns a VArrayCompsite.
 	  if dim.is_a?(Integer) && dim>=0 && dim<files.rank
 	    files = files[ *([0]*dim+[true]+[0]*(files.rank-dim-1)) ]
@@ -564,13 +564,34 @@
 	  raise TypeError, "argument files: not a NetCDF, String, NArray, Array, or Regexp"
 	end
 	[files, ncvar0]
       end
 
+      #--
+      #  re-arrange the array of files and check some inconsistencies
+      #
+      #  e.g.)
+      #     files = NArray[ file0.nc, file1.nc, file2.nc, file3.nc ]
+      #     The files have a variable, u(lat,lev,time),
+      #     which is tiled to four parts as the followings:
+      #                |    lat  |     lev    | time |
+      #       file0.nc | -90...0 | 1000...500 | 0..1 |
+      #       file1.nc |   0..90 | 1000...500 | 0..1 |
+      #       file2.nc | -90...0 |  500..100  | 0..1 |
+      #       file3.nc |   0..90 |  500..100  | 0..1 |
+      #
+      #     Then the output is
+      #        NArray[ [ [ file0.nc, file1.nc ],
+      #                  [ file2.nc, file3.nc ] ] ].
+      #
+      #  The order of the elements of the input array is arbitrary.
+      #
       def __files_dim_matching( files, varname )
 	# files: NArray of NetCDF
 
+	return files if files.length == 1
+
 	# < read the first file and check its rank >
 
 	file0 = files[0]
 	ncvar0 = file0.var(varname)
 	if files.rank > ncvar0.rank
@@ -579,47 +600,64 @@
 	convention = NetCDF_Conventions.find( ncvar0.file )
 	axposnames = convention::coord_var_names(ncvar0)
 
 	#< find correspoding dimensions >
 
-	j2i = Array.new
+        rank = axposnames.length
+        cvals = Array.new(rank){ Array.new }
+        sign = Array.new(rank)
+        units = Array.new(rank)
+        files.each do |file|
+          axposnames.each_with_index do |axposname, i|
+            nvax = file.var(axposname)
+            raise "No coordinate variable: #{axposname}' not found" if !nvax
+            axpos = VArrayNetCDF.new( nvax )
+            axv0 = axpos.val[0]
+            if units[i]
+              un = axpos.units
+              axv0 = un.convert2(axv0,units[i]) unless un == units[i]
+            else # the first file
+              units[i] = axpos.units
+              sign[i] = (axpos.val[-1] >= axv0 ? 1 : -1)
+            end
+            cvals[i].push axv0
+          end
+        end
+        # For the above example,
+        # cvals = [ [  -90,    0, -90,   0 ],
+        #           [ 1000, 1000, 500, 500 ],
+        #           [    0,    0,   0,   0 ] ]
+        # sign =  [ 1, -1, 1 ]
 
-	for ifl in 0...files.rank
-	  for jdt in 0...ncvar0.rank
-	    axpos_last = nil
-	    for kfl in 0...files.shape[ifl]
-	      nvax = files[*([0]*ifl+[kfl]+[0]*(files.rank-ifl-1))].var(axposnames[jdt])
-	      raise "No coordinate variable: #{axposname}' not found" if !nvax
-	      axpos = VArrayNetCDF.new( nvax )
-	      if kfl > 0
-		if axpos.length != axpos_last.length
-		  # not equal ==> this is the dimension looking for
-		  j2i[jdt] = ifl
-		  break
-		else
-		  axv0 = axpos.val[0]
-		  axv_last0 = axpos_last.val[0]
-		  if axv0 != axv_last0
-		    # not equal ==> this is the dimension looking for
-		    j2i[jdt] = ifl
-		    break
-		  end
-		end
-	      end
-	      axpos_last = axpos
-	    end
-	    break if j2i.include?(ifl)
-	  end
-	  if files.shape[ifl]>1 && !j2i.include?(ifl)
-	    raise "No dimension correpodence found for #{ifl}th dim" 
-	  end
-	end
-
-	for d in files.rank...ncvar0.rank
-	  j2i[ (j2i.index(nil) || files.rank) ] = d
-	  files.newdim!(d)
-	end
-	files = files.transpose(*j2i)
+        axs = Array.new(rank)
+        cvals.each_with_index do |cval, i|
+          axs[i] = cval.uniq
+          axs[i].sort!{|a,b| (a <=> b)*sign[i] }
+        end
+        # For the above example,
+        # axs   = [ [ -90, 0 ], [ 1000, 500 ], [ 0 ] ]
+        files_new = NArray.object( *axs.map{|ax| ax.length} )
+        # For the above example,
+        # shape of the files_new is [2, 2, 1]
+        files.length.times do |n|
+          file = files[n]
+          idx = Array.new(rank)
+          rank.times do |i|
+            axv0 = cvals[i][n]
+            idx[i] = axs[i].index(axv0)
+          end
+          # For the above example,
+          # idx =
+          #       [ 0, 0, 0 ] for n==0
+          #       [ 1, 0, 0 ] for n==1
+          #       [ 0, 1, 0 ] for n==2
+          #       [ 1, 1, 0 ] for n==3
+          files_new[*idx] = file
+        end
+        if files_new.eq(nil).count_true > 0
+          raise "No dimension correspondence found"
+        end
+        files = files_new
 
 	files
       end
 
       def __to_na_of_netcdf( files )