-
-
- You can’t perform that action at this time.
-
-
-
-
- You signed in with another tab or window. Reload to refresh your session.
- You signed out in another tab or window. Reload to refresh your session.
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/test/meshfiles/shp/multipoly/multipoly.cpg b/test/meshfiles/shp/multipoly/multipoly.cpg
new file mode 100644
index 000000000..57decb481
--- /dev/null
+++ b/test/meshfiles/shp/multipoly/multipoly.cpg
@@ -0,0 +1 @@
+ISO-8859-1
diff --git a/test/meshfiles/shp/multipoly/multipoly.dbf b/test/meshfiles/shp/multipoly/multipoly.dbf
new file mode 100644
index 000000000..24fa8a9d3
Binary files /dev/null and b/test/meshfiles/shp/multipoly/multipoly.dbf differ
diff --git a/test/meshfiles/shp/multipoly/multipoly.shp b/test/meshfiles/shp/multipoly/multipoly.shp
new file mode 100644
index 000000000..095e5e032
Binary files /dev/null and b/test/meshfiles/shp/multipoly/multipoly.shp differ
diff --git a/test/meshfiles/shp/multipoly/multipoly.shx b/test/meshfiles/shp/multipoly/multipoly.shx
new file mode 100644
index 000000000..4d26d5647
Binary files /dev/null and b/test/meshfiles/shp/multipoly/multipoly.shx differ
diff --git a/test/test_geopandas.py b/test/test_geopandas.py
new file mode 100644
index 000000000..3c7e441e7
--- /dev/null
+++ b/test/test_geopandas.py
@@ -0,0 +1,57 @@
+import os
+import numpy as np
+
+from unittest import TestCase
+from pathlib import Path
+
+import uxarray as ux
+import matplotlib.pyplot as plt
+
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+
+class TestGeopandas(TestCase):
+
+ shp_filename = current_path / "meshfiles" / "shp" / "cb_2018_us_nation_20m" / "cb_2018_us_nation_20m.shp"
+ shp_filename_5poly = current_path / "meshfiles" / "shp" / "5poly/5poly.shp"
+ shp_filename_multi = current_path / "meshfiles" / "shp" / "multipoly/multipoly.shp"
+ geojson_filename = current_path / "meshfiles"/ "geojson"/ "sample_chicago_buildings.geojson"
+
+ def test_read_shpfile(self):
+ """Read a shapefile."""
+
+ uxgrid = ux.Grid.from_file(self.shp_filename)
+ assert (uxgrid.validate())
+
+ def test_read_shpfile_multi(self):
+ """Read a shapefile, that consists of multipolygons."""
+
+ uxgrid = ux.Grid.from_file(self.shp_filename_multi)
+ assert (uxgrid.validate())
+
+ def test_read_shpfile_5poly(self):
+ """Read a shapefile, that consists of 5 polygons of different
+ shapes."""
+
+ uxgrid = ux.Grid.from_file(self.shp_filename_5poly)
+ assert (uxgrid.validate())
+
+ def test_read_geojson(self):
+ """Read a geojson file with a few of Chicago buildings.
+
+ Number of polygons: 10
+ Polygon 1: 26 sides
+ Polygon 2: 36 sides
+ Polygon 3: 29 sides
+ Polygon 4: 10 sides
+ Polygon 5: 30 sides
+ Polygon 6: 8 sides
+ Polygon 7: 7 sides
+ Polygon 8: 9 sides
+ Polygon 9: 7 sides
+ Polygon 10: 19 sides
+ """
+
+ uxgrid = ux.Grid.from_file(self.geojson_filename)
+ assert (uxgrid.n_face == 10)
+ assert (uxgrid.n_max_face_nodes == 36)
diff --git a/test/test_grid.py b/test/test_grid.py
index 6b80aa87a..a31223ffe 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -281,7 +281,6 @@ def test_read_scrip(self):
# Test read from scrip and from ugrid for grid class
grid_CSne8 = ux.open_grid(gridfile_CSne8) # tests from scrip
- pass
class TestOperators(TestCase):
@@ -926,7 +925,6 @@ def test_from_dataset(self):
xrds = xr.open_dataset(self.gridfile_scrip)
uxgrid = ux.Grid.from_dataset(xrds)
- pass
def test_from_face_vertices(self):
single_face_latlon = [(0.0, 90.0), (-180, 0.0), (0.0, -90)]
@@ -961,7 +959,35 @@ def test_populate_bounds_GCA_mix(self):
face_bounds = bounds_xarray.values
nt.assert_allclose(grid.bounds.values, expected_bounds, atol=ERROR_TOLERANCE)
- def test_opti_bounds(self):
- import uxarray
- uxgrid = ux.open_grid(gridfile_CSne8)
- bounds = uxgrid.bounds
+ def test_populate_bounds_MPAS(self):
+ uxgrid = ux.open_grid(self.gridfile_mpas)
+ bounds_xarray = uxgrid.bounds
+
+
+class TestNormalizeExistingCoordinates(TestCase):
+ gridfile_mpas = current_path / "meshfiles" / "mpas" / "QU" / "mesh.QU.1920km.151026.nc"
+ gridfile_CSne30 = current_path / "meshfiles" / "ugrid" / "outCSne30" / "outCSne30.ug"
+
+ def test_non_norm_initial(self):
+ """Check the normalization of coordinates that were initially parsed as
+ non-normalized."""
+ from uxarray.grid.validation import _check_normalization
+ uxgrid = ux.open_grid(self.gridfile_mpas)
+
+ # Make the coordinates not normalized
+ uxgrid.node_x.data = 5 * uxgrid.node_x.data
+ uxgrid.node_y.data = 5 * uxgrid.node_y.data
+ uxgrid.node_z.data = 5 * uxgrid.node_z.data
+ assert not _check_normalization(uxgrid)
+
+ uxgrid.normalize_cartesian_coordinates()
+
+ assert _check_normalization(uxgrid)
+
+ def test_norm_initial(self):
+ """Coordinates should be normalized for grids that we construct
+ them."""
+ from uxarray.grid.validation import _check_normalization
+ uxgrid = ux.open_grid(self.gridfile_CSne30)
+
+ assert _check_normalization(uxgrid)
diff --git a/uxarray/constants.py b/uxarray/constants.py
index d55e25b49..e821de4d6 100644
--- a/uxarray/constants.py
+++ b/uxarray/constants.py
@@ -22,3 +22,5 @@
ENABLE_FMA = False
GRID_DIMS = ["n_node", "n_edge", "n_face"]
+
+WGS84_CRS = "EPSG:4326"
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index b1ff3e43c..64d6ee653 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -12,10 +12,10 @@
Union,
)
-
# reader and writer imports
from uxarray.io._exodus import _read_exodus, _encode_exodus
from uxarray.io._mpas import _read_mpas
+from uxarray.io._geopandas import _read_geodataframe
from uxarray.io._ugrid import (
_read_ugrid,
_encode_ugrid,
@@ -38,6 +38,7 @@
_set_desired_longitude_range,
_populate_node_latlon,
_populate_node_xyz,
+ _normalize_xyz,
)
from uxarray.grid.connectivity import (
_populate_edge_node_connectivity,
@@ -72,6 +73,7 @@
_check_connectivity,
_check_duplicate_nodes,
_check_area,
+ _check_normalization,
)
from xarray.core.utils import UncachedAccessor
@@ -158,6 +160,7 @@ def __init__(
# initialize attributes
self._antimeridian_face_indices = None
+ self._ds.assign_attrs({"source_grid_spec": self.source_grid_spec})
# initialize cached data structures and flags (visualization)
self._gdf = None
@@ -175,6 +178,9 @@ def __init__(
self._ball_tree = None
self._kd_tree = None
+ # flag to track if coordinates are normalized
+ self._normalized = None
+
# set desired longitude range to [-180, 180]
_set_desired_longitude_range(self._ds)
@@ -221,7 +227,9 @@ def from_dataset(
elif source_grid_spec == "ICON":
grid_ds, source_dims_dict = _read_icon(dataset, use_dual=use_dual)
elif source_grid_spec == "Shapefile":
- raise ValueError("Shapefiles not yet supported")
+ raise ValueError(
+ "Use ux.Grid.from_geodataframe(>> import uxarray as ux
+ >>> grid = ux.Grid.from_file("path/to/file.shp")
+
+ Note
+ ----
+ All formats supported by `geopandas.read_file` can be used.
+ See more at: https://geopandas.org/en/stable/docs/reference/api/geopandas.read_file.html#geopandas-read-file
+ """
+
+ # determine grid/mesh specification
+ if backend == "geopandas":
+ if str(filename).endswith(".shp"):
+ source_grid_spec = "Shapefile"
+ elif str(filename).endswith(".geojson"):
+ source_grid_spec = "GeoJSON"
+ else:
+ source_grid_spec = "OtherGeoFormat"
+
+ grid_ds, source_dims_dict = _read_geodataframe(filename)
+
+ elif backend == "xarray":
+ grid_ds, source_dims_dict = cls.from_dataset(filename)
+
+ else:
+ raise ValueError("Backend not supported")
+
+ return cls(grid_ds, source_grid_spec, source_dims_dict)
+
@classmethod
def from_topology(
cls,
@@ -1420,6 +1475,38 @@ def compute_face_areas(
return self._face_areas, self._face_jacobian
+ def normalize_cartesian_coordinates(self):
+ """Normalizes Cartesian coordinates."""
+
+ if _check_normalization(self):
+ # check if coordinates are already normalized
+ return
+
+ if "node_x" in self._ds:
+ # normalize node coordinates
+ node_x, node_y, node_z = _normalize_xyz(
+ self.node_x.values, self.node_y.values, self.node_z.values
+ )
+ self.node_x.data = node_x
+ self.node_y.data = node_y
+ self.node_z.data = node_z
+ if "edge_x" in self._ds:
+ # normalize edge coordinates
+ edge_x, edge_y, edge_z = _normalize_xyz(
+ self.edge_x.values, self.edge_y.values, self.edge_z.values
+ )
+ self.edge_x.data = edge_x
+ self.edge_y.data = edge_y
+ self.edge_z.data = edge_z
+ if "face_x" in self._ds:
+ # normalize face coordinates
+ face_x, face_y, face_z = _normalize_xyz(
+ self.face_x.values, self.face_y.values, self.face_z.values
+ )
+ self.face_x.data = face_x
+ self.face_y.data = face_y
+ self.face_z.data = face_z
+
def to_xarray(self, grid_format: Optional[str] = "ugrid"):
"""Returns a xarray Dataset representation in a specific grid format
from the Grid object.
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 81829a68a..21daab6dd 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -21,7 +21,7 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True):
Cartesian coordinates of the first GCA.
gca2_cart : [n, 3] np.ndarray where n is the number of intersection points
Cartesian coordinates of the second GCA.
- fma_disabled : bool, optional (default=False)
+ fma_disabled : bool, optional (default=True)
If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
Returns
@@ -116,7 +116,7 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=True):
def gca_constLat_intersection(
- gca_cart, constZ, fma_disabled=False, verbose=False, is_directed=False
+ gca_cart, constZ, fma_disabled=True, verbose=False, is_directed=False
):
"""Calculate the intersection point(s) of a Great Circle Arc (GCA) and a
constant latitude line in a Cartesian coordinate system.
@@ -130,7 +130,7 @@ def gca_constLat_intersection(
gca_cart : [2, 3] np.ndarray Cartesian coordinates of the two end points GCA.
constZ : float
The constant latitude represented in cartesian of the latitude line.
- fma_disabled : bool, optional (default=False)
+ fma_disabled : bool, optional (default=True)
If True, the FMA operation is disabled. And a naive `np.cross` is used instead.
verbose : bool, optional (default=False)
If True, the function prints out the intermediate results.
diff --git a/uxarray/grid/validation.py b/uxarray/grid/validation.py
index 4f74249a5..b42b8cbb3 100644
--- a/uxarray/grid/validation.py
+++ b/uxarray/grid/validation.py
@@ -1,12 +1,11 @@
import numpy as np
from warnings import warn
-
from uxarray.constants import ERROR_TOLERANCE
# validation helper functions
-def _check_connectivity(self):
+def _check_connectivity(grid):
"""Check if all nodes are referenced by at least one element.
If not, the mesh may have hanging nodes and may not a valid UGRID
@@ -15,28 +14,28 @@ def _check_connectivity(self):
# Check if all nodes are referenced by at least one element
# get unique nodes in connectivity
- nodes_in_conn = np.unique(self.face_node_connectivity.values.flatten())
+ nodes_in_conn = np.unique(grid.face_node_connectivity.values.flatten())
# remove negative indices/fill values from the list
nodes_in_conn = nodes_in_conn[nodes_in_conn >= 0]
# check if the size of unique nodes in connectivity is equal to the number of nodes
- if nodes_in_conn.size == self.n_node:
+ if nodes_in_conn.size == grid.n_node:
print("-All nodes are referenced by at least one element.")
return True
else:
warn(
"Some nodes may not be referenced by any element. {0} and {1}".format(
- nodes_in_conn.size, self.n_node
+ nodes_in_conn.size, grid.n_node
),
RuntimeWarning,
)
return False
-def _check_duplicate_nodes(self):
+def _check_duplicate_nodes(grid):
"""Check if there are duplicate nodes in the mesh."""
- coords1 = np.column_stack((np.vstack(self.node_lon), np.vstack(self.node_lat)))
+ coords1 = np.column_stack((np.vstack(grid.node_lon), np.vstack(grid.node_lat)))
unique_nodes, indices = np.unique(coords1, axis=0, return_index=True)
duplicate_indices = np.setdiff1d(np.arange(len(coords1)), indices)
@@ -53,9 +52,9 @@ def _check_duplicate_nodes(self):
return True
-def _check_area(self):
+def _check_area(grid):
"""Check if each face area is greater than our constant ERROR_TOLERANCE."""
- areas = self.face_areas
+ areas = grid.face_areas
# Check if area of any face is close to zero
if np.any(np.isclose(areas, 0, atol=ERROR_TOLERANCE)):
warn(
@@ -66,3 +65,47 @@ def _check_area(self):
else:
print("-No face area is close to zero.")
return True
+
+
+def _check_normalization(grid):
+ """Checks whether all the cartesiain coordinates are normalized."""
+
+ if grid._normalized is True:
+ # grid is already normalized, no need to run extra checks
+ return grid._normalized
+
+ if "node_x" in grid._ds:
+ if not (
+ np.isclose(
+ (grid.node_x**2 + grid.node_y**2 + grid.node_z**2),
+ 1.0,
+ atol=ERROR_TOLERANCE,
+ )
+ ).all():
+ grid._normalized = False
+ return False
+ if "edge_x" in grid._ds:
+ if not (
+ np.isclose(
+ (grid.node_x**2 + grid.node_y**2 + grid.node_z**2),
+ 1.0,
+ atol=ERROR_TOLERANCE,
+ )
+ ).all():
+ grid._normalized = False
+ return False
+ if "face_x" in grid._ds:
+ if not (
+ np.isclose(
+ (grid.node_x**2 + grid.node_y**2 + grid.node_z**2),
+ 1.0,
+ atol=ERROR_TOLERANCE,
+ )
+ ).all():
+ grid._normalized = False
+ return False
+
+ # set the grid as normalized
+ grid._normalized = True
+
+ return True
diff --git a/uxarray/io/_geopandas.py b/uxarray/io/_geopandas.py
new file mode 100644
index 000000000..4852a6d29
--- /dev/null
+++ b/uxarray/io/_geopandas.py
@@ -0,0 +1,242 @@
+import geopandas as gpd
+import xarray as xr
+import numpy as np
+from uxarray.conventions import ugrid
+from uxarray.constants import INT_DTYPE, INT_FILL_VALUE, WGS84_CRS
+
+
+def _read_geodataframe(filepath, driver=None, **kwargs):
+ """Read geospatial data, using geopandas.
+
+ Parameters
+ ----------
+ filepath : str
+ Filepath to geospatial data.
+ driver : str, optional
+ Driver to use, by default None, in which case geopandas will try to infer the driver from the file extension.
+ **kwargs
+ Keyword arguments to pass to geopandas.read_file().
+
+ Returns
+ -------
+ xr.Dataset
+ ugrid aware xarray.Dataset.
+ """
+ grid_ds = xr.Dataset()
+
+ gdf, max_coord_size = _gpd_read(filepath, driver=driver, **kwargs)
+
+ node_lon, node_lat, connectivity = _extract_geometry_info(gdf, max_coord_size)
+
+ grid_ds["node_lon"] = xr.DataArray(
+ data=node_lon, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LON_ATTRS
+ )
+ grid_ds["node_lat"] = xr.DataArray(
+ data=node_lat, dims=ugrid.NODE_DIM, attrs=ugrid.NODE_LAT_ATTRS
+ )
+ grid_ds["face_node_connectivity"] = xr.DataArray(
+ data=connectivity,
+ dims=ugrid.FACE_NODE_CONNECTIVITY_DIMS,
+ attrs=ugrid.FACE_NODE_CONNECTIVITY_ATTRS,
+ )
+
+ return grid_ds, None
+
+
+def _gpd_read(filepath, driver=None, **kwargs):
+ """Read a geospatial data using geopandas.
+
+ Parameters
+ ----------
+ filepath : str
+ Filepath to geospatial data.
+ driver : str, optional
+ Driver to use, by default None, in which case geopandas will try to infer the driver from the file extension.
+ **kwargs
+ Keyword arguments to pass to geopandas.read_file().
+
+ Returns
+ -------
+ gpd.GeoDataFrame
+ GeoDataFrame with geometries.
+ int
+ Maximum number of nodes in a polygon/multipolygon.
+ """
+ try:
+ gdf = gpd.read_file(filepath, driver=driver, **kwargs)
+ gdf = _set_crs(gdf)
+ except Exception as e:
+ print(f"An error occurred while reading the geospatial data: {e}")
+
+ max_polygon_nodes = gdf["geometry"].apply(_get_num_nodes).max()
+
+ return gdf, max_polygon_nodes
+
+
+def _set_crs(gdf):
+ """Set CRS for GeoDataFrame if not already set.
+
+ Parameters
+ ----------
+ gdf : gpd.GeoDataFrame
+ GeoDataFrame to set CRS for.
+
+ Returns
+ -------
+ gpd.GeoDataFrame
+ GeoDataFrame with CRS set.
+ """
+ if gdf.crs is None:
+ gdf = gdf.set_crs(WGS84_CRS)
+ print("Original CRS: None\nAssigned CRS:", gdf.crs)
+
+ if gdf.crs != WGS84_CRS:
+ gdf = gdf.to_crs(WGS84_CRS)
+ print("Transformed CRS:", gdf.crs)
+
+ return gdf
+
+
+def _extract_geometry_info(gdf, max_coord_size):
+ """Extract node and connectivity information from GeoDataFrame.
+
+ Parameters
+ ----------
+ gdf : gpd.GeoDataFrame
+ GeoDataFrame with geometries.
+ max_coord_size : int
+ Maximum number of nodes in a polygon/multipolygon.
+
+ Returns
+ -------
+ np.ndarray
+ Array of longitudes.
+ np.ndarray
+ Array of latitudes.
+ np.ndarray
+ Connectivity array.
+ """
+ node_lon_list = []
+ node_lat_list = []
+ connectivity = np.empty((0, max_coord_size - 1), dtype=INT_DTYPE)
+
+ node_index = 0
+
+ for _, row in gdf.iterrows():
+ geometry = row["geometry"]
+ if geometry.geom_type == "Polygon":
+ node_lat, node_lon, connectivity, node_index = _read_polygon(
+ geometry, node_lat_list, node_lon_list, connectivity, node_index
+ )
+ elif geometry.geom_type == "MultiPolygon":
+ node_lat, node_lon, connectivity, node_index = _read_multipolygon(
+ geometry, node_lat_list, node_lon_list, connectivity, node_index
+ )
+ else:
+ print(f"Unsupported geometry type: {geometry.geom_type}")
+
+ # Convert lists to numpy arrays at the end
+ node_lon = np.array(node_lon_list)
+ node_lat = np.array(node_lat_list)
+
+ return node_lon, node_lat, connectivity
+
+
+def _get_num_nodes(geom):
+ """Get number of nodes in a polygon/multipolygon.
+
+ Parameters
+ ----------
+ geom : gpd.GeoSeries
+ GeoPandas geometry object.
+
+ Returns
+ -------
+ int
+ Maximum number of nodes in a polygon or multipolygon.
+ """
+ if geom.geom_type == "Polygon":
+ return len(geom.exterior.coords)
+ elif geom.geom_type == "MultiPolygon":
+ return max(len(polygon.exterior.coords) for polygon in geom.geoms)
+ else:
+ return 0 # Not a polygon or multipolygon
+
+
+def _read_multipolygon(geometry, node_lat, node_lon, connectivity, node_index):
+ """Read a multipolygon.
+
+ Parameters
+ ----------
+ geometry : gpd.GeoSeries
+ GeoPandas geometry object.
+ node_lat : list
+ List of latitudes.
+ node_lon : list
+ List of longitudes.
+ connectivity : np.ndarray
+ Connectivity array.
+ node_index : int
+ Index of the current node.
+
+ Returns
+ -------
+ np.ndarray
+ Updated array of latitudes.
+ np.ndarray
+ Updated array of longitudes.
+ np.ndarray
+ Updated connectivity array.
+ int
+ Updated node index.
+ """
+ for polygon in geometry.geoms:
+ node_lon, node_lat, connectivity, node_index = _read_polygon(
+ polygon, node_lat, node_lon, connectivity, node_index
+ )
+ return np.array(node_lat), np.array(node_lon), connectivity, node_index
+
+
+def _read_polygon(polygon, node_lat, node_lon, connectivity, node_index):
+ """Read a polygon and append its coordinates.
+
+ Parameters
+ ----------
+ polygon : gpd.GeoSeries
+ GeoPandas geometry object.
+ node_lat : list
+ List of latitudes.
+ node_lon : list
+ List of longitudes.
+ connectivity : np.ndarray
+ Connectivity array.
+ node_index : int
+ Index of the current node.
+
+ Returns
+ -------
+ np.ndarray
+ Updated array of latitudes.
+ np.ndarray
+ Updated array of longitudes.
+ np.ndarray
+ Updated connectivity array.
+ int
+ Updated node index.
+ """
+ # Append polygon coordinates to node lists
+ node_lon.extend(polygon.exterior.coords.xy[0][:-1])
+ node_lat.extend(polygon.exterior.coords.xy[1][:-1])
+
+ coord_size_polygon = len(polygon.exterior.coords.xy[0][:-1])
+ max_coord_size = connectivity.shape[1]
+
+ new_row = np.array(range(node_index, node_index + coord_size_polygon))
+ padding_length = max_coord_size - len(new_row)
+ padding_array = np.full(padding_length, INT_FILL_VALUE)
+ new_row = np.concatenate((new_row, padding_array))
+
+ connectivity = np.vstack((connectivity, new_row))
+ node_index += coord_size_polygon
+
+ return node_lat, node_lon, connectivity, node_index
diff --git a/uxarray/io/_shapefile.py b/uxarray/io/_shapefile.py
deleted file mode 100644
index d841072ce..000000000
--- a/uxarray/io/_shapefile.py
+++ /dev/null
@@ -1,10 +0,0 @@
-def _read_shpfile(filepath):
- """Read shape file.
-
- Not implemented.
- """
- raise RuntimeError(
- "Function not implemented yet. FYI, attempted to read SHAPE file: "
- + str(filepath)
- )
- # TODO: create ds
-