Plot shapefiles in X-Y coordinates

Q&A
1

Hi there,

I’m not used to work with shapefile format and I have a bunch I’d like to process for a geologic map.

There are actually 8 file formats associated to the data :

  • .cpg
  • .dbf
  • .lyr
  • .prj
  • .sbn
  • .sbx
  • .shp
  • .shx

I think only shp and shx are useful for GMT, the others being for ArcGIS, I suspect. Again, I’m not sure… I don’t even know how to extract information properly from those.

gmt info file_2154.shp
    file_2154.shp: N = 24303	<84725.096698/1201589.26104>	<6032821.09236/7106074.34002>

As you can see, the units are in meter (UTM). From the documentation I know that the projection is Lambert II extended (-Jl) centred at {3ºE ; 46º30’N} (or 700000/6600000) , using Clarke-1880-IGN ellipsoïd.

In order to have a better grip on GMT, I’d like to convert the file to full WGS84 lon/lat system. When I explore the mapproject module, I got a bunch of errors, either because of the units (-R[...]+u1) or because of the projection

area="-R-22841.3905/1712212.6192/1087335.9605/2703971.1254+u1"
projection="-Jl700000/6600000/44/49/1:50000"
gmt mapproject $area $projection $data -C -E -I -F --PROJ_ELLIPSOID=Clarke-1880-IGN > test.txt
    Option -R: The +u<unit> modifier only applies to projected coordinates. Your unit 1 is ignored

Ultimately, I’d like to plot all the files, associated colours to strata and everything. But for that, I need to understand how to deal with them.

Any help ?

I can provide some files if needed (shp ~20Mb , shx ~20Kb)

Guillaume

1 Like
2

Did you try +ue for meter? +u1 is not a recognized thing.

3

When I need to plot shapefile I use ogr2ogr to convert the shape to a GMT multiple segment file

for example

ogr2ogr -f "GMT" $gmt_format_file $shape_file
gmt mapproject $gmt_format_file ...
gmt psxy $gmt_format_file ...

ogr2ogr can be found with gdal library

Hope it can help you

4

I tried, it doesn’t solve the issue. My text file has a first column filed with 80-ish values and a second column of zeros.

I think GMT does that under the hood.
I’m currently able to plot the shape file “as is”, but I don’t know how to handle the informations / values and so on, so I can’t associate a colour to specific entries

5

  • All of those .xxx are part of the shapefile (no ArcGis involved) though the .cpg and .lyr are new to me.

  • mapproject requires -R in geogs, even when doing inverse projections. That’s a quite annoying thing, I agree.

  • doing it with Julia should be pretty straightforward.

6

You can do the projection change with GDAL also. It will get the input projection information from the shapefile data automatically. Look at the options for “ogr2ogr” that can both project and convert to GMT vectors.

7

Check .shp metadata using ogrinfo

ogrinfo -so file_2154.shp

it’ll show layer names. then run the ogrinfo with the layer name appended to get that layer metadata.
-so means “summary only”, without it ogrinfo will dump all the file/layer content on your terminal.

8

Hehehe. Always.

area="-R10/20/33/54+ue"
projection="-Jl3/46.5/44/49/1:50000"
gmt mapproject $area $projection $data -C -E -I -F  > test.txt

test.txt stays with 1 column of 80-ish values and 1 column of zeros.

Running this command (debug) :

gmt plot -W -Bafg file_2154.shp -png test -Vd

I got this :

gmtinfo [DEBUG]: Running ogr2ogr -mapFieldType Integer64=Integer -skipfailures -f "OGR_GMT" "gmt_ogr_5387.gmt" "file_2154.shp"

Revised options: -W -Bafg file_2154.shp -Vd -R40000/1280000/6000000/7160000 -JX15c/15c

Plotting segment 0
...
Plotting segment 80072

(if I plot all the files, the figure is quite nice even without colouring. This is “just” one file :

test
test1959×1835 730 KB

All it gives me is :

INFO: Open of `file_2154.shp'
      using driver `ESRI Shapefile' successful.
1: file_2154 (Line String)

Nothing else

9

then run the ogrinfo with the layer name appended to get that layer metadata:
ogrinfo -so file_2154.shp file_2154

“Line string” means contour lines if i remember correctly

10

Ok, that gave the information I can read from the .prj file (projection, ellipsoid …)

If I remove the “summary”, I got the polygon coordinates for each element (>80k). Now I need to find what these elements are supposed to be.

11

There should be something else right under the layer projection info, a list of contour feature parameters (set at some values for every contour).

Can be easier to run ogr2ogr exactly as gmt does under the hood and study the header of the output file and few first records.
ogr2ogr -mapFieldType Integer64=Integer -skipfailures -f "OGR_GMT" "gmt_ogr_5387.gmt" "file_2154.shp"

12

Yep, I got that :

# @NCODE|DESCR|TERRE_MER|CODE_LEG|NOM_SYMB|WT_SYMB|C_SYMB|M_SYMB|J_SYMB|N_SYMB
# @Tinteger|string|string|integer|string|double|double|double|double|double

Where CODE_LEG | NOM_SYMB | … refer to the spatial database somehow

13

Now you actually need to check documentation hopefully accompanying your shapefile data. Otherwise it is going to be very hard to guess the meaning of all these parameters:
CODE|DESCR|TERRE_MER|CODE_LEG|NOM_SYMB|WT_SYMB|C_SYMB|M_SYMB|J_SYMB|N_SYMB

Well, DESCR must be some sort of description string. Others are not clear. Every contour header should have all these parameters set to some values (according to the value types). Once you know what exactly you want, you can extract it with something like
ogr2ogr -where CODE_LEG=some_value -f GMT "gmt_ogr_5387.gmt" "file_2154.shp"

14

Thanks @mkononets :slight_smile:
The doc talk about the structure of the database, but not so much about the parameters themselves. Hopefully the “description” will help !

15

Ok ! So another file contains actual polygons (surfaces) with explicit names. There’s a code associated starting with @ (e.g. @Dq3).
Is it possible to list all those then use [..] -where CODE=@Dq3 [...] to split the file into small bricks ?

16

I think you need to exclude @D, which indicates the beginning of the parameter list
try -where CODE=123 if CODE is a number (type real or integer) or -where CODE='123' if CODE is a string type

17

Hm, I tried several ways … it oddly doesn’t work.

I switched to an easier test. Using ogrinfo -so I’ve isolated this one :

TERRE_MER: String (1.0)

Which has “T” or “M” as value

ogr2ogr -where TERRE_MER='T' -f GMT "file_2154.gmt" "file_2154.shp"
    ERROR 1: "T" not recognised as an available field.
    ERROR 1: SetAttributeFilter(TERRE_MER=T) on layer 'file_2154' failed.

'T', T or "T" doesn’t change anything

18

try putting the expression after -where in double quotes:
ogr2ogr -where "TERRE_MER='T'" -f GMT "file_2154.gmt" "file_2154.shp"

19

Alright ! That throw a bunch of warning due to Integer conversion. But it works.

Thanks again @mkononets. Now I have to define what I want and how to convert it :slight_smile:

20

gmt uses -mapFieldType Integer64=Integer to avoid Integer64 to Real conversion warnings
there’s also -skipfailures but this might not let you see if something went radically wrong.

BTW you do not necessarily have to convert to GMT format. If you leave it out and give some other .shp output filename, the result will also be in .shp format (with no potentially incompatible type conversions). However .gmt files are human-readable, which is helpful sometimes.

1 Like