search

GenericMappingTools / GMT.jl

Generic Mapping Tools Library Wrapper for Julia
Other
194 stars 28 forks source link

ERROR: GMT: Failure to open virtual file #1558

Open stepanoslejsek opened 21 hours ago

stepanoslejsek commented 21 hours ago

I encountered the same error as in https://github.com/GenericMappingTools/GMT.jl/issues/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 commented 19 hours ago

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 commented 18 hours ago

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)