Psxy file to grd?

Hi all!

I have a typical psxy input file with the sample points (contours of the bathymetry for Lake Baikal). They are in UTM 48 but I need to transform this into a grid to replace the bathymetry for some gravity reductions.

Here is a sample of the file. Any idea how to achieve this without many hiccups?

>
# @DPolyline|GSLAYER|-1450.00000|0.00000|3|25
657122.762641825 5847026.93765868
657168.335094928 5847083.85795265
657215.085401375 5847139.81487531
657262.99320167 5847194.78405775
657312.037632235 5847248.7415612
657362.197334495 5847301.66388745
657413.450464183 5847353.52798915
657465.774700847 5847404.31127975
657519.147257577 5847453.99164343
657573.544890923 5847502.54744466
657628.943911022 5847549.95753768
657685.320191911 5847596.20127563
657742.649182038 5847641.25851965
657800.905914948 5847685.10964755
657860.065020161 5847727.73556241
657920.100734219 5847769.11770088
657980.986911906 5847809.23804128
658042.697037634 5847848.07911142
658105.204236988 5847885.62399624
658168.481288436 5847921.85634515
658232.500635175 5847956.76037916
658297.23439714 5847990.32089777
658362.654383141 5848022.52328555
658428.732103142 5848053.35351855
658495.438780667 5848082.79817037
658562.745365331 5848110.84441803
658630.622545495 5848137.48004756
658699.040761029 5848162.69345927
659904.352562425 5848440.67732444
>
# @DPolyline|GSLAYER|-1350.00000|0.00000|252|24
648909.413787187 5841537.2303288
649078.642558269 5841522.01745956
>
# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
666774.617023622 5878993.91187628
667346.012602095 5879079.1589866
667971.864322102 5879112.09855081
668597.71604211 5879112.09855081
669245.527471591 5878980.34029397
669684.758866205 5878531.81821512
>
# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
666774.617023622 5878993.91187628
666664.800231739 5878983.12294987

Well, no escape. You will have to grid them surface datafile -R... -I... -G...
Note, wrap the numbers under triple back ticks instead of using > to quote.

Thanks for the tip. It is just that the file looks like this (without the |) and the Z values is in the header:

|# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
|669919.124812163 5876184.29223821
|670004.014773802 5876796.34931265
|>
|# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
|670004.014773802 5876796.34931265
|670009.763995602 5877451.90931406
|669911.600581159 5878060.54038019
|669684.758866205 5878531.81821512
|>
|# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
|666774.617023622 5878993.91187628
|667346.012602095 5879079.1589866
|667971.864322102 5879112.09855081
|668597.71604211 5879112.09855081
|669245.527471591 5878980.34029397
|669684.758866205 5878531.81821512
|>
|# @DPolyline|GSLAYER|-1650.00000|0.00000|212|27
|666774.617023622 5878993.91187628
|666664.800231739 5878983.12294987
|>

Normally you’ll have to use -a2=NAME but I do not see a @N record in you file that defines the NAME of the Z value, so out of luck?

This files comes from transforming a shape file through ogr2ogr so… maybe it is kind of broken. I guess I will have to make a weird script to fix it.

Try reading the original shapefile with GMT.jl

D = ogr2ogr("shape...")

and see if D[1].attrib is not empty.

Not empty at all.
How can I access the table and the value saved in “ELEVATION”?

So there you have your z. Loop over the D[k] elements and use that to create a 3rth column with the z’s
D[1].attrib["ELEVATION"]

Edit: … and don’t tell anything to @Andreas.

2 Likes

It works great!. I will leave this script here in case someone has a similar problem:

using GMT
D = ogr2ogr("bathymetry.shp")
 
DATA=zeros(32743,3)
c=1
for i=1:length(D)
     z = parse(Float64,D[i].attrib["ELEVATION"])
     for j=1:size(D[i],1)
         DATA[c,1]=D[i].data[j,1]
         DATA[c,2]=D[i].data[j,2]
         DATA[c,3]=z
         c=c+1;
     end
end
3 Likes

ogr2ogr has an option -zfield <field_name>

Uses the specified field to fill the Z coordinate of geometries.

(https://gdal.org/programs/ogr2ogr.html#cmdoption-ogr2ogr-zfield)

1 Like

Thanks @mkononets

@marianoarnaiz you can try to see if it works

D = ogr2ogr("bathymetry.shp", ["-zfield","ELEVATION"])