Shape file to GMT

Q&A
1

Dear all,

I know the similar issue in another question..

But, my question is a little bit different.

I would like to map the fault as well, but I followed this document to open a shape file in the GMT.

What’s the problem in my GMT script?

I made my GMT script below:

#!/changwan/bin/env bash

region="-R127/130/34/38"
#projection
#input_path="/home/changwan/GMT/Fault_reference_line/Fault/YSF_reference line.shp"

gmt begin YFS png 

ogr2ogr -f OGR_GMT test_shape.gmt 'YSF_reference line'.shp

gmt convert test_shape.gmt -bo2f > test_shape.bf2

gmt select test_shape.bf2 -bi2f -bo2f $region > test_shape_final.bf2


         gmt coast $region -Gdarkseagreen2 -Scornflowerblue 
#        gmt plot test_shape_final.bf2 -bi2f -Wthick,red -Gred
#        gmt psxy 'YSF_reference line.shp' -bi2f -Wthick,red -Gred 
#        gmt plot test_shape.gmt -Wthick,red -Gred 
  
         gmt basemap --MAP_FRAME_TYPE=fancy -Ba2f1g1 -BWSne+t"Yangsan Fault System"

gmt end show

And I got this error:

ERROR 4: Unable to open YSF_reference line.shx or YSF_reference line.SHX. Set SHAPE_RESTORE_SHX config option to YES to restore or create it.
FAILURE:
Unable to open datasource `YSF_reference line.shp' with the following drivers.
  -> `PCIDSK'
  -> `netCDF'
  -> `PDS4'
  -> `JP2OpenJPEG'
  -> `JPEG2000'
  -> `PDF'
  -> `MBTiles'
  -> `EEDA'
  -> `ESRI Shapefile'
  -> `MapInfo File'
  -> `UK .NTF'
  -> `OGR_SDTS'
  -> `S57'
  -> `DGN'
  -> `OGR_VRT'
  -> `REC'
  -> `Memory'
  -> `BNA'
  -> `CSV'
  -> `NAS'
  -> `GML'
  -> `GPX'
  -> `LIBKML'
  -> `KML'
  -> `GeoJSON'
  -> `GeoJSONSeq'
  -> `ESRIJSON'
  -> `TopoJSON'
  -> `Interlis 1'
  -> `Interlis 2'
  -> `OGR_GMT'
  -> `GPKG'
  -> `SQLite'
  -> `OGR_DODS'
  -> `ODBC'
  -> `WAsP'
  -> `PGeo'
  -> `MSSQLSpatial'
  -> `OGR_OGDI'
  -> `PostgreSQL'
  -> `MySQL'
  -> `OpenFileGDB'
  -> `XPlane'
  -> `DXF'
  -> `CAD'
  -> `Geoconcept'
  -> `GeoRSS'
  -> `GPSTrackMaker'
  -> `VFK'
  -> `PGDUMP'
  -> `OSM'
  -> `GPSBabel'
  -> `SUA'
  -> `OpenAir'
  -> `OGR_PDS'
  -> `WFS'
  -> `WFS3'
  -> `HTF'
  -> `AeronavFAA'
  -> `Geomedia'
  -> `EDIGEO'
  -> `GFT'
  -> `SVG'
  -> `CouchDB'
  -> `Cloudant'
  -> `Idrisi'
  -> `ARCGEN'
  -> `SEGUKOOA'
  -> `SEGY'
  -> `XLS'
  -> `ODS'
  -> `XLSX'
  -> `ElasticSearch'
  -> `Walk'
  -> `Carto'
  -> `AmigoCloud'
  -> `SXF'
  -> `Selafin'
  -> `JML'
  -> `PLSCENES'
  -> `CSW'
  -> `VDV'
  -> `GMLAS'
  -> `MVT'
  -> `TIGER'
  -> `AVCBin'
  -> `AVCE00'
  -> `NGW'
  -> `HTTP'
plot [WARNING]: File test_shape_final.bf2 is empty!

like the previous question.

Then, I open the GMT file. test_shape.gmt. I can see this below:

gmt_file
gmt_file1634×1256 445 KB

I think ogr2ogr is working very well because we can see coordinate and other information.

I tried plot and psxy module to open shape file without using ogr2gor, but I failed, and only the cursor blinked without any map…

The last warning said the file is empty, and then the problem is gmt select module?
What’s the problem in my GMT script?

Regards,
Changwan

2

Yes, it looks like your ogr2ogr command is working. The error message either means that you don’t have the “.shx” file that goes with your “.shp” file or the programs are getting confused by the file name that contains a space. I see that you put quotes around the file name to help it read the file name, but the “ogr2ogr” program may still not work completely with that file name. I would try renaming all the files of your Shapefile (usually there is datafile.shp, datafile.shx, datafile.prj, and datafile.dbf) to a file name without a space.

The second problem is that it seems your fault is located at longitude 102-103 E, 48 N, which is outside of the region that you specified. That means that there are no points after you do the gmt select.

3

I appreciate that.
I will try with your suggestion and upload the result again!
Thanks a lot!

Changwan

4

Yeah! Your solution was successful!
My mistake in *.gmt file, the name of the file, and its address!
In the address, there should have been *.shp file and other files for ESRI shape file including *shx.

I got the new gmt file below:

gmt_file_YSF
gmt_file_YSF1425×1288 346 KB

And, I corrected the script line in the GMT plot module.

gmt plot test_shape_final.bf2 -bi2f -Wthick,red

Finally, I got a neat fault line!

I hope this topic would be helpful to others.
Thank you so much!
Changwan

5

I got a good result, but I have another error and warning:

ERROR 1: PROJ: proj_identify: Cannot find proj.db
Warning 1: The output driver does not seem to natively support Integer64 type for field id. Converting it to Real instead. -mapFieldType can be used to control field type conversion.

In the input address, there are *.shp, *.shx, and *.prj files.
Why did this error show? and How can I solve the warning?

6

The error means that the Proj library that GDAL is using for coordinate transformations does not find a file named “proj.db”. That file is a SQLite database where all the parameters needed for the transformations are stored. Environment variable “PROJ_LIB” should point to the directory where the proj.db file is. You have some issue with your GDAL or Proj installation. If you do not change the coordinate system then the error may not affect you.

Warning means that you have some field of type 64 bit integer in your shapefile but GMT format does not supports that datatype so it is changed into datatype “real”.
You can avoid the warning by defining the datatype mapping yourself as documented in https://gdal.org/programs/ogr2ogr.html
-mapFieldType Integer64 Real.