Thank you for this repo - it is really useful.

I have come across an issue with the creation of the xyz points from the zmap file. When creating the xy arrays the points are incorrectly located. I would expect to see the output points in the cell centers. From inspecting the code I think this occurs when inputting the cell edge min_x, max_x etc. values into np.linspace rather than the cell center min and max values. Hopefully the example below illustrates the point and provides a fix.

Here is an example code that reads a zmap file, creates a dataframe of points using .to_pandas() and creating a shapefile from that to view the point locations

from zmapio import ZMAPGrid
import geopandas as gpd
import numpy as np
from shapely.geometry import Point

z_file = ZMAPGrid(os.path.join(environ['USERPROFILE'],'Desktop','UK_DEM.dat')) # read file

df=z_file.to_pandas() # make dataframe
df.dropna(inplace=True) # drop nulls
df['geometry']=df.apply(lambda x: Point(x['X'],x['Y']),axis=1) # make point geoms
df=gpd.GeoDataFrame(df,crs='epsg:4326') # make geodataframe and export as shapefile
df.to_file(os.path.join(environ['USERPROFILE'],'Desktop','UK_DEM_POINTS.shp')) # export shapefile

When plotting the shapefile over the zmap in ArcGIS/QGIS we can see that the points are located in the cell top left in the top left corner:

Near the center in the middle of the zmap:

Point in the bottom right of cell toward the bottom right corner of the zmap:

To make these consistently cell center when creating the x and y meshgrids the code would need to look something like this:

# get cell size
dx=(max_x-min_x)/no_cols
dy=(max_y-min_y)/no_rows
# make arrays of x and y
x = np.linspace(min_x+(dx/2), max_x-(dx/2), no_cols) # shift min value right by 1/2 cell and max value left by 1/2 cell
y = np.linspace(max_y+(dy/2), min_y-(dy/2), no_rows)
# meshgrid
X,Y=np.meshgrid(x,y)
# export dataframe to view points location compared with zmap
center_points=gpd.GeoDataFrame({'X':X.ravel(),'Y':Y.ravel(),'geometry':map(Point,zip(X.ravel(),Y.ravel()))},crs='epsg:4326')
center_points.to_file(os.path.join(environ['USERPROFILE'],'Desktop','UK_DEM_POINTS_CENTER.shp'))

The resultant shapefile of points plotted over the zmap shows the points now placed in cell center. Again the image below is the top left of the zmap.

I hope that is clear. If not let me know.

Many Thanks,

Here's the minimal reproducable example

from zmapio import ZMAPGrid

zmap = ZMAPGrid('D:/tmp/test_pet_sq.txt')

out = ZMAPGrid(z_values = zmap.z_values, 
                min_x = zmap.min_x,
                max_x = zmap.max_x,
                min_y = zmap.min_y,
                max_y = zmap.max_y
                )

out.comments = ['test', 'zmapio']
out.nodes_per_line = 3
out.field_width = 15
out.decimal_places = 4
out.name = 'zmapio_test'

out.write('D:/tmp/test_zmapio_sq.txt')

Files for testing and output generated by this code are attached. As you can see the actual data in an output file differs from an input file.

If I change out generation by the following:

out = ZMAPGrid(z_values = zmap.z_values.swapaxes(0, 1), 

generated file looks correct.

I'm not sure if it's input or output related issue but it doesn't look like a correct behaviour.

test_pet_sq.txt
test_zmapio_sq.txt

It's a very simple issue, RMS adds the following string

+Grid data starts after this line

after the end of header section (so basically after closing @). Zmapio tries to interpret this line as a z value and fails. You need to skip this line if occurs.

Great library for reading and writing zmap files. It would be great to add CRS transformations to this library. To do it correctly, every XY location needs transformation, not just the xMin, xMax, yMin, yMax. I got this far and it works using a GeoDataFrame for the transformation. However, I cannot represent the coordinates in the z_file since it only contains the min/max and coordinates are re-created by x = np.linspace(self.min_x, self.max_x, self.no_cols). The results are a slightly shifted/rotated grid. Perhaps it will work if I replace linespace by a loop of transformed values, but then that would not work when importing into other software.

If I use a CSV where every X,Y,Z is represented the transformation is perfect.

Maybe I'm doing something wrong, suggestions welcome.

The goal is to grid data in lat/long coordinates using griddata() and export to a zMap file in another specified CRS.

def writeZmap(xi, yi, zim, exportCRS):
success = False
try:
z = ma.filled(zim, np.nan)
x, y = np.meshgrid(xi, yi)
dat = np.column_stack([x.flatten(), y.flatten(), z.flatten(),])
df = pd.DataFrame(dat, columns=["X", "Y", "Z"])
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.X, df.Y, df.Z))
gdf = gdf.set_crs('EPSG:4326') # World Geodetic System
gdf = gdf.to_crs(exportCRS)
x_tx = np.unique(gdf.geometry.x)
y_tx = np.unique(gdf.geometry.y).T
z_tx = df['Z'].values.reshape(z.shape[0], z.shape[1]).swapaxes(0, 1)
z_file = ZMAPGrid(z_values=z_tx, min_x=x_tx.min(), max_x=x_tx.max(),
min_y=y_tx.max(), max_y=y_tx.min())
success = True
print('Zmap file created....')
return success, z_file
except Exception as e:
success = False
print('Error creating Zmap file: ' + str(e))
return success, None

Thank you for this library; I'm using it in gio.

The default null value of 1e+30 is passed to a format string like "{0:>{1}.{2}f}".format(1e30, 12, 4) (see line in code) and produces '1000000000000000019884624838656.0000' instead of '1e+30 or possibly '1.0000e+30'.

Could fix by using the str() of the null value for NaNs, and the current format string for 'normal' values. Or you could change the format string to something like "{0:>{1}.{2}e}" (which would use scientific notation for everything, not sure if that's desirable).

I will submit a PR for a fix; you can decide what to do with it!