Hi Gery for some weird reason your gmt is conecting polygons but your code works great just have look I have GMT 6 as well
This issue has been fixed in GMT master (reading the .shp in psxy) but you would need to build GMT yourself.
Iām sorry I missed this question and thread back in July! The GMT Julia system sounds quite intriguing.
I recently learned a few new .shp shapefile tricks and wanted to share it with the group, in case anyone needs to know this esoteric stuff!
Joaquimās last post (7/21) mentions a fix in GMT to use a shapefile directly within a psxy command. That sounds very cool, but I havenāt tried it.
I typically use GDALās ogr2ogr command, much like what Gery offered (above).
Lately, I needed to access specific polygons from within a huge global shapefile, 2GB in size. Originally, I had created a huge .gmt counterpart (>4GB in size) with ogr2ogr, then, was reading through the entire huge file to map features for much smaller regions.
Recently, I discovered the -fid option in the ogr2ogr command which permits a very quick way to dump just one specific polygon (I typically use the -f OGR_GMT option to dump as a text file). It reduced run times that originally took several minutes down to a few seconds, simply by dynamically dumping polygons as needed, rather than reading through the entire dumped shapefile.
In my case, the next step might be to use ogr2ogr to create temporary .shp files for needed polygons, then feed those into psxy directly. Joaquim, if you read this, does that sound like a possible way to go? Might speed things up even more to by-pass the OGR_GMT dump-step.
And Iām intrigued with your intrigue 
I donāt remember what was the issue that was fixed in master because using shape files directly (with under the hood conversion with ogr2ogr) has been added quite some time ago.
I didnāt know about the -fid either. How do you use it? I mean, how do you know what is the feature that you want to select? But about your question, if you know what to select and use that knowledge to extract a small shape out of an huge one, than that seems useful. If easy to use we could envision a way to use that from GMT, though Iāve been saying that for a long time, the future is to read the ogr vector files directly. Like one can do from Julia and even restrict the extracted feature by reagion (i.e. using -R).
Thanks for the reply, Juaquim! Yes, Iād be all for some manner to be able to set -R and grab restricted polygons with a -R specification! I built something like that myself. Permit me to describe my situation.
The huge global shapefile I was referring to is a shapefile of water bodies (lakes, reservoirs, rivers, etc.). Based on that shapefile, I used gmtspatial to compute areas and centroids for each water body and dump them out. The dump also has a sequential water body index number. So, this data set acts like my āRosetta Stoneā for determining what is needed when a map is requested.
Once a specific region is defined for mapping (for example, a 3 degree square region, something like -R-97/-94/34/37), then the region is extended by a half degree, to make sure that most lake centroids, near the edges, are selected, too. A gmtselect command is issued, with the extended -R option, to capture all the centroids falling within the 4 degree square region, while retaining the water body index numbers. Based on these water body indices, the polygons IDs are now known.
The -fid option is useful if you want to ātrapā just one polygon at a time. At first, I simply created a do-loop to get each water body polygon inside the extended bounds, then map it. That works well and is pretty fast, but, as it turns out, āfidā is a special field recognized by OGR SQL, which means you can construct a sub-shape file based on a list of polygons.
Say you want to map polygons 5-15. The loop method goes like this:
i=5
while [ $i -le 15 ]; do
ogr2ogr -fid $i t.shp huge_shapefile.shp
psxy t.shp -J -R ā¦
i=$(($i+1))
done
Or, you can create a comma delimited list, and use one ogr2ogr command (very slick):
pshps=ā5,6,7,8,9,10,11,12,13,14,15ā
ogr2ogr -where āfid in ($pshps)ā t.shp huge_shapefile.shp
psxy t.sph -J -R ā¦
It would really be useful to be able to provide a -R specification and have a very fast extraction of required polygon segments, right from within a GMT command. Iām all for that! Until then, there are clunky ways to do it!
I donāt think weāll have a -R with the current ogr2ogr syscall. But you can do it with the Julia wrapper. Example (see how all polygons are extracted as long as they partially cross the BoundingBox)
using GMT
Pt = gmtread("C:\\programs\\compa_libs\\covid19pt\\extra\\mapas\\concelhos\\concelhos.shp", region=(-10,-6,40, 41));
plot(Pt, fmt=:png, show=1)
at this point you can proceed with the analysis in Julia or save it into a multi-segment file and continue with shell GMT
gmtwrite("sub_shape.dat", Pt)
