ERROR: GMT: Failure to open virtual file

[stepanoslejsek]

I encountered the same error as in #59 and https://discourse.julialang.org/t/gmt-error-message-failure-to-open-virtual-file/116843

However, after applying all provided solutions, the error persisted. The code that causes the error is following:

function select_points_near_tram_track(gt_dem, railway_map, distance_threshold=150)
  near_pts_idx = []
  railway_points = pointify(railway_map["route"].geometry)
  tram_route_matrix = hcat([p.coords.lon.val for p in railway_points], [p.coords.lat.val for p in railway_points])
  for i = 1:nrow(gt_dem)
    map_point = [gt_dem.geometry[i].coords.lon.val, gt_dem.geometry[i].coords.lat.val]
    distance_to_track = mapproject(map_point, dist2line=tram_route_matrix).data[3]
    if distance_to_track <= distance_threshold
      push!(near_pts_idx, i)
    end
  end
  return gt_dem[near_pts_idx, :]
end

It basically loops over a huge amount of points and computes a distance to a tram track. After running the loop for 100k times, the error is raised.

ERROR: GMT: Failure to open virtual file
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] GMTJL_Set_Object(API::Ptr{Nothing}, X::GMT.GMT_RESOURCE, ptr::Matrix{Float64}, pad::Int64)
    @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\gmt_main.jl:748
  [3] gmt(::String, ::Matrix{Float64}, ::Vararg{Matrix{Float64}})
    @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\gmt_main.jl:160
  [4] common_grd(::Dict{Symbol, Any}, ::String, ::Matrix{Float64}, ::Matrix{Float64})
    @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\common_options.jl:4063
  [5] common_grd(::Dict{Symbol, Any}, ::String, ::String, ::String, ::Matrix{Float64}, ::Vararg{Matrix{Float64}})
    @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\common_options.jl:4053
  [6] mapproject(cmd0::String, arg1::Vector{Float64}, arg2::Nothing; kwargs::@Kwargs{dist2line::Matrix{Float64}})
    @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\mapproject.jl:106
  [7] mapproject
    @ C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\mapproject.jl:73 [inlined]
  [8] mapproject
    @ C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\mapproject.jl:118 [inlined]
  [9] select_points_near_tram_track(gt_dem::GeoTable{PointSet{🌐, GeodeticLatLon{…}}, @NamedTuple{alt::Vector{…}}}, railway_map::Dict{String, GeoTable}, distance_threshold::Int64)   
    @ Main .\REPL[7]:8
 [10] select_points_near_tram_track(gt_dem::GeoTable{PointSet{🌐, GeodeticLatLon{WGS84Latest, Quantity{…}}}, @NamedTuple{alt::Vector{Union{…}}}}, railway_map::Dict{String, GeoTable})    @ Main .\REPL[7]:2
 [11] top-level scope
    @ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types.

I am using GMT.jl v1.20.0 from master branch, and I have following version

julia> GMT.GMTver
v"6.6.0"

but after first calling the mapproject function, I got warning about GMT version GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.

Here are all version of packages I am using:

(VehicleDataAnalysis) pkg> st
Status `C:\Users\oslejsek\AA4CC\vehicle-data-in-osm\VehicleDataAnalysis\Project.toml`
  [336ed68f] CSV v0.10.14
  [13f3f980] CairoMakie v0.12.14
  [a93c6f00] DataFrames v1.7.0
  [5752ebe1] GMT v1.20.0 `https://github.com/GenericMappingTools/GMT.jl.git#master`
  [f5a160d5] GeoIO v1.18.4
  [61d90e0f] GeoJSON v0.8.1
  [dcc97b0b] GeoStats v0.69.0
  [a98d9a8b] Interpolations v0.15.1
  [0f8b85d8] JSON3 v1.14.0
  [a90b1aa1] LibGEOS v0.9.2
  [e170d443] Tyler v0.2.0
  [ade2ca70] Dates

[joa-quim]

Hi,
I don't think this error has anything to do with the old #59. Now, the install process is very different.
But regarding the error, I need a way to reproduce it. Can you isolate the the first call of this line in a MWE?

distance_to_track = mapproject(map_point, dist2line=tram_route_matrix).data[3]

