Filter tcols in batchextract when some are not in the cube axes

src/Cubes/Cubes.jl

@@ -230,26 +230,30 @@ function batchextract(x,i)
     axinds = map(sch.names) do n
         findAxis(n,x)
     end
-
     tcols = columns(i)
     #Try to find a column denoting new axis name and values
     newaxcol = nothing
+
     if any(isnothing,axinds)
         allnothings = findall(isnothing,axinds)
         if length(allnothings) == 1
             newaxcol = allnothings[1]
         end
+        tcols = (;[p[1:2] for p in zip(keys(tcols), values(tcols), axinds) if !isnothing(last(p))]...)
         axinds = filter(!isnothing,axinds)
     end
+
     allax = 1:ndims(x)
     axrem = setdiff(allax,axinds)
     ai1, ai2 = extrema(axinds)
+
     if !all(diff(sort(collect(axinds))).==1)
         #Axes to be extracted from are not consecutive in cube -> permute
         p = [1:(ai1-1);collect(axinds);filter(!in(axinds),ai1:ai2);(ai2+1:ndims(x))]
         x_perm = permutedims(x,p)
         return batchextract(x_perm,i)
     end
+
     cartinds = map(axinds,tcols) do iax,col
         axcur = caxes(x)[iax]
         map(col) do val

test/Cubes/batchextraction.jl

@@ -50,4 +50,34 @@ r = c[othersites]
 @test r.point.values == [3,5]
 @test c[lon=33.0,time=Date(2000,7,1)][:] == r[point=5][:]
 
+end
+
+@testitem "Batch extraction with site first" begin
+    # This should also be a test case for #194 because it is the same underlying cause.
+    using Dates
+    lons = range(30,35,length=10)
+    lats = range(50,55,step=0.5)
+    times = Date(2000,1,1):Month(1):Date(2000,12,31)
+
+    data = rand(length(lons),length(lats), length(times));
+
+    c = YAXArray([RangeAxis("longitude",lons),RangeAxis("latitude",lats),RangeAxis("time",times)],data)
+    c_perm = permutedims(c,(3,2,1))
+
+    sites_first = [(site=string(i), lon = rand()*5+30, lat = rand()*5+50) for i in 1:200]
+    sites_perm = [(lat = rand()*5+50, lon = rand()*5+30,site = string(i % 5)) for i in 1:200]
+
+    sites_pure = [(lon = n.lon, lat=n.lat) for n in sites_first]
+    lon,lat = sites_pure[10]
+    r = c[sites_first]
+    @test r isa YAXArray
+    @test YAXArrays.Cubes.axname.(caxes(r)) == ["site","time"]
+    @test r.site.values == string.(1:200)
+    @test all(isequal.(c[lon=lon,lat=lat][:], r[10,:]))
+
+    r = c_perm[sites_first]
+    @test r isa YAXArray
+    @test YAXArrays.Cubes.axname.(caxes(r)) == ["time","site"]
+    @test r.site.values == string.(1:200)
+    @test all(isequal.(c[lon=lon,lat=lat][:], r[:,10]))
 end