GDAL librarify

[joa-quim]

Add structure to call GDAL programs via its C API.
Current implementation already allows calling the gdal_grid program using a temporary -A option of nearneighbor. An annoying issue is that, so far, simple input data must be passed in ascii via the GDAL's VRT mechanism.

As example, put the following into a lixo.csv file

2.53857 5.16657 0
2.48365 6.26811 1
8.68883 4.55983 2
4.11104 7.78623 3
1.79704 6.53027 4
7.17493 3.81713 5
3.41052 8.18161 6
8.35062 1.43348 7
8.1706 4.46765 8
5.27815 1.15172 9

create a lixo.vrt file with

<OGRVRTDataSource>
    <OGRVRTLayer name="lixo">
        <SrcDataSource>lixo.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

and run

gmt nearneighbor lixo.vrt -R0/10/0/10 -Gjunk.nc -I0.05 -A"-a nearest"

to get the same as the nat_nn.sc test (but, just realized, UD flipped)

[PaulWessel]

Since you are working in a WIP branch, do we really need to approve commits?

[joa-quim]

Hm, I didn't request a review. It was some github ghost.
But you could give a look at it and even try. Dave's people wanted a fast NN algo. This one is (and multi-threaded though not sure if by default)

[PaulWessel]

Ok, but later this week. We said we would release 6.1 soon (February) so we should discuss a checklist with @serisman and see if there are specific things we need to complete first.

[joa-quim]

@PaulWessel please have a go with this. I added a -Ns option to nearneighbor that does the same as -Nn but using GDAL.
Also added a wrapper module called gmtgdal which is a wrapper to the gmt_gdal_librarified() functions. The name is a bit poor and if ofc open to discussion. As well as its functionality, which for now is

usage: gmt gmtgdal <infile> -A<prog> [-F<gd opts>] [-G<outgrid>] [-M[+r+w]] [-W<outds>] [-R<xmin>/<xmax>/<ymin>/<ymax>[+r][+u<unit>]] [-V[<level>]]
           Specify input file name
        -A Specify the GDAL program name.
        -F list of GDAL options for the selected program in -A wrapped in double quotes.
        -G Name of output grid.
        -M Means to both read and write with GDAL. Use +r to only read or +w to write.
        -W Name of output dataset when writen by GDAL.

Only the GDAL's gdal_grid, gdalinfo and gdal_dem programs are usable (and latter very incomplete)

[PaulWessel]

OK, trying nearneighbor:

gmt nearneighbor lixo.vrt -R0/10/0/10 -Gjunk.nc -I0.05 -Ns
nearneighbor [ERROR]: Distance units must be one of d|m|s|e|f|k|M|n|u

What am I doing wrong?
Also, we cannot have people building vrt files - that would have to happen automatically if -Ns.

[joa-quim]

Ok, do

gmt nearneighbor lixo.cvs -R0/10/0/10 -Gjunk.nc -I0.05 -Ns

nearneighbor does not know how to read .vrt files.

If you were using gmtdal, then

gmt gmtgdal lixo.vrt -R0/10/0/10 -Gjunk.nc -I0.05 -F"-a nearest" -M+r

should work too.

[PaulWessel]

I switched to the GDALlibrify branch, rebuilt everything, and

gmt nearneighbor lixo.csv -R0/10/0/10 -Gjunk.nc -I0.05 -Ns
nearneighbor [ERROR]: Distance units must be one of d|m|s|e|f|k|M|n|u
pwessel@MacBeth:~/UH/RESEARCH/CVSPROJECTS/GMTdev/gmt-dev/dbuild-> gmt gmtgdal lixo.vrt -R0/10/0/10 -Gjunk.nc -I0.05 -F"-a nearest" -M+r

ERROR: No module named gmtgdal was found. This could mean one of four cases:

1. There actually is no such module; please check your spelling.
2. You used a modern mode module name while running in GMT classic mode.
3. Module exists in the GMT supplemental library, but the library could not be found.
4. Module exists in a GMT custom library, but none was specified via GMT_CUSTOM_LIBS.
   Shared libraries must be in standard system paths or set via environmental parameter DYLD_LIBRARY_PATH.

