diff --git a/Dockerfile b/Dockerfile index 5d6b4f9..583973e 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ FROM python:3.11 RUN apt-get update -y ; apt-get install -y nano netcdf-bin -RUN pip install xarray netcdf4 numpy scipy pymongo shapely +RUN pip install xarray netcdf4 numpy scipy pymongo shapely geopy diff --git a/gridtools.py b/gridtools.py index 4453552..3671a5b 100644 --- a/gridtools.py +++ b/gridtools.py @@ -1,16 +1,30 @@ -import scipy, numpy, copy +import scipy, numpy, copy, math def label_features(feature_map, structure=[[1,1,1],[1,1,1],[1,1,1]], connected_poles=True, periodic_dateline=True): # given a 2D numpy array feature_map[latitude][longitude] labeling features with 1 and voids with 0, # label distinct isolated features. # periodic_dateline = True makes a periodic boundary on the inner index # connected_poles = True makes all features touching a pole connected + labeled_map = scipy.ndimage.label(feature_map, structure=structure)[0] + # periodic boundary if periodic_dateline: - for y in range(labeled_map.shape[0]): - if labeled_map[y, 0] > 0 and labeled_map[y, -1] > 0: - labeled_map[labeled_map == labeled_map[y, -1]] = labeled_map[y, 0] + if structure == [[0,1,0],[1,1,1],[0,1,0]]: + # no diag + for y in range(labeled_map.shape[0]): + if labeled_map[y, 0] > 0 and labeled_map[y, -1] > 0: + labeled_map[labeled_map == labeled_map[y, -1]] = labeled_map[y, 0] + elif structure == [[1,1,1],[1,1,1],[1,1,1]]: + # diagonally connected + for y in range(labeled_map.shape[0]): + if labeled_map[y, 0] > 0 and labeled_map[y, -1] > 0: + labeled_map[labeled_map == labeled_map[y, -1]] = labeled_map[y, 0] + elif labeled_map[y, 0] > 0 and labeled_map[max(y-1,0), -1] > 0: + labeled_map[labeled_map == labeled_map[max(y-1,0), -1]] = labeled_map[y, 0] + elif labeled_map[y, 0] > 0 and labeled_map[min(y+1,labeled_map.shape[0]-1), -1] > 0: + labeled_map[labeled_map == labeled_map[min(y+1,labeled_map.shape[0]-1), -1]] = labeled_map[y, 0] + # connected poles if connected_poles: spole_features = [x for x in numpy.unique(labeled_map[0]) if x > 0] @@ -24,15 +38,21 @@ def label_features(feature_map, structure=[[1,1,1],[1,1,1],[1,1,1]], connected_p return labeled_map -def trace_shape(labeled_map, label, nlatsteps=361): - # trace the shape labeled with