GeoData: "array has wrong ordering"

General Usage
,
1

Dear All,

I am a new user of Julia following this Example. I got an error using the below code, any suggestions or comments will be appreciated :

using GeophysicalModelGenerator, GeoDatasets,DelimitedFiles

Load data:

data=readdlm(“peakdelay_3_Degrees_Hz.txt”,‘,’,‘\n’, skipstart=0,header=false)
lon = data[:,1];
lat = data[:,2];
depth = 1*data[:,3];
peakdelay = data[:,4];

resolution = (length(unique(lon)), length(unique(lat)), length(unique(depth)))

Lon = reshape(lon, resolution);
Lat = reshape(lat, resolution);
Depth = reshape(depth, resolution);
peakdelay = reshape(peakdelay, resolution);

EQ_Data = GeoData(Lon,Lat,Depth,(Peakdelay=peakdelay,));

I am getting this error while running the last command line, I have tried to double check each variable :

ERROR: It appears that the lon array has a wrong ordering
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] GeoData(lon::Array{Float64, 3}, lat::Array{Float64, 3}, depth::Array{Float64, 3}, fields::@NamedTuple{Peakdelay::Array{Float64, 3}}, atts::Nothing)
@ GeophysicalModelGenerator ~/.julia/packages/GeophysicalModelGenerator/AiGOf/src/data_types.jl:175
[3] GeoData(lon::Array{Float64, 3}, lat::Array{Float64, 3}, depth::Array{Float64, 3}, fields::@NamedTuple{Peakdelay::Array{Float64, 3}})
@ GeophysicalModelGenerator ~/.julia/packages/GeophysicalModelGenerator/AiGOf/src/data_types.jl:158
[4] top-level scope
@ REPL[19]:1

The file input is in this format :

104.347,29.347,1500,0,0
104.347,29.347,433.333333333333,0,0
104.347,29.347,-633.333333333333,0,0
104.347,29.347,-1700,0,0
104.347,29.347,-2766.66666666667,0,0
104.347,29.347,-3833.33333333333,0,0
104.347,29.347,-4900,0,0
104.347,29.347,-5966.66666666667,0,0
104.347,29.347,-7033.33333333333,0,0
104.347,29.347,-8100,0,0
104.347,29.3705384615385,1500,0,0
104.347,29.3705384615385,433.333333333333,0,0
104.347,29.3705384615385,-633.333333333333,0,0
104.347,29.3705384615385,-1700,0,0
104.347,29.3705384615385,-2766.66666666667,0,0
104.347,29.3705384615385,-3833.33333333333,0,0
104.347,29.3705384615385,-4900,0,0

2

(I don’t personally use these packages, just trying to help.) That error comes from here in the source, specifically from line 175.

I’m trying to mentally parse that if condition, and it reads to me like it expects an array with a particular orientation, e.g. row-major vs column-major ordering, and will throw an error if the provided array is the opposite orientation.

1 Like
3

Can you perhaps upload the “peakdelay_3_Degrees_Hz.txt” file, such that I can reproduce it?
In this example, Lon, Lat,Depth should be 3D arrays with Lon changing in the first dimension as in:

julia> using GeophysicalModelGenerator
julia> Lon,Lat,Depth=XYZGrid(0:2,1:4,5:6);
julia> Lon
3×4×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0

[:, :, 2] =
 0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0

I suspect that your Lon changes in the second dimension, which would throw the error you see.

4

Thank you for your kind response. I have send you a file in person with some details.

5

ok, for general info, here some worked-out examples:
Correct

julia> Lon,Lat,Depth=XYZGrid(0:2,1:4,5:6);
julia> Lon
3×4×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0

[:, :, 2] =
 0.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0
 2.0  2.0  2.0  2.0
