Python map_coordinates() equivalent in Julia

What would be julia equivalent of that python code ? I am using map_coordinates — SciPy v1.14.0 Manual function. Here vals_r ,r_grid_stream_pixand theta_grid_stream_pix are of type Matrix{Float64}.

vals_r = map_coordinates(vals_r, (r_grid_stream_pix, theta_grid_stream_pix),
            order=1, cval=np.nan)

I couldn’t understand solution given in topic Equivalent of map_coordinates in Julia.

For a Matrix{T} as input array, the function map_coordinates is defined as follows:

using Interpolations
import StaticArrays: SVector
function map_coordinates(input::Matrix{T}, coordinates::Vector{SVector{N, T}}; method=BSpline(Linear())) where {N,T}
    output = T[]
    itp = interpolate(input, method, OnCell())
    for coord in coordinates
        push!(output, itp[coord...])
    end
    output
end

input is a Matrix{T} of size (m,n), and its elements
are the the values of a real function, f, defined on the planar rectangle [1, m] x [1,n], at the integer coordinates (i, j), i=1…m, j=1, …n.
i.e. input[i,j]=f(i, j).
coordinates is a vector of planar coordinates (x,y), with x\in[1,m], y\in [1, n].
map_coordinates evaluates the function f at these coordinates and pushes it values to output.

const Sv=SVector
input = [1.65 -1.35 -0.98;
           2.4  1.2 0.7;
           4.3  3.36 2.54;
          3.2   2.86  1.74]

Let us get the values at the indices for second column

julia> map_coordinates(input, [Sv(1.0,2), Sv(2,2.0), Sv(3,2), Sv(4,2)])
4-element Vector{Float64}:
 -1.35
  1.2
  3.36
  2.86

Retrieve the values at coordinates inside the rectangle on which the function is defined:

julia> coordinates=[Sv(1.1, 1.5), Sv(1.2, 1.89), Sv(3,2)]
julia> map_coordinates(input, coordinates)
3-element Vector{Float64}:
  0.31500000000000006
 -0.5495999999999999
  3.36

If at least a coordinate is outside this rectangle we get an error. The scipy function assigns a prescribed value, cval=0, to such coordinates. (we can add this to the above definition) .

julia>  map_coordinates(input, [Sv(0.97, 1.2)])
ERROR: BoundsError: attempt to access 4×3 interpolate(::Matrix{Float64}, BSpline(Linear())) with element type Float64 at index [0.97, 1.2]
1 Like

I am getting error as shown below.

ERROR: LoadError: MethodError: no method matching (::var"#map_coordinates#3"{var"#map_coordinates#1#4"})(::Matrix{Float64}, ::Tuple{Matrix{Float64}, Matrix{Float64}}; order::Int64, cval::Float64)

Closest candidates are:
  (::var"#map_coordinates#3")(::Matrix{T}, ::Array{SVector{N, T}, 1}; method) where {N, T} got unsupported keyword arguments "order", "cval"
   @ Main ~/Fluxes.jl:257

Stacktrace:
 [1] main()
   @ Main ~/Fluxes.jl:269
 [2] top-level scope
   @ ~/Fluxes.jl:349
in expression starting at /home/raman/Fluxes.jl:349

How to include order and cval values in function?

The error is caused by cval. In my definition of ‘map_coordinates’ there is no kwarg, cval. You can modify the function code to assign to output a prescribed value, cval, corresponding to an outside coord.

1 Like

I tried adding keywords but i still gets error.

ERROR: LoadError: MethodError: no method matching (::var"#map_coordinates#3"{var"#map_coordinates#1#4"})(::Matrix{Float64}, ::Tuple{Matrix{Float64}, Matrix{Float64}})

Closest candidates are:
  (::var"#map_coordinates#3")(::Matrix{T}, ::Array{SVector{N, T}, 1}; method, order, cval) where {N, T}
   @ Main ~/Fluxes.jl:257

