Gmtselect using -e pattern

Q&A
1

Dear masters of the GMT multiverse,

I’ve been counting points within a set of polygons with gmtselect in this way:

grd2xyz mync.nc | gmtselect -Fpolygons.shp | wc -l

I am now wondering if it is possible to modify the above call to choose a specific polygon with gmtselect. I went through the gmtselect man page and did not get how to select a polygon based on its name (“ChooseMe” in the 6th column called “names”), so, using the -e option. I tried the following without success:

grd2xyz mync.nc | gmtselect -Fpolygons.shp -e"ChooseMe" | wc -l

or

grd2xyz mync.nc | gmtselect -Fpolygons.shp -a6=names -e"ChooseMe" | wc -l

Not sure if the -e option has been implemented to do this or I am doing something wrong here. Any support is much appreciated.

All the best,

Gery

2

Mmm, maybe you can use grdmask that allow to use -a and -e (instead of grdxyz and gmt select).

Maybe this could work:

gmt grdmask polygons.shp -Rmync.nc  -a6=names -e"ChooseMe" -GMask.nc -NNAN/1/1
gmt grdmath Mask.nc mync.nc = new_nc.nc 
gmt grd2xyz new_nc.nc | wc -l

Does your grid (mync.nc) has holes?

3

Thanks Esteban for the idea, unfortunately it did not work, grdmask says:

grdmask [ERROR]: Without -S, we expect to read polygons but none found

it seems to me that the “ChooseMe” polygon is not being identified by the -e option, so it did not select that polygon. I have the impression that I am misusing or misunderstanding the purpose of the -e flag (or its combination with the -a flag). I thought the problem was the shapefile itself, but this gmtseslect call works as expected (ie. all points falling within/intersecting the whole polygon are selected):

 grd2xyz mync.nc | gmtselect -Fpolygons.shp

Btw, I am using v6.4.0 :slightly_smiling_face:

4

Could you share you shape file (or a sample of it)?

5

as far as I can understand, -a only allows numeric values, not string literals like “ChooseMe”
see examples by @Esteban82 below

6

I have used -e with string. But it has to be part of the record (not the header)

For example:

#Long Lat  String
-22      -44    ChooseMe
-21      -44    ChooseMe
-20      -44    ChooseMe
7

I think it is possible to extract a specific polygon from a multipolygon shapefile using gdal’s ogr2ogr, something like:

ogr2ogr -f ChooseMe.shp -dialect sqlite -sql "SELECT * from polygons WHERE name='ChooseMe'" polygons.shp

a shorter version using -where should also work:

ogr2ogr -f ChooseMe.shp -where "name='ChooseMe'" polygons.shp
8

-a is for aspatial parameters, that is, those from the headers.

9

Many thanks for all the answers, I just found that @Esteban82 and @mkononets were right, “ChooseMe” is in the header because a shp-to-gmt conversion with ogr2ogr showed me that. So, the -e flag does not apply here as pointed by @Esteban82. Not sure if a specific polygon can be selected with the -a flag as I expected, though.

10

How to use -a

Check this test file (capitals.gmt)

It has this format (it seems to me that it was a shp file converted to gmt format).

# @VGMT1.0 @GMULTIPOINT
# @R1.35/359.83/-51.72/64.15                     
# @Jp"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "
# @Jw"GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]"
# @Npopulation|name
# @Tinteger|string
# FEATURE_DATA
>
# @D1424400|"KABUL"
69.17	34.50
>
# @D244153|"TIRANA"
19.82	41.33
>
# @D1740461|"ALGER"
3.00	36.75
>
# @D3519|"PAGO PAGO"
189.28	-14.31
>

Use -a

You can use -a to add a third column to the records base on the info of the header.

For example, to add a second column (0 and 1 are for X and Y) with the populutation data you can use this command:

gmt select @capitals.gmt -a2=population
>
69.17   34.5    1424400
>
19.82   41.33   244153
>
3       36.75   1740461
>
189.28  -14.31  3519
>
1.5     42.5    16151

You could also add the name with:

gmt select @capitals.gmt -a2=population,3=name
>
69.17   34.5    1424400 KABUL
>
19.82   41.33   244153  TIRANA
>
3       36.75   1740461 ALGER
>
189.28  -14.31  3519    PAGO PAGO
>
1.5     42.5    16151   ANDORRA LA VELLA

Filter data

Now, you could use -e to filter your records

gmt select @capitals.gmt -a2=population,3=name | gmt select -e"KABUL" 
>
69.17   34.5    1424400 KABUL

BUG??

I try to do this but I got an empty answer. I think it is a bug.
gmt select @capitals.gmt -a2=population,3=name -e"KABUL"

2 Likes
11

cool, it works! I had problems with -a at some point before. Just as you wrote above, empty output or nothing plotted when I called gmt text and tried to populate text strings on the fly using -a. Now seems working totally fine.

12

Thanks @Esteban82 “muchas gracias genio”, that one did the trick. In case someone else is trying something similar, I got it working first with the good ogr2ogr like this:

ogr2ogr -f "ESRI Shapefile" selected.shp poi.shp -dialect sqlite -sql "select point.geometry from poi point, 'pol.shp'.pol polygon where st_intersects(poi.geometry, pol.geometry) and pol.names like 'ChooseMe'"

but GMT is like 1-min faster in a 12th Gen Intel© Core™ i7-1255U × 10 16 RAM.

Hope this helps someone else out there.