Do you have any other settings in cmake that needs to be applied?

[joa-quim]

Ah, forgot about it. I am using the set (EXTRA_MODULES gmtgdal.c) in ConfigUser.cmake. With it, any of

gmt gmtgdal lixo.vrt -R0/10/0/10 -Gjunk.nc -I0.05 -F"-a nearest" -M+r
and
gmt gmtgdal lixo.csv -R0/10/0/10 -Gjunk.nc -I0.05 -F"-a nearest"

work. But something screwed meanwhile with nearneighbor -Ns. Will have to debug and fix.

[joa-quim]

Wait, don't know why I thought it hadn't work. It does

$ gmt nearneighbor lixo.csv -R0/10/0/10 -Gjunk.nc -I0.05 -Ns -V
nearneighbor [INFORMATION]: Reading Data Table from File lixo.csv
nearneighbor [INFORMATION]: Cartesian input grid
nearneighbor [INFORMATION]: Cartesian input grid
nearneighbor [INFORMATION]: gdal options used: -ot Float32 -txe 0.000000 10.000000 -tye 0.000000 10.000000 -outsize 201 201 -a nearest:radius1=0.000000:radius2=0.000000:nodata=NaN -of MEM
nearneighbor [INFORMATION]: Writing grid to file junk.nc

[PaulWessel]

Since your gdal_vector calls GMT_Read_Data I think you need to add support for -i, -h, -q. All you need is to add them to the OPTIONS and the usage part and the rest is automatic. Plus man page.

[joa-quim]

It already has. So it misses -q #define THIS_MODULE_OPTIONS "->RVbhi" From: Paul Wessel <[email protected]> Sent: Sunday, February 9, 2020 10:57 PM To: GenericMappingTools/gmt <[email protected]> Cc: Joaquim Manuel Freire Luís <[email protected]>; Author <[email protected]> Subject: Re: [GenericMappingTools/gmt] WIP: GDAL librarify (#2583) Since your gdal_vector calls GMT_Read_Data I think you need to add support for -i, -h, -q. All you need is to add them to the OPTIONS and the usage part and the rest is automatic. Plus man page. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#2583?email_source=notifications&email_token=AAEDF2ODHR5PVSWFFY5SAFDRCCC4BA5CNFSM4KM4BVOKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELG2YJA#issuecomment-583904292>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAEDF2LGVR5I522FOYR662LRCCC4BANCNFSM4KM4BVOA>.

[joa-quim]

OK, the wrapper program name. I named it gmtgdal as a lower case version of GMTGDAL. It's true that gmt-proggys tend to use table data (but several do not), whilst grd-proggys deal with grids, but several deal with grids and do not have grd in the name.

So far gmtgdal deals almost exclusively with rasters (grids and images), so a grdgdal name would be more appropriated. On the other hand, the gmt prefix has own advantages. We can call them as

gmt gdal ...

and at a certain point I think it would be useful to make it write OGR formats out of our table data files (thus avoiding the need of system calls to ogr2ogr).

All of these to say. I have no strong opinion on a gmtgdal vs grdgdal vs someothername

Opinions?
(@seisman)

[joa-quim]

ping @PaulWessel @seisman

[PaulWessel]

A couple of thoughts: I think for the command line it does not matter. However, I know for the API it is better if a module only returns one type of thing. We are even considering to break some old modules into several to be simpler, and maybe we will do that with modules that returns different things depending on options. Going forward with new modules, I think we should try hard to only make them produce one thing. So I think grdgdal sounds like a much better name for now, and if you end up making a gmtgdal that deals with tables down the line then that just is another well-named module.

So I think grdgdal for the current module would be better.

[joa-quim]

Right. Agree with that. A gmtgdal might (or not) see the light of day further down the road.

[PaulWessel]

Looks good, but the example in the man page uses -M+W instead of -M+w.