Stacktrace:
 [1] main()
   @ Main ~/Fluxes.jl:277
 [2] top-level scope
   @ ~/Fluxes.jl:357
in expression starting at /home/raman/Fluxes.jl:357

I also tried this modified code: :point_down:

    function map_coordinates(input::Matrix{T}, coordinates::Vector{SVector{N, T}}; method=BSpline(Linear()), order=1, cval=NaN) where {N,T}
        output = T[]
        itp = interpolate(input, method, OnCell())
        for coord in coordinates
            try
                push!(output, itp[coord...])
            catch e
                if isa(e, BoundsError)
                    push!(output, cval)
                else
                    rethrow(e)
                end
            end
        end
        output
    end

order is a keyword for the Python version. Here it isn’t necessary. Also BoundsError is thrown only when method=BSpline(Linear()). For method=BSpline(Cubic()) it associates an interpolated value:

output=map_coordinates(input, [Sv(0.89, 1.2)];method=BSpline(Cubic()))
1-element Vector{Float64}:
 0.7674784494767151

Hence you should modify the code according to how you want to set the output (no matter which interpolation method is used), corresponding to a coord, exterior to the rectangle of definition.

@empet I replaced Cubic and also removed order=1, cval=nan but i still have error :yawning_face:

ERROR: LoadError: MethodError: no method matching (::var"#map_coordinates#3"{var"#map_coordinates#1#4"})(::Matrix{Float64}, ::Tuple{Matrix{Float64}, Matrix{Float64}})

Closest candidates are:
  (::var"#map_coordinates#3")(::Matrix{T}, ::Array{SVector{N, T}, 1}; method) where {N, T}
   @ Main ~/Fluxes.jl:257

Stacktrace:
 [1] main()
   @ Main ~/Fluxes.jl:277
 [2] top-level scope
   @ ~/Fluxes.jl:357
in expression starting at /home/raman/Fluxes.jl:357

I feel like you’ve been pointed to this in almost every single thread you’ve ever opened on this forum, but I’ll still do it again:

