How to image dataset in the text file with GMT

[C.SUN]

Hi, folks

I would like to open GPR radargram (b_scan_image) made by *.txt format with GMT.
I got this text file (30 by 301 table) from a GPR binary file with a certain format by using my own Fortran90 program that have signal processing algorithm. I selected GMT to image my final GPR data.

I used xyz2grd to change *.txt file into *.grd file, and then used grdimage to image *grd file.

My scripts is below:

gmt begin radargran.png
gmt xyz2grd B_SCAN_IMAGE_GPR.txt -R/0/30/0/301 -Gtext.grd -I1/1 -An
gmt grdimage text.grd -R0/30/0/301 -JX5i -Ccolors.cpt
gmt end show

I got this error message:grdimage [ERROR]:

Passing zmax <= zmin prevents automatic CPT generation!
grdimage [ERROR]: Failed to read CPT colors.cpt.
psconvert [ERROR]: Unable to decode BoundingBox file /home/changwan/.gmt/sessions/gmt_session.21647/psconvert_21657c.bb (maybe no non-white features were plotted?)

Ultimately, I want to make three dimensional diagram of GPR data like these:


(reference: https://www.bradyflinchum.com/)

As a first step, I would like to image b_scan_image.
I attached my text fileB_SCAN_IMAGE_GPR.txt (150.2 KB) .
Please help me out!

Changwan

[pwessel]

The error message has to do with your colors.cpt file which it says it failed to read, so everything falls apart after that (such as psconvert). How did you make the colors.cpt file and what does it look like?

[C.SUN]

Thanks. Professor Paul.

I made a simpler table (11 by 10)B_SCAN_IMAGE_GPR.txt (1.8 KB) , and extracted minimum and maximum values in the text file by using python code.

I used the both values to make the colors.cpt
I change my gmt script below:

gmt makecpt -Crainbow -T-52416.0/32640.0/1000.0 > colors.cpt

gmt begin b_scan png
       gmt xyz2grd B_SCAN_IMAGE_GPR.txt -R0/11/0/10 -Gtext.grd -ZTLa -I1/1
       gmt grdimage text.grd -R0/11/0/10 -JX5i -Ccolors.cpt
gmt end show

Nevertheless, I just got the errors:

After grdimage, the final result shows this :

I think I am using xyz2grd in a wrong way. In the table of text file, I know the number of numerical values–110. BUT, the first error says " found 33 records" Why?

Which option do I have to add in the xyz2grd module?

Regards,
Changwan

[pwessel]

I think you will need to use upper case A in -ZTLA. Read the documentation for the difference between a and A. Also, while I cannot see it I believe you need to add -r to your xyz2grd command to say these are pixel-registered and not gridline-registered grid.

[C.SUN]

I really appreciate that.
I got this image:

But, I expected this kind of image :

This is imaged by python. I need more works.
Anyway, I first got an image from GMT by reading the text file at this moment.
I will question again after reading a GMT documentation.

Thanks a lot
Changwan

[Joaquim]

I obtain similar to your python example (except ud-flipped) if I read your file like this in Julia

D = gmtread("B_SCAN_IMAGE_GPR.txt");
size(D[1]) = 
(301, 30)

imshow(D[1].data, figsize=2, aspect=1, shade=true, fmt=:png)

Which means that the data is stored by rows. First row has 30 columns, next row the next 30 columns …
You have to convert this info in your xyz2grd command.

[C.SUN]

Hi, Joaquim. Thanks for your comment. I should have given more information about my data in the first question.

The real data have the structure 4096 by 30 obtained from GPR exploration, but I cut it small for usefulness.

The uploaded data structure is 301 by 30. As you can see the image from python that I uploaded, we can find a unambiguous image than that from my GMT script.

The column, 30, is the ground track–the direction of the exploration. The row, 301, is the depth.
The interval of the column is the 0.05 m, and the interval of the depth is 0.005 m.

Thank you for your time and effort for my problem. Please give me constructive comments to update my GMT script.

[C.SUN]

Yeah. You are right.

I used xyz2grd in my script:

gmt makecpt -Crainbow -T-3373312.0/2784704.0/100.0 > colors.cpt

gmt begin b_scan png
     gmt xyz2grd B_SCAN_IMAGE_GPR.txt  -R0/300/0/29 -Gtext.grd -ZTLA -I1/1 
     gmt grdimage text.grd -R0/300/0/29 -B -JX5i -Ccolors.cpt
gmt end show

With this script, I got the image above.

[Joaquim]

Well, one of the reasons I wrote those Julia programs is to not have to deal with data reformating anymore

However, visiting the xyz2grd I don’t think it is the right tool for this. That program is for data with x,y,z or just z but what you have is already a matrix. The simplest is probably to stick it a surfer6 ascii header or an ESRI bil header, and read the file directly in GMT. You might still face data transposition issues but grdmath can deal with it if I remember well.

EDIT: OK, so I was wrong (again). Good

[C.SUN]

OK. Great. I will try other modules in the GMT documentation and then ask again in this forum!!

Thanks a lot
Changwan