How to extract attribute values from shapefile for cpt file generation

Q&A
1

Aloha.
I’m trying to figure out how to extract attributed data from a shapefile so I can construct a cpt file (makecpt). I have successfully plotted a map from the shapefile using:

gmt psxy $OUTPUT_GMT $REGION $PROJECTION -W0.05+cf -aZ=KAcreFt -Crainbow.cpt -K > $OUTPUT_PS

but cannot figure out how to extract the KAcreFt attribute into a file that I can use to generate a colormap. I’ve tried gmt convert but seem to end up with all of the vertices of all of the polygons with the attribute value as the 3-rd column. Any advice is welcome.

Mahalo.
J.

2

May be you can refer to Coloring OGR_GMT polygon files based on attribute column. The key point is -G+z that is provided by GMT 6.x

3

Aloha.
I see I had this problem before and I don’t think it was ever solved. This answer is not helpful as I’ve tried a wide-range of makecpt and psxy arguments.

The problem seems to be that makecpt is using one of the geometry fields to generate the cpt. The attribute arguments of psxy don’t seem to have any effect and there is a cryptic comment about changing the selection of input variables in makecpt to somehow use i ; whatever that means.

Any advice is welcome.
J.

4

More experimentation shows that gmtinfo only sees the geometry of the shapefile and none of the attributes. Is there any method to read shapefiles that preserves attributes? psxy seems to be the only one.

5

gdal’s ogrinfo can extract all attribute values that exist in a shapefile. Assuming your attribute name is KAcreFt and shapefile name is shapefile_name.shp:

ogrinfo -sql "SELECT DISTINCT KAcreFt FROM shapefile_name" shapefile_name.shp
6

I know the answer. Ah but I don’t know Julia, but I have to show it anyway. Reading shapes and their attributes with GMT.jl is trivial

julia> D = gmtread("country_borders.shp")
        This file has islands (holes in polygons).
        Use `gmtread(..., no_islands=true)` to ignore them.
Vector{GMTdataset} with 4287 segments
Showing first segment. To see other segments just type its element number. e.g. D[2]

Attribute table
┌──────┬───────────────────────────┬────────────────────┬─────────────────────────┬──────────────────────────────┬───────────
│  Row │ NAME_SV                   │ NAME_ZH            │ NAME_ES                 │ SUBUNIT                      │ NAME_UR  ⋯
├──────┼───────────────────────────┼────────────────────┼─────────────────────────┼──────────────────────────────┼───────────
│    1 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    2 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    3 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    4 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    5 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    6 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    7 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    8 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│    9 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│   10 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│   11 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│   12 │ Indonesien                │ 印度尼西亚         │ Indonesia               │ Indonesia                    │ انڈونیشی ⋯
│  ⋮   │             ⋮             │         ⋮          │            ⋮            │              ⋮               │          ⋱
│ 4277 │ Wake                      │ 威克岛             │ Isla Wake               │ Spratly Islands              │ جزیرہ وی ⋯
│ 4278 │ Wake                      │ 威克岛             │ Isla Wake               │ Spratly Islands              │ جزیرہ وی ⋯
│ 4279 │ Wake                      │ 威克岛             │ Isla Wake               │ Spratly Islands              │ جزیرہ وی ⋯
│ 4280 │ Wake                      │ 威克岛             │ Isla Wake               │ Spratly Islands              │ جزیرہ وی ⋯
│ 4281 │ Clippertonön              │ 克利珀顿岛         │ Isla Clipperton         │ Clipperton Island            │ جزیرہ کل ⋯
│ 4282 │ Macao                     │ 澳门               │ Macao                   │ Macao S.A.R                  │ مکاؤ     ⋯
│ 4283 │ Macao                     │ 澳门               │ Macao                   │ Macao S.A.R                  │ مکاؤ     ⋯
│ 4284 │ Ashmore- och Cartieröarna │ 阿什莫尔和卡捷群岛 │ Islas Ashmore y Cartier │ Ashmore and Cartier Islands  │ جزائر ای ⋯
│ 4285 │ Bajo Nuevo                │ 巴霍努埃沃浅滩     │ Bajo Nuevo              │ Bajo Nuevo Bank (Petrel Is.) │ باخو نوئ ⋯
│ 4286 │ Serranilla Bank           │ 塞拉尼拉浅滩       │ Isla Serranilla         │ Serranilla Bank              │ سیرانیلا ⋯
│ 4287 │ Scarboroughrevet          │ 黄岩岛             │ Bajo de Masinloc        │ Scarborough Reef             │ سکیربورو ⋯
└──────┴───────────────────────────┴────────────────────┴─────────────────────────┴──────────────────────────────┴───────────
                                                                                            165 columns and 4264 rows omitted
BoundingBox: [117.70360790395524, 117.93482506600003, 4.03156159100007, 4.163414542001791]
Global BoundingBox: [-179.99999999999991, 180.0, -89.99999999999994, 83.63410065300008, 0.0, 0.0]
PROJ: +proj=longlat +datum=WGS84 +no_defs
WKT: GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

17×2 GMTdataset{Float64, 2}
 Row │     Lon      Lat
─────┼──────────────────
   1 │ 117.704  4.16341
   2 │ 117.704  4.16341
   3 │ 117.738  4.15724
7

Mahalo. The ogrinfo solution works although the syntax is confusing (layers, filenames, …). It took a while to get it to work since SQL syntax error are amazingly terse.

I also found this approach to work well. That is, writing a CSV file and then extracting whatever field is wanted. In the example below, there is only one field but gawk could be used to pick from many.

ogr2ogr -f CSV ./output.csv shapefile.shp
gmt6 makecpt cpt_input.csv -E50 -Ccyclic > ./kgal_cyclic.cpt

I find this to be much more flexible.
J

8

Second line should be:

gmt6 makecpt output.csv -E50 -Ccyclic > ./kgal_cyclic.cpt

9

And I will take a look at julia.