CRS transformation?

[baylejs]

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

[baylejs]

I got it working. Basically have to use scipy.griddata() to re-grid input on a regularly spaced output grid for zmap file format. In the following code, xi, yi are in lat/long coordinates, zim is 2D mask array of values at those locations. geopandas is used for CRS transformation.

def writeZmap(xi, yi, zim, exportCRS, gridResolution):
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)

	# create regular spaced grid in new CRS
	xmin = gdf.geometry.x.min()
	ymin = gdf.geometry.y.min()

	xmax = gdf.geometry.x.max()
	ymax = gdf.geometry.y.max()

	xii = np.linspace(xmin, xmax, num=xi.shape[0])
	yii = np.linspace(ymin, ymax, num=yi.shape[0])

	# meshgrid of transformed locations
	xx = gdf.geometry.x.values.reshape(yi.shape[0], xi.shape[0])
	yy = gdf.geometry.y.values.reshape(yi.shape[0], xi.shape[0])

	# sample values from CRS1 input array to CRS2 regular spaced zmap grid
	zi = griddata((xx.flatten(), yy.flatten()), z.flatten(), (xii[None, :], yii[:, None]), method='nearest')

	# rotate z array for z_file export
	zi = np.rot90(zi, k=3)

	z_file = ZMAPGrid(z_values=zi, min_x=min(xii), max_x=max(xii),
	                                    min_y=min(yii), max_y=max(yii))

	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