That said, your approach will be extremely slow. Calling mapproject a _huge number of times, once per point, _is not the way. The way would be to call it once with all the points.

Also, you don't need all that complication introduced by reading the data with GeoIO, GeoTables, etc. gmtread should do all that is needed with the plus that data comes already inside a GMTdaset.

But still furthermore. If you want to isolate the points that are at a certain distance to the railway, use the buffer function to obtain that polygon and use it with inpolygon to isolate those points.

[stepanoslejsek]

I knew that there should be more efficient way, how to do such thing. Since I am not that much familiar with GMT.jl. The reason that I use GeoTables is the ease of performing a lot of dataprocessing and geometry manipulation on files representing the tram track (also a .geojson files), since it is essentially an 'extended DataFrame'. Also, the next step is going to do some spatial interpolation, which is supported by GeoStats.jl. I'll consider using GMT.jl and GMTdatasets.

Regarding the MWE for the first call, I could reproduce it with the code from https://discourse.julialang.org/t/how-to-calculate-the-nearest-point-on-a-line-to-a-given-point/105743/6

using GMT
line = [146.54294815408372 -36.04985956921628; 146.5457684472937 -36.04924988190973; 146.54706548053542 -36.04879871026288; 146.54919201177802 -36.047176909745176; 146.55074543531066 -36.04585993436855]
pt2 = mapproject([146.54376257030356 -36.04616479173944], dist2line=line)

[stepanoslejsek]

When I tried buffer function, the following set of errors occurred. The whole code was:

julia> using GMT

julia> elev_map = gmtread("./map_data/elevation_map/zabaged_laz/zabaged_laz.geojson")
Vector{GMTdataset} with 3081064 segments
Showing first segment. To see other segments just type its element number. e.g. D[2]

Attribute table
┌─────────┬────────┬────────────┐
│     Row │ alt    │ Feature_ID │
├─────────┼────────┼────────────┤
│       1 │ 274.03 │ 1          │
│       2 │ 274.01 │ 2          │
│       3 │ 273.27 │ 3          │
│       4 │ 273.15 │ 4          │
│       5 │ 273.03 │ 5          │
│    ⋮    │   ⋮    │     ⋮      │
│ 3081060 │ 226.75 │ 3081060    │
│ 3081061 │ 227.47 │ 3081061    │
│ 3081062 │ 226.39 │ 3081062    │
│ 3081063 │ 227.29 │ 3081063    │
│ 3081064 │ 227.29 │ 3081064    │
└─────────┴────────┴────────────┘
             3081054 rows omitted
BoundingBox: [18.14125131335861, 18.14125131335861, 49.80814636486823, 49.80814636486823]
Global BoundingBox: [18.14125131335861, 18.254640543109993, 49.76887306370226, 49.843758859141765, 0.0, 0.0]
PROJ: +proj=longlat +datum=WGS84 +no_defs
WKT: GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

1×2 GMTdataset{Float64, 2}
 Row │     Lon      Lat 
─────┼──────────────────
   1 │ 18.1413  49.8081

julia> rail = gmtread("./map_data/railway_map/railway_tram7.geojson")
Vector{GMTdataset} with 281 segments
Showing first segment. To see other segments just type its element number. e.g. D[2]

