What is your favorite mapping package?

j_u · October 29, 2021, 2:00pm

It seems that I have only xdg:

whereis xdg
xdg: /etc/xdg

I am reading that “xdg-open is for use inside a desktop session only. It is not recommended to use xdg-open as root.” I am using remote machine accessing it via SSH and VS Code and that’s the way I am usually using Julia. I will try with GNOME desktop soon. Anyway, its not a big problem.

How can I direct GMT to save the file with the plot into desired location?

A.
In case of data = gmtread("/path/to/gpx"); I am getting:

ERROR: LoadError: Must select one input data type (grid, image, dataset, cmap or ps)

However, this syntax seems to work correctly with geojson, and semi correctly with kml. (As for geojson and kml please see the plots below as this is the same situation).

B.
In case of data = gmtread("/path/to/gpx", ogr=true); I am getting:

Vector{GMTdataset} with 0 segments

With geojson everything seems to work ok, however, 180 meridian seems to be not automaticly recognized. Kml also seems to work but the plot is only showing half of my data. Enclosed please find the files with plots.

geojson:

kml:

C.
What would be the correct syntax to prepare the plot with both netCDF and geojson (gpx) data combined?

joa-quim · October 29, 2021, 3:30pm

Use savefig="figname" Image format is picked from figname extension (default saves it in current dir)

The gmtread case. This function automatically detects many, raster and vector, formats. But not all. When it doesn’t you get that error message telling user that he/she must help a bit by specifying what type of data is expected to be read.

The kml and json cases. I would probably have to gave a look to the data (if we go into this detail level, please open an issue in GMT.jl so we don’t bother others with package imperfections). But generically, GMT knows how to wrap data around the the globe when it’s told that data is geographic and one asks for a projection. Otherwise, dateline limit comes into play. There are configuration option in GMT to tell were we break the sphere. Default is dateline but that can be set with the FORMAT_GEO_OUT parameter that one would set probably as par=(FORMAT_GEO_OUT="+D") option in imshow or plot modules.
Using projection=:guess often works pretty well and gets you a nice projected map (when the referencing system is known from input data).

The gpx case. Now here things look pretty more complicated. I made a couple of attempts and I’m puzzled because those files do not seem to be read very well by GDAL and violate the GDAL vector model. That is why you get that

Vector{GMTdataset} with 0 segments

For example, if you try with (unzip attached file)

julia> gmtread("25708.gpx", ogr=true)
Vector{GMTdataset} with 0 segments

But if you convert it to, for example, geopackage with

ogr2ogr("25708.gpx", dest="25708.gpkg")

Then we can do (again, this can be made more user friendly if I can teach imshow to recognize the file format)

D = gmtread("25708.gpkg")     # No need to tell it the format because ``.gpkg`` are automatically detected
imshow(D, proj=:guess, savefig="tricked_gpx.png")

We can continue this in another thread. Basically, use grdimage(), plot!(...), etc. Remember GMT builds maps by layers.

P.S. Sorry, Discourse doesn’t let me attach a .zip or .gps file

j_u · October 29, 2021, 7:36pm

Thanks a lot for this intro to GMT! I really appreciate the time you took to provide the advice.

Just one more question, as it is not fully clear to me. You wrote:

“If you have a GDAL built with expat (you do if you installed GMT.jl on Windows)”

and on GDAL webpage [GPX - GPS Exchange Format — GDAL documentation] it is written:

“OGR has support for GPX reading (if GDAL is build with expat library support) and writing.”

Do I have this read and write OGR GPX support if I am on Linux? It seems that ogr2ogr gpx to gpkg conversion is working, however, I wanted to double check, thus the question. Is there any (other) way to check it?

joa-quim · October 29, 2021, 8:17pm

You likely have the GDAL binaries as well so you should be able to do ogr2ogr --formats and see the list of available drivers.
With the wrapper I don’t know how to list them all but you can do this.

julia> GMT.Gdal.GDALGetDriverCount()
220

julia> GMT.Gdal.OGRGetDriverCount()
83

julia> getdriver("gpx")       # will be NULL if driver does not exist or the name is wrong
Driver: GPX/GPX

BTW, I think I fixed the gpx issue. If you update to GMT.jl#master the gpx format should be detected automatically and data loaded from it.

j_u · October 30, 2021, 12:17am

Thank you for the info.

I removed GMT v0.38.0 and added GMT.jl#master which resulted in GMT v0.38.0 https://github.com/GenericMappingTools/GMT.jl.git#master in my status.

On GMT v0.38.0 I was receiving:

julia> GMT.Gdal.GDALGetDriverCount()
0

julia> GMT.Gdal.OGRGetDriverCount()
0

julia> getdriver("gpx")
NULL Driver

On GMT v0.38.0 https://github.com/GenericMappingTools/GMT.jl.git#master I am receiving:

julia> GMT.Gdal.GDALGetDriverCount()
215

julia> GMT.Gdal.OGRGetDriverCount()
80

julia> getdriver("gpx")
Driver: GPX/GPX

Among impressive list of supported formats ogr2ogr --formats brings:

GPX -vector- (rw+v): GPX

so I guess GPX is supported.

However, when
data = gmtread("path/to/gpx", ogr=true)
I still see
Vector{GMTdataset} with 0 segments
and it is not possible to plot it due to
ERROR: Input has zero segments where it can't be.

I did some additional tests, however, I was not able to overcome the limit of datetime line. If you find it suitable, I think that it might be better if you take a look at the underling data by yourself. Below please find the links to gpx and the same data in geojson. Should I be of any help please let me know.

