mapproject(): "GMT: Failure to open virtual file"

[droodman]

I'm getting the above error after I call mapproject() exactly 250,000 times in a row. I would like to call it many more times!

If I call it again from the REPL after my program crashes, it gives the same error. I have to restart Julia to get rid of the message.

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 16 × Apple M4 Max
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 12 default, 0 interactive, 6 GC (on 12 virtual cores)
Environment:
  JULIA_NUM_THREADS = auto
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1

This happens in both GMT.jl 1.31.1 and 1.32.0.

[joa-quim]

Can you show the command used?

Note that there is an overhead of 1e-3 to 1e-4 sec to call the GMT C lib via the gmt(...) command, so doing that large number of calls will be slow anyway. The advised way is to do a low number of calls but with larger amount of data.

Using 1.32.0 will not change nothing in this regard.

[droodman]

The exact line is

distance = mapreduce(g->mapproject([long lat], dist2line=g)[3], min, gcountries[(row.ISO2,)].geometry_GMT)

where long and lat are Float64 and the gcountries[(row.ISO2,)].geometry_GMT items are constructed with logic like

ls2mat(ls) = [fn(pt) for pt ∈ GI.getpoint(ls), fn ∈ (GI.x, GI.y)]  # convert linestring to nx2 matrix
[mat2ds(ls2mat(ls), wkt="EPSG:4326") for ls ∈ GI.getgeom(mls)]  # represent a multilinestring as a vector of these objects

And mls is a MultiLineString. In other words, I have converted MultiLineStrings into vectors of outputs from mat2ds(). Then given long and lat, I'm finding the minimum distance to that MultiLineString, which actually represents a country's fiber optic network.

Advice welcome on how to do this efficiently. I'm running a Monte Carlo simulation in which I perturb some 4500 coordinates and recompute their distances from the fiber network each time. I figure I can introduce Distributed to save some time.

[joa-quim]

I would need something that I can try, but even then is not assured. That error smells OS related. Note that the input to mapproject can be a simple matrix. It doesn't need to be a GMTdataset (EPSG:4326 is the default expected by mapproject)

And mls is a MultiLineString

Where do you get those? Can't you read the data with GMT into GMTdatasets and thus avoid all those calls to GDAL? (each getgeom does one GDAL call).

[joa-quim]

Note, resetGMT() should avoid having to kill the Julia session.

[droodman]

I'm starting with the AfTerFibre shapefile representing the fiber network in Africa. For each of the 8 African countries I'm working with, I intersect that continent-scale map with the country's national territory, as represented by another shapefile. Then, given a coordinate from a given country, I only need to compute the minimum distance to the part of the network in that country, which should save time. Actually, if there's a fiber line just over the border, that counts too when determining nearness to the network. So I form 0.2-degree buffers around each country before intersecting with the network map. The result is one MultiLineString for each of the 8 countries.

Can GMT do those calculations? I have been using GeoInterface, GeometryOps, and GeoDataFrames for this initial processing, then converting to GMT.

Anyway the time cost of the initial processing is small relative to the cost of a ~million calls to mapproject. Tomorrow, I can pull together and share the actual example, with data sets.

The advised way is to do a low number of calls but with larger amount of data.

How would I do that here? Can I make mapproject()'s first argument an nx2 matrix of coordinates?

[joa-quim]

Can GMT do those calculations? I have been using GeoInterface, GeometryOps, and GeoDataFrames for this initial processing, then converting to GMT.

Yes. It has buffers, intersections, etc... Actually it has access to the entire GDAL. But if you have that MultiLineString you should be able to convert it to a vector of GMTdasets, which mapproject can ingest. It doesn't need to be on point-by-point basis. It should be able to do it all at once. OFC, this needs to be put to test.

Can I make mapproject()'s first argument an nx2 matrix of coordinates?

Absolutely.

country's national territory, as represented by another shapefile.

D = coast(DCW=:ZA, dump=true), gets you South Africa. Same for all world countries (pass in the 2 char country code)

[droodman]

OK, here is an archive with data and code. On my machine it crashes after ~10 minutes. And if I check the value of the `calls` in the REPL after the crash I get exactly 250,000. https://openphilanthropy-my.sharepoint.com/:u:/g/personal/david_roodman_openphilanthropy_onmicrosoft_com/Ed04dYF-XLVBlI2Z7-7IWJYBa9AFvtdrpy0xiFddwMdBjQ?e=ohl8iO

[joa-quim]

OK, I got it. Will try later this afternoon on Windows. BTW, what OS do you use?

[joa-quim]

And if I check the value of the calls in the REPL after the crash I get exactly 250,000.

And memory was exhausted? I mean, could it e due to a memory leak?

[droodman]

I'm in macOS. I can try Windows and Linux under Windows later today. I haven't checked memory usage just because to stop exactly at 250,000 for that reason would be a coincidence.

[droodman]

Just ran it again: No memory pressure on my Mac.

[joa-quim]

Findings.

First, this, epsg=4326, is wrong and caused me many GDAL errors. If you want to pass an EPSG it must be an integer, like: epsg=4326

I can reproduce the error on Windows after ~3 min and, curiously, calls = 249998. So clearly a recursive limit in C compilers is being hit. I have a vague idea of the problem but to solve it fully I need to find out exactly what limit is being reached. The error shows up in GMT.jl because all C states (from the GMT C lib) are being kept in a struct that is kept in a memory shared with Julia (The G_API[1] pointer). So, to workaround the problem for now what you can do is this:

global limits = 100000; insert before the loop starts
...
      global calls += length(gcountries[(row.ISO2,)].geometry_GMT)
      (calls > limit) && (global limit += 100000; resetGMT())

Note, I had to add the global due to the ultra-irritating Julia scope rules.

With that your example run in:

julia> @time include("c:/v/mapproject_crash/crash_demo.jl")
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
GMT [WARNING]: GMT_COMPATIBILITY: Expects values from 6 to 6; reset to 6.
444.479178 seconds (237.69 M allocations: 22.512 GiB, 0.57% gc time, 12.75% compilation time: 5% of which was recompilation)

You may not have that annoying warning as my GMT build is local and uses MSVC instead of cross compiled GCC.

Now, it looks that you are doing lots of unnecessary iterations. Namely your clusters have points in all countries, which you loop for each country (I think). Have a look at the function randinpolygon that let you generate points inside each country only. I'm 99.99% sure that you could do all of it with GMT only and thus simplify the procedure even further but that would require more work to test.

[joa-quim]

And, what does your system prints with? (or is this linux only?)

cat /proc/sys/fs/file-max

[joa-quim]

Could resist half-trying it. A small example for Namibia that I visited recently and got astonished by its beauties

Dna = filter(coast(DCW=:NA, dump=true), _unique=1);    # Get the largest polygon of Namibea

Dnet = gmtread("c:/v/mapproject_crash/terrestrial.shp");

Di = GMT.intersection(Dna, Dnet);       # Needed to prefix it with GMT because LibGEOS was loaded too.

viz(Dna, plot=(data=Di,lc=:red), proj=:guess)

Now, only generate points inside and compute distances.