abduhbm / zmapio Goto Github PK
View Code? Open in Web Editor NEWPython library for reading and writing map gridded data using ZMAP Plus ASCII Grid format
License: MIT License
Python library for reading and writing map gridded data using ZMAP Plus ASCII Grid format
License: MIT License
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.
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!
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,
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.