Revise SpatialReference properties, add geotransform #367
Conversation
add missing autotests for get_extent
- remove unused set_origin (assumed to not be needed?) - remove theta property, which was used in equations that assumed CW rotations (we use CCW elsewhere) - reverse sign in xll and xul to use CCW rotation convention - add generic cos_sin(deg), which snaps to right angles - use shorthand logic "x or y" instead of "x if x is not None else y"
Use this for export_array with .tif, since it is more robust with older versions of affine (before v.2, rotation was clock-wise, and after it is counter clock-wise. E.g. the added autotest previously failed with affine 1.2.0 packaged with Ubuntu).
- Remove restriction that a uniform grid is needed; use dx and dy fields, as is supported by GDAL and Surfer, but perhaps not others. - Raise a warning when rotation is used, as it currently modifies grid values via spline interpolation. Add TODO to consider removing this feature, perhaps making an example in the notes or elsewhere. Rewrite test_export_array: - First do some simpler tests on unrotated exports. - More compact reading/testing of .asc - Use np.testing.assert_almost_equal list/tuple/array checks that are not 100% identical. Other changes: - Move dx/dy to properties, which raise warnings if not constant. - Rewrite part of 'Notes' section for export_array - Raise ValueError if filename extension not supported - Add some other tests for dx/dy in test_dynamic_xll_yll
|
Mike, a few of us are currently working on refactoring some of the SpatialReference code. We are generalizing the code so that the existing SpatialReference, exporting, and plotting code can be used with MODFLOW6. Specifically, we are adding a ModelGrid class with StructuredModelGrid and VertexModelGrid subclasses that handle model grid specific behavior. SpatialReference will work with the different types of model grids to perform model grid specific operations, like get_xcenter_array. That said, I briefly looked at your changes and I believe they can all be incorporated into what we are doing. Andy Leaf may also have comments on this since he is the one currently working on the SpatialReference class. Given your interest and knowledge of the SpatialReference code I think it would be valuable if you reviewed our changes once we have completed them. Your input would be greatly appreciated. |
|
@spaulins-usgs happy to review the mf6 spatial reference code, as these are critical components for exporting a variety of formats (rasters, netcdf, vtk, etc.). And yes, it does make sense to have subclasses, such as StructuredModelGrid, to handle mesh-specific properties such as the ones that I've touched in this PR. |
|
@mwtoews, I think this looks reasonable. I agree with @spaulins-usgs that all of these changes can be integrated into the upcoming revisions (see also #368 and let us know what you think). My only concerns are potentially adding to an already somewhat cluttered list of public attributes (but maybe this isn't that big of a deal) and also dx/dy when delr/delc have variable spacings. Instead of using median values, I think they should just be |
|
As for dx/dy/geotransform returning Also, considering #368, should the names dx/dy be renamed? In most other contexts, it seems that (e.g.) "x"/"y" means a projected coordinate, while here it means local Cartesian in model units. I can alternatively "hide" these properties as |
- replace dx/dy with resolution with length multiplier applied - return None rather than median of delr/delc - geotransform raises ValueError for non-uniform grids - check full bounds for recent rasterio only
|
I removed Any other comments? |
|
Closing this PR for now as it's based on deprecated methods. It may re-appear in a different form, as it would be useful to have an Affine GeoTransform property for exporting regular grids to custom raster formats. |
The commit messages have more details, but the aim of this PR is to revise a few properties in the SpatialReference class. Here are the main changes:
dx,dy,xlength,ylengthandgeotransformproperties to SpatialReference, which are used for array exports.set_originandtheta, as they don't appear to be needed (anymore?) and have no tests. They can be re-added, but some of the trigonometry will need to be revised, and tests added.thetawas removed as it was using -ve angles for clockwise rotations. The correct rotation direction is counter clockwise (i.e. positive angles).export_array, as it is supported by Surfer and GDAL (and hence QGIS). However, not all software support dx/dy in Arc Ascii grid formats. A warning is added if either delr and/or delc are not uniform.export_arraywill do a spline interpolation. This modifies the grid dimensions and values, which are often extrapolated beyond the grid's min/max values. I've added a warning to be shown, as it's not typical that an export method would do this processing. Consider removing this feature, and moving the recipe to the Notes section or other place.In case anyone is wondering where this is all going, I'd like to eventually revise
export_netcdfto write netCDF files that are more compliant with the GDAL family of GIS software, which is best done with ageotransform.Happy to discuss any of the changes, as there is a lot going on.