Attribute table
┌─────┬────────┬──────────────┬──────────────────┬────────┬─────────┬─────────────────────────────────┬───────────────────┬──────────────────────────────────────┬──────────────┬─────│ Row │ oneway │ crossing_ref │ public_transport │ bridge │ covered │ network:wikipedia               │ operator:wikidata │ name                                 │ electrified  │ op ⋯├─────┼────────┼──────────────┼──────────────────┼────────┼─────────┼─────────────────────────────────┼───────────────────┼──────────────────────────────────────┼──────────────┼─────│   1 │        │              │                  │        │         │ cs:Tramvajová doprava v Ostravě │ Q11323689         │ Tram 7 Poruba, Vřesinská - Výškovice │              │ cs ⋯│   2 │        │              │                  │        │         │ cs:Tramvajová doprava v Ostravě │ Q11323689         │ Tram 7 Poruba, Vřesinská - Výškovice │              │ cs ⋯│   3 │        │              │                  │        │         │ cs:Tramvajová doprava v Ostravě │ Q11323689         │ Tram 7 Výškovice - Poruba, Vřesinská │              │ cs ⋯│   4 │        │              │                  │        │         │ cs:Tramvajová doprava v Ostravě │ Q11323689         │ Tram 7 Výškovice - Poruba, Vřesinská │              │ cs ⋯│   5 │        │              │ platform         │        │         │                                 │                   │                                      │              │    ⋯│  ⋮  │   ⋮    │      ⋮       │        ⋮         │   ⋮    │    ⋮    │                ⋮                │         ⋮         │                  ⋮                   │      ⋮       │    ⋱│ 277 │        │              │                  │        │         │                                 │                   │                                      │              │    ⋯│ 278 │        │              │                  │        │         │                                 │                   │                                      │              │    ⋯│ 279 │        │              │                  │        │         │                                 │                   │                                      │              │    ⋯│ 280 │        │              │                  │        │         │                                 │                   │                                      │              │    ⋯│ 281 │        │              │                  │        │         │                                 │                   │                                      │              │    ⋯└─────┴────────┴──────────────┴──────────────────┴────────┴─────────┴─────────────────────────────────┴───────────────────┴──────────────────────────────────────┴──────────────┴─────                                                                                                                                                       31 columns and 271 rows omittedBoundingBox: [18.158548, 18.24835, 49.7725848, 49.8348159]
Global BoundingBox: [18.1582647, 18.248401, 49.7725848, 49.834852, 0.0, 0.0]
PROJ: +proj=longlat +datum=WGS84 +no_defs
WKT: GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

452×2 GMTdataset{Float64, 2}
 Row │     Lon      Lat 
─────┼──────────────────
   1 │ 18.1585  49.8222
   2 │ 18.1586  49.8222
   3 │ 18.1587  49.8222
   4 │ 18.1588  49.8222
   5 │ 18.1589  49.8222
   6 │ 18.159   49.8222
  ⋮  │    ⋮        ⋮
 448 │ 18.2248  49.7727
 449 │ 18.2248  49.7727
 450 │ 18.2249  49.7727
 451 │ 18.2249  49.7727
 452 │ 18.2249  49.7728
        441 rows omitted

julia> buffer(rail, 150)
ERROR 1: IllegalArgumentException: point array must contain 0 or >1 elements

ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryType'.

ERROR 1: IllegalArgumentException: point array must contain 0 or >1 elements

ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryType'.


ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryType'.


ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetCoordinateDimension'.

ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetPointCount'.

ERROR 10: Pointer 'hGeom' is NULL in 'OGR_G_GetGeometryCount'.

String[]

Is there something I'm missing?

[joa-quim]

Vector{GMTdataset} with 3081064 segments

Wow, that's a lot of segments, and each with only 1 point. It would be much better if you had a single segments with 308104 points (a multi-point geometry). I deduce from the file name, zabaged_laz.geojson that this is LAZ data. If you have that data in LAZ format you can read it too with gmtread.

The problem in buffer(rail, 150) is that the buffer function is a cartesian one and you are passing the coordinates in geographical and the threshold in meters. For geographical data, use the buffergeo() function. e.,g.:

buffergeo(rail, width=150)

And still about the previous error with mapproject. I would like to understand why it failed. Could you create the MWE that I asked?

[stepanoslejsek]

The problem in buffer(rail, 150) is that the buffer function is a cartesian one and you are passing the coordinates in geographical and the threshold in meters. For geographical data, use the buffergeo() function. e.,g.:

I provided the MWE here: #1558 (comment). At the first call, this warning occurs GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6. After a_big_number_of_calls of mapproject it crashes.

[joa-quim]

Sorry, I got an email notification for your second post, not the one referring the MWE. If it crashes only after a long number of calls, that makes it much more difficult to catch. I would need the full example to try to reproduce it (the GMT [WARNING] is harmless).