please read it carefully and take it to heart. Make sure that you post minimial, reproducible code snippets that yield the errors you are asking about. @empet has demonstrated how to use the function successfully, you are clearly running something else (the error comes from a main() function you are running, and shows that you are calling a function with a tuple of matrices as second argument when a vector of SVectors is required per the method definition.

Just blindly making superficial changes suggested by people on here into some different bit of code that you’re running but not showing really isn’t a good use of time for anyone.

7 Likes

Input variables are of type Matrix{Float64}, but i gets error.

using Interpolations, StaticArrays
vals_r= rand(288,258)
r_grid_stream_pix= rand(1000,50)
theta_grid_stream_pix= rand(1000,50)

function map_coordinates(input::Matrix{T}, coordinates::Vector{SVector{N, T}}; method=BSpline(Cubic())) where {N,T}
        output = T[]
        itp = interpolate(input, method, OnCell())
        for coord in coordinates
            try
                push!(output, itp[coord...])
            catch e
                if isa(e, BoundsError)
                    push!(output, cval)
                else
                    rethrow(e)
                end
            end
        end
        output
    end
    vals_r = map_coordinates(vals_r, (r_grid_stream_pix, theta_grid_stream_pix))

I am getting error :point_right:

ERROR: MethodError: no method matching map_coordinates(::Matrix{Float64}, ::Tuple{Matrix{Float64}, Matrix{Float64}})

Closest candidates are:
  map_coordinates(::Matrix{T}, ::Array{SVector{N, T}, 1}; method) where {N, T}
   @ Main REPL[5]:1

Stacktrace:
 [1] top-level scope
   @ REPL[6]:1

You get that error because instead of coordinates which must be a Vector of of SVector(x,y), you pass a tuple of two matrices:

(r_grid_stream_pix, theta_grid_stream_pix)

Look at the examples given in my first answer to see how the coordinates are defined.

2 Likes

I understand that it may be confusing as i am working on a lengthy script. My input variables are defined like this:

    theta_grid_stream_pix = ((theta_grid_stream_coord .+ theta[1]) ./ (2.0 * π .+ 2.0 * theta[1])) .* (2 * ntheta + 1)

    r_grid_stream_pix = similar(r_grid_stream_coord)
    for (i, j) in Tuple.(CartesianIndices(r_grid_stream_coord))
        r_val = r_grid_stream_coord[i, j]
        index = sum(r .< r_val) - 1
        if index < 1
            r_grid_stream_pix[i, j] = -1
        elseif index < nr - 1
            r_grid_stream_pix[i, j] = index + (r_val - r[index]) / (r[index + 1] - r[index])
        else
            r_grid_stream_pix[i, j] = nr
        end
    end

What should i do make them vector instead of Tuple as required by function? See attached file for detailed code. :point_down:
Fluxes.jl (11.7 KB)

As long as we don’t know what problem is to be solved through the posted code it is difficult to suggest how to transform your data.

map_coordinates is used in image processing for resampling images at non-integer position, i.e. for reconstruction of
continuous intensity image from a discrete one.
If the input array is interpreted as an image, then at each pixel (i,j) is stored its intensity, input[i,j].
Coordinates is vector of non-integer positions where the intensity is to be evaluated through interpolation.

If your component matrices, in the tuple you posted, represent respectively the first and the second coordinate, then you should flatten them.
If A is a matrix, vec(A) is its flattened version. Now if x is the vector representation of the first matrix, and y of the second one, then
coordinates=[SVector(a, b) for (a,b) in zip(x,y)].

@empet This folder contains library needed by my Julia script
Fluxes.jl (11.7 KB).
Please run it by modifying directory location to understand my problem. I want to get python map_coordinates — SciPy v1.14.0 Manual function equivalent in Julia. I tried as far as i know with no success.

The python version of map_coordinates works with input matrix of indices 0,1, \ldots, m-1, respectively 0,1\ldots,n-1.
Hence the function f is evaluated at points within the rectangle [0,m-1]\times[0, n-1].

The coordinates at which the function f
is to be evaluated are given in an odd way: coordinates = np.array([[x1, x2], [y1, y2]]), i.e. in the first row are x-coords and in the second one the y-coords.

Since Julia matrix indexing starts from 1, the function f is defined on [1,m]\times [1,n]
Hence coordinates =np.array([[x1, x2], [y1, y2]) must be converted to Julia
coordinates =[Sv(x1+1, y1+1), Sv(x2+1, y2+1)]
Let us check it. Take the second example from scipy manual:

from scipy.ndimage import map_coordinates
import numpy as np
a = np.arange(12.).reshape((4, 3))
inds = np.array([[0.5, 2], [0.5, 4]])
map_coordinates(a, inds, order=1, cval=-33.3)
array([  2. , -33.3])

Note that here the function f is evaluated at the points
(x,y)=(0.5, 0.5)\in [0, 2], and (x,y)=(2,4), with x=2\in [0,2], but y=4 \notin [0,2].
That’s why the value at the last point is set on cval=-33.3.

JUlia corresponding code:

const Sv=SVector
a=[0.0 1 2;
   3 4 5;
   6 7 8]

inds = [Sv(0.5+1, 0.5+1), Sv(2.0+1, 4+1)] 
map_coordinates(a, inds;  cval=-33.3) # map_coordinates is the Julia version that tests for BoundsError
2-element Vector{Float64}:
   2.0
 -33.3

The kwarg order from python version is replaced by method in Julia.

1 Like

If you are interested only in translating to Julia, a Python code involving map_coordinates then it’s better to define
the input matrix as an OffsetMatrix:

function map_coordinates(input::OffsetMatrix{T, Matrix{T}}, coordinates::Vector{SVector{2, T}}; 
        method=BSpline(Linear()),  cval=0.0) where T<:Real

and call it in the case of the given example from scipy manual, as follows:

using OffsetArrays
const Sv=SVector

a=[0.0 1 2;
3 4 5;
6 7 8]
m, n = size(a)
inds = [Sv(0.5, 0.5), Sv(2.0, 4)]
map_coordinates(OffsetArray(a, 0:m-1, 0:n-1), inds;  cval=-33.3)  
2-element Vector{Float64}:
   2.0
 -33.3

With input matrix defined as an OffsetMatrix with indices from 0, it works similar to python, without adding 1 to each coordinate, because the rectangle of definition of the function resulted from interpolation is [0,m-1]\times[0,n-1].

@empet I am still getting error:

ERROR: LoadError: MethodError: no method matching (::var"#map_coordinates#3"{var"#map_coordinates#1#4"})(::Matrix{Float64}, ::Tuple{Matrix{Float64}, Matrix{Float64}})
Stacktrace:
 [1] main()
   @ Main ~/Fluxes.jl:278
 [2] top-level scope
   @ ~/Fluxes.jl:358
in expression starting at /home/raman/Fluxes.jl:358

and my modified code is :point_down:

#### Define map_coordinates() in julia
    function map_coordinates(input::OffsetMatrix{T, Matrix{T}}, coordinates::Vector{SVector{2, T}}; 
        method=BSpline(Linear()),  cval=0.0) where T<:Real
        output = T[]
        itp = interpolate(input, method, OnCell())
        for coord in coordinates
            try
                push!(output, itp[coord...])
            catch e
                if isa(e, BoundsError)
                    push!(output, cval)
                else
                    rethrow(e)
                end
            end
        end
        output
    end
    
    # Join vector data through boundaries
    vals_r = hcat(vals_r_left[:, 1:1], vals_r_right, vals_r_left[:, end:-1:1], vals_r_right[:, 1:1])
    print(typeof(theta_grid_stream_pix))
    vals_r = map_coordinates(vals_r, (r_grid_stream_pix, theta_grid_stream_pix))

And it is also still the same error message.
You pass as a second argument a Tuple but the method expects a Vector. That is the same error as 5 days ago. Not that emmets second argument inds is a vector of SVectors.
Your second argument – note the round brackets – (r_grid_stream_pix, theta_grid_stream_pix) is a Tuple

As mentioned above and in several other of your threads. Please put in a bit more effort than just copying code that is answered to just post the error message to et us continue completely coding your problem or to quote from above

Especially when it is still the same error message and that has been answered!

1 Like

I have tried many time wrapping those input variables with Vector() and collect() but still getting same Tuple error. These input variables are defined as:

    theta_grid_stream_coord = π .- atan.(x_grid_stream', -z_grid_stream')
    theta_grid_stream_pix = ((theta_grid_stream_coord .+ theta[1]) ./ (2.0 * π .+ 2.0 * theta[1])) .* (2 * ntheta + 1)

    r_grid_stream_pix = similar(r_grid_stream_coord)
    for (i, j) in Tuple.(CartesianIndices(r_grid_stream_coord))
        r_val = r_grid_stream_coord[i, j]
        index = sum(r .< r_val) - 1
        if index < 1
            r_grid_stream_pix[i, j] = -1
        elseif index < nr - 1
            r_grid_stream_pix[i, j] = index + (r_val - r[index]) / (r[index + 1] - r[index])
        else
            r_grid_stream_pix[i, j] = nr
        end
    end

Please tell where is problem?

You always only post snippets so I can only guess.

Those lines generate the two variables that You put into a tuple when you call the function. But I can only guess, since you never post reproducible running code – I think I can not help you, because I am only repeating the requests everyone else does on – please read the answers, please make a reasonable effort yourself, please post reproducible examples that we can run.

Please see this Python map_coordinates() equivalent in Julia - #13 by raman_kumar it contains everything.