From e61f2dd78f8f79be96a903523bc5d966bbe00051 Mon Sep 17 00:00:00 2001 From: Nicholas Delli Carpini Date: Mon, 13 Nov 2023 14:21:29 -0500 Subject: [PATCH 1/5] fix regular, hycom, and fvcom grids --- xpublish_wms/grid.py | 252 +++++++++++++++++++++------ xpublish_wms/utils.py | 16 ++ xpublish_wms/wms/get_feature_info.py | 13 +- xpublish_wms/wms/get_map.py | 7 +- 4 files changed, 222 insertions(+), 66 deletions(-) diff --git a/xpublish_wms/grid.py b/xpublish_wms/grid.py index a6c17d1..158ae57 100644 --- a/xpublish_wms/grid.py +++ b/xpublish_wms/grid.py @@ -11,7 +11,7 @@ import xarray as xr from sklearn.neighbors import BallTree -from xpublish_wms.utils import strip_float, to_mercator +from xpublish_wms.utils import strip_float, to_mercator, lnglat_to_mercator class RenderMethod(Enum): @@ -130,7 +130,10 @@ def sel_lat_lng( parameters, ) -> Tuple[xr.Dataset, list, list]: """Select the given dataset by the given lon/lat and optional elevation""" + + subset = self.mask(subset) subset = subset.cf.interp(longitude=lng, latitude=lat) + x_axis = [strip_float(subset.cf["longitude"])] y_axis = [strip_float(subset.cf["latitude"])] return subset, x_axis, y_axis @@ -157,10 +160,43 @@ def crs(self) -> str: return "EPSG:4326" def project(self, da: xr.DataArray, crs: str) -> xr.DataArray: - if crs == "EPSG:4326": - da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) - elif crs == "EPSG:3857": - da = da.rio.reproject("EPSG:3857") + da = self.mask(da) + + coords = dict() + # need to convert longitude and latitude to x and y for the mesh to work properly + # regular grid doesn't have x & y as dimensions so have to remake the whole data array + for coord in da.coords: + if coord != da.cf["longitude"].name and coord != da.cf["latitude"].name: + coords[coord] = da.coords[coord] + + # build new x coordinate + coords["x"] = ( + "x", + da.cf["longitude"].values, + da.cf["longitude"].attrs + ) + # build new y coordinate + coords["y"] = ( + "y", + da.cf["latitude"].values, + da.cf["latitude"].attrs + ) + # build new data array + da = xr.DataArray( + data=da, + dims=("y", "x"), + coords=coords, + name=da.name, + attrs=da.attrs, + ) + + # convert to mercator + if crs == "EPSG:3857": + lng, lat = lnglat_to_mercator(da.cf["longitude"], da.cf["latitude"]) + + da = da.assign_coords({"x": lng, "y": lat}) + da = da.unify_chunks() + return da @@ -188,6 +224,8 @@ def crs(self) -> str: return "EPSG:4326" def project(self, da: xr.DataArray, crs: str) -> xr.DataArray: + da = self.mask(da) + if crs == "EPSG:4326": da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) elif crs == "EPSG:3857": @@ -221,6 +259,9 @@ def sel_lat_lng( parameters, ) -> Tuple[xr.Dataset, list, list]: """Select the given dataset by the given lon/lat and optional elevation""" + + subset = self.mask(subset) + subset = sel2d( subset[parameters], lons=subset.cf["longitude"], @@ -258,14 +299,21 @@ def mask( da: Union[xr.DataArray, xr.Dataset], ) -> Union[xr.DataArray, xr.Dataset]: mask = self.ds[f'mask_{da.cf["latitude"].name.split("_")[1]}'] - mask = mask.cf.isel(time=0).squeeze(drop=True).cf.drop_vars("time") + if "time" in mask.cf.coords: + mask = mask.cf.isel(time=0).squeeze(drop=True).cf.drop_vars("time") + else: + mask = mask.cf.squeeze(drop=True).cf + mask[:-1, :] = mask[:-1, :].where(mask[1:, :] == 1, 0) mask[:, :-1] = mask[:, :-1].where(mask[:, 1:] == 1, 0) mask[1:, :] = mask[1:, :].where(mask[:-1, :] == 1, 0) mask[:, 1:] = mask[:, 1:].where(mask[:, :-1] == 1, 0) + return da.where(mask == 1) def project(self, da: xr.DataArray, crs: str) -> xr.DataArray: + da = self.mask(da) + if crs == "EPSG:4326": da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) elif crs == "EPSG:3857": @@ -311,10 +359,11 @@ def sel_lat_lng( " ", ) + new_selected_ds = self.mask(subset[[parameter]]) new_selected_ds = sel2d( - subset, - lons=subset.cf[lng_coord], - lats=subset.cf[lat_coord], + new_selected_ds, + lons=new_selected_ds.cf[lng_coord], + lats=new_selected_ds.cf[lat_coord], lon0=lng, lat0=lat, ) @@ -371,34 +420,94 @@ def bbox(self, da: xr.DataArray) -> Tuple[float, float, float, float]: float(da.cf["latitude"].max()), ) + def mask( + self, + da: Union[xr.DataArray, xr.Dataset], + ) -> Union[xr.DataArray, xr.Dataset]: + # mask values where the longitude is >500 bc of RTOFS-Global https://oceanpython.org/2012/11/29/global-rtofs-real-time-ocean-forecast-system/ + lng = da.cf["longitude"] + + # add missing dimensions to lng mask (ie. elevation & time) + if len(da.dims) > len(lng.dims): + for dim in set(da.dims).symmetric_difference(lng.dims): + lng = lng.expand_dims({dim: da[dim]}) + + mask = lng > 500 + + np_mask = np.ma.masked_where(mask == True, da.values)[:] + np_mask[:-1, :] = np.ma.masked_where(mask[1:, :] == True, np_mask[:-1, :])[:] + np_mask[:, :-1] = np.ma.masked_where(mask[:, 1:] == True, np_mask[:, :-1])[:] + np_mask[1:, :] = np.ma.masked_where(mask[:-1, :] == True, np_mask[1:, :])[:] + np_mask[:, 1:] = np.ma.masked_where(mask[:, :-1] == True, np_mask[:, 1:])[:] + + return da.where(np_mask.mask == False) + def project(self, da: xr.DataArray, crs: str) -> Any: - # TODO: Figure out global coords + da = self.mask(da) + + # create 2 separate DataArrays where points lng>180 are put at the beginning of the array + mask_0 = xr.where(da.cf["longitude"] <= 180, True, False) + temp_da_0 = da.where(mask_0.compute() == True, drop=True) + da_0 = xr.DataArray( + data=temp_da_0, + dims=temp_da_0.dims, + name=temp_da_0.name, + coords=temp_da_0.coords, + attrs=temp_da_0.attrs + ) + + mask_1 = xr.where(da.cf["longitude"] > 180, True, False) + temp_da_1 = da.where(mask_1.compute() == True, drop=True) + temp_da_1.cf["longitude"][:] = temp_da_1.cf["longitude"][:] - 360 + da_1 = xr.DataArray( + data=temp_da_1, + dims=temp_da_1.dims, + name=temp_da_1.name, + coords=temp_da_1.coords, + attrs=temp_da_1.attrs, + ) + + # put the 2 DataArrays back together in the proper order + da = xr.concat([da_1, da_0], dim="X") + if crs == "EPSG:4326": da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) elif crs == "EPSG:3857": - x, y = to_mercator.transform(da.cf["longitude"], da.cf["latitude"]) - x_chunks = ( - da.cf["longitude"].chunks if da.cf["longitude"].chunks else x.shape - ) - y_chunks = da.cf["latitude"].chunks if da.cf["latitude"].chunks else y.shape - - da = da.assign_coords( - { - "x": ( - da.cf["longitude"].dims, - dask_array.from_array(x, chunks=x_chunks), - ), - "y": ( - da.cf["latitude"].dims, - dask_array.from_array(y, chunks=y_chunks), - ), - }, - ) + lng, lat = lnglat_to_mercator(da.cf["longitude"], da.cf["latitude"]) + da = da.assign_coords({"x": lng, "y": lat}) da = da.unify_chunks() + return da + def sel_lat_lng( + self, + subset: xr.Dataset, + lng, + lat, + parameters, + ) -> Tuple[xr.Dataset, list, list]: + """Select the given dataset by the given lon/lat and optional elevation""" + + for parameter in parameters: + subset[parameter] = self.mask(subset[parameter]) + + subset.cf["longitude"][:] = xr.where(subset.cf["longitude"] > 180, subset.cf["longitude"] - 360, subset.cf["longitude"])[:] + + subset = sel2d( + subset[parameters], + lons=subset.cf["longitude"], + lats=subset.cf["latitude"], + lon0=lng, + lat0=lat, + ) + + x_axis = [strip_float(subset.cf["longitude"])] + y_axis = [strip_float(subset.cf["latitude"])] + return subset, x_axis, y_axis + + class FVCOMGrid(Grid): def __init__(self, ds: xr.Dataset): self.ds = ds @@ -478,6 +587,8 @@ def sel_lat_lng( ) -> Tuple[xr.Dataset, list, list]: """Select the given dataset by the given lon/lat and optional elevation""" + subset = self.mask(subset) + lng_rad = np.deg2rad(subset.cf["longitude"]) lat_rad = np.deg2rad(subset.cf["latitude"]) @@ -485,7 +596,7 @@ def sel_lat_lng( tree = BallTree(stacked, leaf_size=2, metric="haversine") idx = tree.query( - [[np.deg2rad((360 + lng) if lng < 0 else lng), np.deg2rad(lat)]], + [[np.deg2rad(360 + lng), np.deg2rad(lat)]], return_distance=False, ) @@ -534,38 +645,51 @@ def select_by_elevation( return da + def project(self, da: xr.DataArray, crs: str) -> Any: - # fvcom nodal variables have data on both the faces and edges + da = self.mask(da) + + data = da.values + # create new data by getting values from the surrounding edges if "nele" in da.dims: - elem_count = self.ds.ntve.isel(time=0).values - neighbors = self.ds.nbve.isel(time=0).values + elem_count = self.ds.ntve.isel(time=0).values if "time" in self.ds.ntve.coords else self.ds.ntve.values + neighbors = self.ds.nbve.isel(time=0).values if "time" in self.ds.nbve.coords else self.ds.nbve.values mask = neighbors[:, :] > 0 - new_values = ( - np.sum(da.values[neighbors[:, :] - 1], axis=0, where=mask) / elem_count - ) - da = xr.DataArray( - data=new_values, - dims=da.dims, - name=da.name, - attrs=da.attrs, - coords=dict( - lonc=( - da.cf["longitude"].dims, - self.ds.lon.values, - da.cf["longitude"].attrs, - ), - latc=( - da.cf["latitude"].dims, - self.ds.lat.values, - da.cf["latitude"].attrs, - ), - time=da.coords["time"], - ), + data = ( + np.sum(data[neighbors[:, :] - 1], axis=0, where=mask) / elem_count ) + coords = dict() + # need to create new x & y coordinates with dataset values while dropping the old ones + # can't keep the original values or else da.cf will have 2 lng/lat arrays + for coord in da.coords: + if coord != da.cf["longitude"].name and coord != da.cf["latitude"].name: + coords[coord] = da.coords[coord] + + # build new x coordinate + coords["x"] = ( + da.cf["longitude"].dims, + self.ds.lon.values, + da.cf["longitude"].attrs + ) + # build new y coordinate + coords["y"] = ( + da.cf["latitude"].dims, + self.ds.lat.values, + da.cf["latitude"].attrs + ) + # build new data array + da = xr.DataArray( + data=data, + dims=da.dims, + name=da.name, + attrs=da.attrs, + coords=coords, + ) + if crs == "EPSG:4326": - da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) + da.__setitem__(da.cf["longitude"].name, da.cf["longitude"] - 360) elif crs == "EPSG:3857": x, y = to_mercator.transform(da.cf["longitude"], da.cf["latitude"]) x_chunks = ( @@ -578,10 +702,12 @@ def project(self, da: xr.DataArray, crs: str) -> Any: "x": ( da.cf["longitude"].dims, dask_array.from_array(x, chunks=x_chunks), + da.cf["longitude"].attrs, ), "y": ( da.cf["latitude"].dims, dask_array.from_array(y, chunks=y_chunks), + da.cf["latitude"].attrs, ), }, ) @@ -590,10 +716,15 @@ def project(self, da: xr.DataArray, crs: str) -> Any: return da def tessellate(self, da: xr.DataArray) -> np.ndarray: + nv = self.ds.nv + if len(nv.shape) > 2: + for i in range(len(nv.shape) - 2): + nv = nv[0] + return tri.Triangulation( da.cf["longitude"], da.cf["latitude"], - self.ds.nv[0].T - 1, + nv.T - 1, ).triangles @@ -660,6 +791,8 @@ def sel_lat_lng( ) -> Tuple[xr.Dataset, list, list]: """Select the given dataset by the given lon/lat and optional elevation""" + subset = self.mask(subset) + lng_rad = np.deg2rad(subset.cf["longitude"]) lat_rad = np.deg2rad(subset.cf["latitude"]) @@ -727,6 +860,8 @@ def select_by_elevation( return da def project(self, da: xr.DataArray, crs: str) -> Any: + da = self.mask(da) + if crs == "EPSG:4326": da = da.assign_coords({"x": da.cf["longitude"], "y": da.cf["latitude"]}) elif crs == "EPSG:3857": @@ -753,10 +888,15 @@ def project(self, da: xr.DataArray, crs: str) -> Any: return da def tessellate(self, da: xr.DataArray) -> np.ndarray: + ele = self.ds.ele + if len(ele.shape) > 2: + for i in range(len(ele.shape) - 2): + ele = ele[0] + return tri.Triangulation( da.cf["longitude"], da.cf["latitude"], - self.ds.ele[0].T - 1, + ele.T - 1, ).triangles diff --git a/xpublish_wms/utils.py b/xpublish_wms/utils.py index bf55f29..d1932f1 100644 --- a/xpublish_wms/utils.py +++ b/xpublish_wms/utils.py @@ -1,4 +1,5 @@ import logging +import math from typing import Union import numpy as np @@ -79,6 +80,21 @@ def lnglat_to_cartesian(longitude, latitude): z = R * np.sin(lat_rad) return np.column_stack((x, y, z)) +def lnglat_to_mercator(longitude, latitude): + """ + Converts data array with cf standard lng/lat to mercator coordinates + """ + constant = 20037508.34 / 180 + + longitude = xr.where(longitude == 180, longitude - 0.000001, longitude) + longitude = xr.where(longitude == -180, longitude + 0.000001, longitude) + longitude = longitude * constant + + latitude = xr.where(latitude == 90, latitude - 0.000001, latitude) + latitude = xr.where(latitude == -90, latitude + 0.000001, latitude) + latitude = (np.log(np.tan((90 + latitude) * math.pi / 360)) / (math.pi / 180)) * constant + + return longitude, latitude to_lnglat = Transformer.from_crs(3857, 4326, always_xy=True) diff --git a/xpublish_wms/wms/get_feature_info.py b/xpublish_wms/wms/get_feature_info.py index 8ec27c1..c8ed66d 100644 --- a/xpublish_wms/wms/get_feature_info.py +++ b/xpublish_wms/wms/get_feature_info.py @@ -147,7 +147,9 @@ def get_feature_info(ds: xr.Dataset, query: dict) -> Response: y_coord = np.linspace(bbox[1], bbox[3], height) if any_has_time_axis: - if len(times) == 1: + if len(selected_ds.cf["time"].shape) == 1 and selected_ds.cf["time"].shape[0] == 1: + selected_ds = selected_ds.cf.isel(time=0) + elif len(times) == 1: selected_ds = selected_ds.cf.interp(time=times[0]) elif len(times) > 1: selected_ds = selected_ds.cf.sel(time=slice(times[0], times[1])) @@ -164,8 +166,6 @@ def get_feature_info(ds: xr.Dataset, query: dict) -> Response: selected_ds = ds.gridded.select_by_elevation(selected_ds, elevation) try: - # Apply masking if necessary - selected_ds = ds.gridded.mask(selected_ds) selected_ds, x_axis, y_axis = ds.gridded.sel_lat_lng( subset=selected_ds, lng=x_coord[x], @@ -176,9 +176,10 @@ def get_feature_info(ds: xr.Dataset, query: dict) -> Response: raise HTTPException(500, f"Error with grid type {ds.gridded.name}: {e}") # When none of the parameters have data, drop it - time_coord_name = selected_ds.cf.coordinates["time"][0] - if any_has_time_axis and selected_ds[time_coord_name].shape: - selected_ds = selected_ds.dropna(time_coord_name, how="all") + if any_has_time_axis: + time_coord_name = selected_ds.cf.coordinates["time"][0] + if selected_ds[time_coord_name].shape: + selected_ds = selected_ds.dropna(time_coord_name, how="all") if not any_has_time_axis: t_axis = None diff --git a/xpublish_wms/wms/get_map.py b/xpublish_wms/wms/get_map.py index dbee864..0060385 100644 --- a/xpublish_wms/wms/get_map.py +++ b/xpublish_wms/wms/get_map.py @@ -122,14 +122,14 @@ def ensure_query_types(self, ds: xr.Dataset, query: dict): self.time = pd.to_datetime(self.time_str).tz_localize(None) else: self.time = None - self.has_time = "time" in ds[self.parameter].cf.coordinates + self.has_time = self.TIME_CF_NAME in ds[self.parameter].cf.coords self.elevation_str = query.get("elevation", None) if self.elevation_str: self.elevation = float(self.elevation_str) else: self.elevation = None - self.has_elevation = "vertical" in ds[self.parameter].cf.coordinates + self.has_elevation = self.ELEVATION_CF_NAME in ds[self.parameter].cf.coords # Grid self.crs = query.get("crs", None) or query.get("srs") @@ -180,7 +180,7 @@ def select_time(self, da: xr.DataArray) -> xr.DataArray: """ if self.time is not None: da = da.cf.sel({self.TIME_CF_NAME: self.time}, method="nearest") - else: + elif self.TIME_CF_NAME in da.cf.coords: da = da.cf.isel({self.TIME_CF_NAME: -1}) return da @@ -240,7 +240,6 @@ def render( # TODO: FVCOM and other grids # return self.render_quad_grid(da, buffer, minmax_only) projection_start = time.time() - da = ds.gridded.mask(da) da = ds.gridded.project(da, self.crs) logger.debug(f"Projection time: {time.time() - projection_start}") From 2bf9fe97df87034481d7f647124ceff284fdefb0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 Nov 2023 19:26:24 +0000 Subject: [PATCH 2/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xpublish_wms/grid.py | 44 ++++++++++++++-------------- xpublish_wms/utils.py | 6 +++- xpublish_wms/wms/get_feature_info.py | 5 +++- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/xpublish_wms/grid.py b/xpublish_wms/grid.py index 158ae57..83d2263 100644 --- a/xpublish_wms/grid.py +++ b/xpublish_wms/grid.py @@ -11,7 +11,7 @@ import xarray as xr from sklearn.neighbors import BallTree -from xpublish_wms.utils import strip_float, to_mercator, lnglat_to_mercator +from xpublish_wms.utils import lnglat_to_mercator, strip_float, to_mercator class RenderMethod(Enum): @@ -170,17 +170,9 @@ def project(self, da: xr.DataArray, crs: str) -> xr.DataArray: coords[coord] = da.coords[coord] # build new x coordinate - coords["x"] = ( - "x", - da.cf["longitude"].values, - da.cf["longitude"].attrs - ) + coords["x"] = ("x", da.cf["longitude"].values, da.cf["longitude"].attrs) # build new y coordinate - coords["y"] = ( - "y", - da.cf["latitude"].values, - da.cf["latitude"].attrs - ) + coords["y"] = ("y", da.cf["latitude"].values, da.cf["latitude"].attrs) # build new data array da = xr.DataArray( data=da, @@ -453,7 +445,7 @@ def project(self, da: xr.DataArray, crs: str) -> Any: dims=temp_da_0.dims, name=temp_da_0.name, coords=temp_da_0.coords, - attrs=temp_da_0.attrs + attrs=temp_da_0.attrs, ) mask_1 = xr.where(da.cf["longitude"] > 180, True, False) @@ -480,7 +472,6 @@ def project(self, da: xr.DataArray, crs: str) -> Any: return da - def sel_lat_lng( self, subset: xr.Dataset, @@ -493,7 +484,11 @@ def sel_lat_lng( for parameter in parameters: subset[parameter] = self.mask(subset[parameter]) - subset.cf["longitude"][:] = xr.where(subset.cf["longitude"] > 180, subset.cf["longitude"] - 360, subset.cf["longitude"])[:] + subset.cf["longitude"][:] = xr.where( + subset.cf["longitude"] > 180, + subset.cf["longitude"] - 360, + subset.cf["longitude"], + )[:] subset = sel2d( subset[parameters], @@ -645,20 +640,25 @@ def select_by_elevation( return da - def project(self, da: xr.DataArray, crs: str) -> Any: da = self.mask(da) data = da.values # create new data by getting values from the surrounding edges if "nele" in da.dims: - elem_count = self.ds.ntve.isel(time=0).values if "time" in self.ds.ntve.coords else self.ds.ntve.values - neighbors = self.ds.nbve.isel(time=0).values if "time" in self.ds.nbve.coords else self.ds.nbve.values + elem_count = ( + self.ds.ntve.isel(time=0).values + if "time" in self.ds.ntve.coords + else self.ds.ntve.values + ) + neighbors = ( + self.ds.nbve.isel(time=0).values + if "time" in self.ds.nbve.coords + else self.ds.nbve.values + ) mask = neighbors[:, :] > 0 - data = ( - np.sum(data[neighbors[:, :] - 1], axis=0, where=mask) / elem_count - ) + data = np.sum(data[neighbors[:, :] - 1], axis=0, where=mask) / elem_count coords = dict() # need to create new x & y coordinates with dataset values while dropping the old ones @@ -671,13 +671,13 @@ def project(self, da: xr.DataArray, crs: str) -> Any: coords["x"] = ( da.cf["longitude"].dims, self.ds.lon.values, - da.cf["longitude"].attrs + da.cf["longitude"].attrs, ) # build new y coordinate coords["y"] = ( da.cf["latitude"].dims, self.ds.lat.values, - da.cf["latitude"].attrs + da.cf["latitude"].attrs, ) # build new data array da = xr.DataArray( diff --git a/xpublish_wms/utils.py b/xpublish_wms/utils.py index d1932f1..e965399 100644 --- a/xpublish_wms/utils.py +++ b/xpublish_wms/utils.py @@ -80,6 +80,7 @@ def lnglat_to_cartesian(longitude, latitude): z = R * np.sin(lat_rad) return np.column_stack((x, y, z)) + def lnglat_to_mercator(longitude, latitude): """ Converts data array with cf standard lng/lat to mercator coordinates @@ -92,10 +93,13 @@ def lnglat_to_mercator(longitude, latitude): latitude = xr.where(latitude == 90, latitude - 0.000001, latitude) latitude = xr.where(latitude == -90, latitude + 0.000001, latitude) - latitude = (np.log(np.tan((90 + latitude) * math.pi / 360)) / (math.pi / 180)) * constant + latitude = ( + np.log(np.tan((90 + latitude) * math.pi / 360)) / (math.pi / 180) + ) * constant return longitude, latitude + to_lnglat = Transformer.from_crs(3857, 4326, always_xy=True) diff --git a/xpublish_wms/wms/get_feature_info.py b/xpublish_wms/wms/get_feature_info.py index c8ed66d..a1edfff 100644 --- a/xpublish_wms/wms/get_feature_info.py +++ b/xpublish_wms/wms/get_feature_info.py @@ -147,7 +147,10 @@ def get_feature_info(ds: xr.Dataset, query: dict) -> Response: y_coord = np.linspace(bbox[1], bbox[3], height) if any_has_time_axis: - if len(selected_ds.cf["time"].shape) == 1 and selected_ds.cf["time"].shape[0] == 1: + if ( + len(selected_ds.cf["time"].shape) == 1 + and selected_ds.cf["time"].shape[0] == 1 + ): selected_ds = selected_ds.cf.isel(time=0) elif len(times) == 1: selected_ds = selected_ds.cf.interp(time=times[0]) From 7ab6cd159e62059c15b87f7e4d6500efcc4b1c53 Mon Sep 17 00:00:00 2001 From: Nicholas Delli Carpini Date: Mon, 13 Nov 2023 14:35:00 -0500 Subject: [PATCH 3/5] formatting --- xpublish_wms/grid.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/xpublish_wms/grid.py b/xpublish_wms/grid.py index 83d2263..3acbeea 100644 --- a/xpublish_wms/grid.py +++ b/xpublish_wms/grid.py @@ -484,11 +484,11 @@ def sel_lat_lng( for parameter in parameters: subset[parameter] = self.mask(subset[parameter]) - subset.cf["longitude"][:] = xr.where( + subset.__setitem__(subset.cf["longitude"].name, xr.where( subset.cf["longitude"] > 180, subset.cf["longitude"] - 360, subset.cf["longitude"], - )[:] + )) subset = sel2d( subset[parameters], @@ -646,18 +646,17 @@ def project(self, da: xr.DataArray, crs: str) -> Any: data = da.values # create new data by getting values from the surrounding edges if "nele" in da.dims: - elem_count = ( - self.ds.ntve.isel(time=0).values - if "time" in self.ds.ntve.coords - else self.ds.ntve.values - ) - neighbors = ( - self.ds.nbve.isel(time=0).values - if "time" in self.ds.nbve.coords - else self.ds.nbve.values - ) - mask = neighbors[:, :] > 0 + if "time" in self.ds.ntve.coords: + elem_count = self.ds.ntve.isel(time=0).values + else: + elem_count = self.ds.ntve.values + if "time" in self.ds.nbve.coords: + neighbors = self.ds.nbve.isel(time=0).values + else: + neighbors = self.ds.nbve.values + + mask = neighbors[:, :] > 0 data = np.sum(data[neighbors[:, :] - 1], axis=0, where=mask) / elem_count coords = dict() From e0555ec5a0ae973ee495b617f2a9e6b55879f1d1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 Nov 2023 19:35:13 +0000 Subject: [PATCH 4/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xpublish_wms/grid.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/xpublish_wms/grid.py b/xpublish_wms/grid.py index 3acbeea..718753f 100644 --- a/xpublish_wms/grid.py +++ b/xpublish_wms/grid.py @@ -484,11 +484,14 @@ def sel_lat_lng( for parameter in parameters: subset[parameter] = self.mask(subset[parameter]) - subset.__setitem__(subset.cf["longitude"].name, xr.where( - subset.cf["longitude"] > 180, - subset.cf["longitude"] - 360, - subset.cf["longitude"], - )) + subset.__setitem__( + subset.cf["longitude"].name, + xr.where( + subset.cf["longitude"] > 180, + subset.cf["longitude"] - 360, + subset.cf["longitude"], + ), + ) subset = sel2d( subset[parameters], From 812a7f06f89c2c6b1ffe41c9b0f557e0a26e8310 Mon Sep 17 00:00:00 2001 From: Nicholas Delli Carpini Date: Mon, 13 Nov 2023 15:05:36 -0500 Subject: [PATCH 5/5] fix mask typing --- xpublish_wms/grid.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/xpublish_wms/grid.py b/xpublish_wms/grid.py index 718753f..899885c 100644 --- a/xpublish_wms/grid.py +++ b/xpublish_wms/grid.py @@ -426,20 +426,20 @@ def mask( mask = lng > 500 - np_mask = np.ma.masked_where(mask == True, da.values)[:] - np_mask[:-1, :] = np.ma.masked_where(mask[1:, :] == True, np_mask[:-1, :])[:] - np_mask[:, :-1] = np.ma.masked_where(mask[:, 1:] == True, np_mask[:, :-1])[:] - np_mask[1:, :] = np.ma.masked_where(mask[:-1, :] == True, np_mask[1:, :])[:] - np_mask[:, 1:] = np.ma.masked_where(mask[:, :-1] == True, np_mask[:, 1:])[:] + np_mask = np.ma.masked_where(mask == 1, da.values)[:] + np_mask[:-1, :] = np.ma.masked_where(mask[1:, :] == 1, np_mask[:-1, :])[:] + np_mask[:, :-1] = np.ma.masked_where(mask[:, 1:] == 1, np_mask[:, :-1])[:] + np_mask[1:, :] = np.ma.masked_where(mask[:-1, :] == 1, np_mask[1:, :])[:] + np_mask[:, 1:] = np.ma.masked_where(mask[:, :-1] == 1, np_mask[:, 1:])[:] - return da.where(np_mask.mask == False) + return da.where(np_mask.mask == 0) def project(self, da: xr.DataArray, crs: str) -> Any: da = self.mask(da) # create 2 separate DataArrays where points lng>180 are put at the beginning of the array - mask_0 = xr.where(da.cf["longitude"] <= 180, True, False) - temp_da_0 = da.where(mask_0.compute() == True, drop=True) + mask_0 = xr.where(da.cf["longitude"] <= 180, 1, 0) + temp_da_0 = da.where(mask_0.compute() == 1, drop=True) da_0 = xr.DataArray( data=temp_da_0, dims=temp_da_0.dims, @@ -448,8 +448,8 @@ def project(self, da: xr.DataArray, crs: str) -> Any: attrs=temp_da_0.attrs, ) - mask_1 = xr.where(da.cf["longitude"] > 180, True, False) - temp_da_1 = da.where(mask_1.compute() == True, drop=True) + mask_1 = xr.where(da.cf["longitude"] > 180, 1, 0) + temp_da_1 = da.where(mask_1.compute() == 1, drop=True) temp_da_1.cf["longitude"][:] = temp_da_1.cf["longitude"][:] - 360 da_1 = xr.DataArray( data=temp_da_1,