Note that you can do a lot of processing with GMT alone, namely grid interpolations. I would be curious to know about typical things that you can do with Geostats that can't be done with GMT. And, BTW, GMT.jl has DataFrames as an extension. Having loaded DataFrames you can convert a GMTdataset into a DF with (sorry, forgot to export that function)
GMT.ds2df(D and vice versa with df2ds. This later one is limited in the sense that only one text column can be retained (but it keeps all the numeric ones)

[stepanoslejsek]

The buffergeo(rail, width=150) does not really work here, causing the following error:

julia> buffergeo(rail, width=150)
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
ERROR: BoundsError: attempt to access 0×0 Matrix{Float64} at index [1:0, 3]
Stacktrace:
 [1] throw_boundserror(A::Matrix{Float64}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Int64})
   @ Base .\abstractarray.jl:737
 [2] checkbounds
   @ .\abstractarray.jl:702 [inlined]
 [3] view(::Matrix{Float64}, ::Function, ::Int64)
   @ Base .\subarray.jl:184
 [4] buffergeo(line::Matrix{…}; width::Int64, unit::Symbol, np::Int64, flatstart::Bool, flatend::Bool, proj::String, epsg::Int64, tol::Float64)
   @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\proj_utils.jl:352
 [5] buffergeo(D::Vector{…}; width::Int64, unit::Symbol, np::Int64, flatstart::Bool, flatend::Bool, epsg::Int64, tol::Float64)
   @ GMT C:\Users\oslejsek\.julia\packages\GMT\yRgLc\src\proj_utils.jl:312
 [6] top-level scope
   @ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

Regarding the GMT.ds2df function, it cannot transform vector of datasets (which I have in my case)

[joa-quim]

ERROR: BoundsError: attempt to access 0×0 Matrix{Float64} at index [1:0, 3]

Well, from the error message it is saying that the rail variable is empty.
If I take the printed numbers from one of your posts above, it works fine

bf = buffergeo([18.1585  49.8222; 18.1586  49.8222; 18.1587  49.8222; 18.1588  49.8222; 18.1589  49.8222; 18.159   49.8222; 18.2248  49.7727; 18.2248  49.7727; 18.2249  49.7727; 18.2249  49.7727; 18.2249  49.7728], width=150)
BoundingBox: [18.15641713521431, 18.2269821907018, 49.77135184822378, 49.82354860980471]

60×2 GMTdataset{Float64, 2}
 Row │   col.1    col.2
─────┼──────────────────
   1 │ 18.227   49.7727
   2 │ 18.227   49.7725
   3 │ 18.2269  49.7723
   4 │ 18.2268  49.7721
...
viz(bf)

[joa-quim]

Can you make that railway_tram7.geojson file available?

[stepanoslejsek]

ERROR: BoundsError: attempt to access 0×0 Matrix{Float64} at index [1:0, 3]

Well, from the error message it is saying that the rail variable is empty. If I take the printed numbers from one of your posts above, it works fine

bf = buffergeo([18.1585  49.8222; 18.1586  49.8222; 18.1587  49.8222; 18.1588  49.8222; 18.1589  49.8222; 18.159   49.8222; 18.2248  49.7727; 18.2248  49.7727; 18.2249  49.7727; 18.2249  49.7727; 18.2249  49.7728], width=150)
BoundingBox: [18.15641713521431, 18.2269821907018, 49.77135184822378, 49.82354860980471]

60×2 GMTdataset{Float64, 2}
 Row │   col.1    col.2
─────┼──────────────────
   1 │ 18.227   49.7727
   2 │ 18.227   49.7725
   3 │ 18.2269  49.7723
   4 │ 18.2268  49.7721
...
viz(bf)

Actually that example is taken from https://discourse.julialang.org/t/how-to-calculate-the-nearest-point-on-a-line-to-a-given-point/105743/6. The railway_tram7.geojson is exported from https://overpass-turbo.eu/s/1SGM

[joa-quim]

But can't you zip it and attached it here?

[stepanoslejsek]

There is the railway_tram7.geojson file
railway_tram7.zip

[joa-quim]

Thanks. I can now reproduce the error. And when I try with the cartesian version I get further info

julia> buffer(rail, 0.00150)
ERROR 1: IllegalArgumentException: point array must contain 0 or >1 elements

So it looks like one, or more, of the 281 segments are empty. I'll try to make the reading and/or buffering functions more resistant to faulty data. But meanwhile this works fine.

bf = buffergeo(rail[1], width=150)

[joa-quim]

OK, the problem is that the file has two lines, one for each rail, plus a lot of other single point segments and it were these single point segments that resulted caused the error in buffergeo. I've fixed it but that doesn't help you. The best is that you use the first rail segment as in my example above.

[joa-quim]

Can we close this?