julia> GeoData(Lon,Lat,Depth,(;Depth))
GeoData 
  size      : (3, 4, 2)
  lon       ϵ [ 0.0 : 2.0]
  lat       ϵ [ 1.0 : 4.0]
  depth     ϵ [ 5.0 : 6.0]
  fields    : (:Depth,)
  attributes: ["note"]

Wrong:

julia> Lat1,Lon1,Depth1=XYZGrid(0:2,1:4,5:6);
julia> Lon1
3×4×2 Array{Float64, 3}:
[:, :, 1] =
 1.0  2.0  3.0  4.0
 1.0  2.0  3.0  4.0
 1.0  2.0  3.0  4.0

[:, :, 2] =
 1.0  2.0  3.0  4.0
 1.0  2.0  3.0  4.0
 1.0  2.0  3.0  4.0
julia> GeoData(Lon1,Lat1,Depth1,(;Depth))
ERROR: It appears that the lon array has a wrong ordering
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] GeoData(lon::Array{Float64, 3}, lat::Array{Float64, 3}, depth::Array{Float64, 3}, fields::@NamedTuple{Depth::Array{Float64, 3}}, atts::Nothing)
   @ GeophysicalModelGenerator ~/.julia/packages/GeophysicalModelGenerator/AiGOf/src/data_types.jl:175
 [3] GeoData(lon::Array{Float64, 3}, lat::Array{Float64, 3}, depth::Array{Float64, 3}, fields::@NamedTuple{Depth::Array{Float64, 3}})
   @ GeophysicalModelGenerator ~/.julia/packages/GeophysicalModelGenerator/AiGOf/src/data_types.jl:158
 [4] top-level scope
   @ REPL[88]:1

Fix:

julia> Lon2 = permutedims(Lon1,(2,1,3));
julia> Lat2 = permutedims(Lat1,(2,1,3));
julia> Depth2 = permutedims(Depth1,(2,1,3));

julia> GeoData(Lon2,Lat2,Depth2,(;Depth2))
GeoData 
  size      : (4, 3, 2)
  lon       ϵ [ 1.0 : 4.0]
  lat       ϵ [ 0.0 : 2.0]
  depth     ϵ [ 5.0 : 6.0]
  fields    : (:Depth2,)
  attributes: ["note"]
6

ok, in your specific example the order is a bit different (depth changes first), so one way to do this is:

using GeophysicalModelGenerator, DelimitedFiles

data=readdlm("1_test_peakdelay.txt",',','\n', skipstart=0,header=false);

lon = data[:,1];
lat = data[:,2];
depth = data[:,3];
peakdelay = data[:,4];
res = (length(unique(lon)), length(unique(lat)), length(unique(depth)))

Lon = zeros(res);
Lat = zeros(res);
Z = zeros(res);
Peakdelay = zeros(res);

num = 1
for i in 1:res[2], j in 1:res[1], k in 1:res[3]
    Lon[i,j,k] = lon[num];
    Lat[i,j,k] = lat[num];
    Z[i,j,k] = depth[num];
    Peakdelay[i,j,k] = peakdelay[num];
    num += 1
end 

EQ_Data = GeoData(Lon,Lat,Z,(; Peakdelay,));
7

or slightly shorter:

using GeophysicalModelGenerator, DelimitedFiles

data=readdlm("1_test_peakdelay.txt",',','\n', skipstart=0,header=false);

lon,lat,depth = data[:,1], data[:,2], data[:,3];
peakdelay = data[:,4];

res = (length(unique(lon)), length(unique(lat)), length(unique(depth)))
Lon = zeros(res);
Lat = zeros(res);
Z = zeros(res);
Peakdelay = zeros(res);

num = 1
for i in 1:res[2], j in 1:res[1], k in 1:res[3]
    I = CartesianIndex(i,j,k)
    Lon[I], Lat[I], Z[I] = lon[num], lat[num], depth[num];
    Peakdelay[I] = peakdelay[num];
    num += 1
end 

EQ_Data = GeoData(Lon,Lat,Z,(; Peakdelay,));
2 Likes