https://ibmb2.s3.eu-de.cloud-object-storage.appdomain.cloud/R1v1.geojson
https://ibmb2.s3.eu-de.cloud-object-storage.appdomain.cloud/R1v1.gpx

joa-quim · October 30, 2021, 12:53am

Hmm, I’m very surprised that you got

julia> GMT.Gdal.GDALGetDriverCount()
0

with 0.38.0 That is one of the CI tests that run on Linux and it always gets a (way) > 0 count. But as I said in another thread, can we move this to the GMT forum? We are now talking on package implementation/(bugs?) details.

j_u · October 30, 2021, 1:26am

I will double check. There is a lot of new things associated with the GMT for me.

Thank you. Please expect a message from me there, however, this might not be immediately. I would like to take some time to at least slightly more familiarize myself with the software. I also mentioned that I have some due things related to this forum. I would like to underline that I am really interested in the GMT package.

j_u · October 30, 2021, 10:34pm

Hi, as I informed I did additional testing. I will try to be concise.

GMT.Gdal.GDALGetDriverCount()
I can not exactly track the reason for julia> GMT.Gdal.GDALGetDriverCount() 0 that I indicated in my post https://discourse.julialang.org/t/what-is-your-favorite-mapping-package/70509/25

I did 2 additional new installs in completely new Julia and Conda environments.

This time I saw julia> GMT.Gdal.GDALGetDriverCount() 0 as well, however, it was on https://github.com/GenericMappingTools/GMT.jl.git#master. It was only once just after the package add. When I exited Julia and started again the number was 215. In the second environment, with GMT v0.38.0. I tried to replicate the situation and executed the command fresh after the package add and I saw 215. When / if I track the reason, I will let you know on GitHub. My current understanding is that all is working fine.

GPX

I can not make gpx work directly. When I follow your suggestion and use gpkg everything works fine (almost, with dateline exception), however, when I try to plot directly from gpx I still see:

Vector{GMTdataset} with 0 segments
Vector{GMTdataset} with 0 segments

When I convert my file to gpkg and set:
imshow(data, proj=:guess, savefig="file_gpkg.png")
I see

psxy [WARNING]: Conical projections with full 360 longitude range require the projection central longitude to be at the mid-point
psxy [ERROR]: Option -R parsing failure. Correct syntax:
        -R<xmin>/<xmax>/<ymin>/<ymax>[/<zmin>/<zmax>]
          Append +r if giving lower left and upper right coordinates
        -Rg or -Rd for global domain
        -R<grdfile> to take the domain from a grid file
psxy [ERROR]: Offending option -R-180/180/34/57
ERROR: LoadError: Something went wrong when calling the module. GMT error number = 72

When I use proj=:Mercator, it results in a plot, however, there is still a problem with dateline.

Note
As I indicated earlier, I am using Julia Version 1.8.0-DEV.829 (I am trying its new BLAS abilities) and PythonCall (due to some reasons it suites me better).

#master
I would like to ask if GMT v0.38.0 and GMT v0.38.0 https://github.com/GenericMappingTools/GMT.jl.git#master is the same thing? The numbers are the same. So far I mostly used registered packages. When I was using #master it happened that the new package version was different to the previous one. Now I see the same version numbers. I was on the GMT v0.38.0, and followed your suggestion to upgrade to #master thus the question.

joa-quim · October 30, 2021, 10:58pm

Confusing, GMT v0.38.0 https://github.com/GenericMappingTools/GMT.jl.git#master should be the master version. And with that it works for me. Note that your data is geographic so if you plot it linearly you see the split you showed before. I made a couple of attempts and I get.

I will tag a new version (will be 0.38.1). Keep and eye on the GMT.jl repo (maybe a 1/2 till it gets done).

julia> D = gdalread("R1v1.gpx");
Warning 1: The output driver does not natively support DateTime type for field time. Misconversion can happen. -mapFieldType can be used to control field type conversion.
Warning 6: Field time create as date field, though DateTime requested.
Warning 6: Normalized/laundered field name: 'geoidheight' to 'geoidheigh'
Warning 6: Normalized/laundered field name: 'ageofdgpsdata' to 'ageofdgpsd'
Warning 6: Normalized/laundered field name: 'route_point_id' to 'route_poin'
Warning 1: The output driver does not natively support DateTime type for field time. Misconversion can happen. -mapFieldType can be used to control field type conversion.
Warning 6: Field time create as date field, though DateTime requested.
Warning 6: Normalized/laundered field name: 'geoidheight' to 'geoidheigh'
Warning 6: Normalized/laundered field name: 'ageofdgpsdata' to 'ageofdgpsd'
Warning 6: Normalized/laundered field name: 'track_seg_id' to 'track_seg_'
Warning 6: Normalized/laundered field name: 'track_seg_point_id' to 'track_se_1'
Warning 1: The output driver does not natively support DateTime type for field time. Misconversion can happen. -mapFieldType can be used to control field type conversion.
Warning 6: Field time create as date field, though DateTime requested.
Warning 6: Normalized/laundered field name: 'geoidheight' to 'geoidheigh'
Warning 6: Normalized/laundered field name: 'ageofdgpsdata' to 'ageofdgpsd'

julia> imshow(D, proj=:eqc, region=(135,240,34,57), par=(MAP_FRAME_TYPE="plain",))

j_u · October 30, 2021, 11:14pm

I like it sharp and simple - no confusion. Thank you for taking the time, for providing solution and for all the explanations. I will look forward to GMT v0.38.1.

joa-quim · October 31, 2021, 12:01am

Done. You can now do ] up (but possibly you need to get out of the #master version)

j_u · October 31, 2021, 3:31pm

Thanks for the info.