Skip to content
Snippets Groups Projects
Commit 210d1079 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Blacken

parent 8d5355a9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 88 additions and 54 deletions
Common/util_eddies.py
+ 88
− 54
View file @ 210d1079
...@@ -7,6 +7,7 @@ import os ...@@ -7,6 +7,7 @@ import os
import shapefile import shapefile
import numpy as np import numpy as np
class SHPC_class: class SHPC_class:
"""We are assuming that eddies are indexed from 1, without jump, for """We are assuming that eddies are indexed from 1, without jump, for
each date. each date.
...@@ -21,7 +22,7 @@ class SHPC_class: ...@@ -21,7 +22,7 @@ class SHPC_class:
""" """
if not os.access(SHPC_dir, os.R_OK): if not os.access(SHPC_dir, os.R_OK):
print(f"{SHPC_dir} not readable", file = sys.stderr) print(f"{SHPC_dir} not readable", file=sys.stderr)
raise OSError raise OSError
self.dir = SHPC_dir self.dir = SHPC_dir
...@@ -29,14 +30,18 @@ class SHPC_class: ...@@ -29,14 +30,18 @@ class SHPC_class:
fname = path.join(SHPC_dir, "n_slices.txt") fname = path.join(SHPC_dir, "n_slices.txt")
try: try:
with open(fname) as f: self.n_slices = int(f.read()) with open(fname) as f:
self.n_slices = int(f.read())
except FileNotFoundError: except FileNotFoundError:
print(fname, "not found. Creating it...") print(fname, "not found. Creating it...")
self.n_slices = 0 self.n_slices = 0
while True: while True:
slice_dir = path.join(SHPC_dir, f"Slice_{self.n_slices}") slice_dir = path.join(SHPC_dir, f"Slice_{self.n_slices}")
if not path.isdir(slice_dir): break
if not path.isdir(slice_dir):
break
self.n_slices += 1 self.n_slices += 1
if self.n_slices == 0: if self.n_slices == 0:
...@@ -46,20 +51,23 @@ class SHPC_class: ...@@ -46,20 +51,23 @@ class SHPC_class:
with open(fname, "w", encoding="utf-8") as f: with open(fname, "w", encoding="utf-8") as f:
f.write(f"{self.n_slices}\n") f.write(f"{self.n_slices}\n")
self._readers = [{"Anticyclones": {}, "Cyclones": {}} self._readers = [
for i in range(self.n_slices)] {"Anticyclones": {}, "Cyclones": {}} for i in range(self.n_slices)
]
self.d_init = [] self.d_init = []
for i_slice, reader in enumerate(self._readers): for i_slice, reader in enumerate(self._readers):
my_shapefile = path.join(SHPC_dir, f"Slice_{i_slice}", def_orient, my_shapefile = path.join(
"extremum") SHPC_dir, f"Slice_{i_slice}", def_orient, "extremum"
)
try: try:
reader[def_orient]["extremum"] = shapefile.Reader(my_shapefile) reader[def_orient]["extremum"] = shapefile.Reader(my_shapefile)
except shapefile.ShapefileException: except shapefile.ShapefileException:
# Use the other orientation for this slice: # Use the other orientation for this slice:
orientation = other_orient(def_orient) orientation = other_orient(def_orient)
my_shapefile = path.join(SHPC_dir, f"Slice_{i_slice}", my_shapefile = path.join(
orientation, "extremum") SHPC_dir, f"Slice_{i_slice}", orientation, "extremum"
)
reader[orientation]["extremum"] = shapefile.Reader(my_shapefile) reader[orientation]["extremum"] = shapefile.Reader(my_shapefile)
else: else:
orientation = def_orient orientation = def_orient
...@@ -67,22 +75,29 @@ class SHPC_class: ...@@ -67,22 +75,29 @@ class SHPC_class:
if len(reader[orientation]["extremum"]) != 0: if len(reader[orientation]["extremum"]) != 0:
d_init_slice = reader[orientation]["extremum"].record(0).date d_init_slice = reader[orientation]["extremum"].record(0).date
else: else:
print(f"__init__: Empty shapefile: {my_shapefile}", print(
file = sys.stderr) f"__init__: Empty shapefile: {my_shapefile}",
file=sys.stderr,
)
d_init_slice = None d_init_slice = None
self.d_init.append(d_init_slice) self.d_init.append(d_init_slice)
self._n_dates = np.zeros(self.n_slices, dtype = int) self._n_dates = np.zeros(self.n_slices, dtype=int)
self._ishape_last = [{"Anticyclones": None, "Cyclones": None} self._ishape_last = [
for i in range(self.n_slices)] {"Anticyclones": None, "Cyclones": None}
for i in range(self.n_slices)
]
# Number of eddies at each date: # Number of eddies at each date:
self._e_max = [{"Anticyclones": None, "Cyclones": None} self._e_max = [
for i in range(self.n_slices)] {"Anticyclones": None, "Cyclones": None}
for i in range(self.n_slices)
]
self._n_shapes = [{"Anticyclones": 0, "Cyclones": 0} self._n_shapes = [
for i in range(self.n_slices)] {"Anticyclones": 0, "Cyclones": 0} for i in range(self.n_slices)
]
def get_n_dates(self, i_slice): def get_n_dates(self, i_slice):
"""Return the number of dates in a given slice, including missing """Return the number of dates in a given slice, including missing
...@@ -94,20 +109,24 @@ class SHPC_class: ...@@ -94,20 +109,24 @@ class SHPC_class:
if self._n_dates[i_slice] == 0: if self._n_dates[i_slice] == 0:
# Find out the number of dates in the slice: # Find out the number of dates in the slice:
if self._ishape_last[i_slice]["Anticyclones"] is not None: if self._ishape_last[i_slice]["Anticyclones"] is not None:
self._n_dates[i_slice] \ self._n_dates[i_slice] = self._ishape_last[i_slice][
= self._ishape_last[i_slice]["Anticyclones"].size "Anticyclones"
].size
elif self._ishape_last[i_slice]["Cyclones"] is not None: elif self._ishape_last[i_slice]["Cyclones"] is not None:
self._n_dates[i_slice] \ self._n_dates[i_slice] = self._ishape_last[i_slice][
= self._ishape_last[i_slice]["Cyclones"].size "Cyclones"
].size
else: else:
try: try:
self._n_dates[i_slice] \ self._n_dates[i_slice] = self.get_ishape_last(
= self.get_ishape_last(i_slice, self.def_orient).size i_slice, self.def_orient
).size
except AttributeError: except AttributeError:
# Try the other orientation: # Try the other orientation:
orientation = other_orient(self.def_orient) orientation = other_orient(self.def_orient)
self._n_dates[i_slice] \ self._n_dates[i_slice] = self.get_ishape_last(
= self.get_ishape_last(i_slice, orientation).size i_slice, orientation
).size
return self._n_dates[i_slice] return self._n_dates[i_slice]
...@@ -119,13 +138,15 @@ class SHPC_class: ...@@ -119,13 +138,15 @@ class SHPC_class:
""" """
if self._ishape_last[i_slice][orientation] is None: if self._ishape_last[i_slice][orientation] is None:
fname = path.join(self.dir, f"Slice_{i_slice}", orientation, fname = path.join(
"ishape_last.txt") self.dir, f"Slice_{i_slice}", orientation, "ishape_last.txt"
)
print("Reading", fname, "...") print("Reading", fname, "...")
try: try:
self._ishape_last[i_slice][orientation] \ self._ishape_last[i_slice][orientation] = np.loadtxt(
= np.loadtxt(fname, dtype = int, ndmin = 1) fname, dtype=int, ndmin=1
)
except OSError: except OSError:
print("Could not read", fname) print("Could not read", fname)
print("Creating it...") print("Creating it...")
...@@ -133,31 +154,32 @@ class SHPC_class: ...@@ -133,31 +154,32 @@ class SHPC_class:
if self._ishape_last[i_slice][orientation] is not None: if self._ishape_last[i_slice][orientation] is not None:
assert self._ishape_last[i_slice][orientation].size != 0 assert self._ishape_last[i_slice][orientation].size != 0
self._e_max[i_slice][orientation] \ self._e_max[i_slice][orientation] = np.ediff1d(
= np.ediff1d(self._ishape_last[i_slice][orientation], self._ishape_last[i_slice][orientation],
to_begin = self._ishape_last[i_slice]\ to_begin=self._ishape_last[i_slice][orientation][0] + 1,
[orientation][0] + 1) )
return self._ishape_last[i_slice][orientation] return self._ishape_last[i_slice][orientation]
def get_reader(self, i_slice, orientation, layer): def get_reader(self, i_slice, orientation, layer):
if layer not in self._readers[i_slice][orientation]: if layer not in self._readers[i_slice][orientation]:
my_shapefile = path.join(self.dir, f"Slice_{i_slice}", orientation, my_shapefile = path.join(
layer) self.dir, f"Slice_{i_slice}", orientation, layer
)
try: try:
self._readers[i_slice][orientation][layer] \ self._readers[i_slice][orientation][layer] = shapefile.Reader(
= shapefile.Reader(my_shapefile) my_shapefile
)
except shapefile.ShapefileException: except shapefile.ShapefileException:
print(f"get_reader: {my_shapefile} not found", print(f"get_reader: {my_shapefile} not found", file=sys.stderr)
file = sys.stderr)
self._readers[i_slice][orientation][layer] = None self._readers[i_slice][orientation][layer] = None
return self._readers[i_slice][orientation][layer] return self._readers[i_slice][orientation][layer]
def get_n_shapes(self, i_slice, orientation): def get_n_shapes(self, i_slice, orientation):
if self._n_shapes[i_slice][orientation] == 0: if self._n_shapes[i_slice][orientation] == 0:
reader = self.get_reader(i_slice, orientation, layer = "extremum") reader = self.get_reader(i_slice, orientation, layer="extremum")
self._n_shapes[i_slice][orientation] = len(reader) self._n_shapes[i_slice][orientation] = len(reader)
return self._n_shapes[i_slice][orientation] return self._n_shapes[i_slice][orientation]
...@@ -170,13 +192,13 @@ class SHPC_class: ...@@ -170,13 +192,13 @@ class SHPC_class:
""" """
last_date = self.d_init[i_slice] last_date = self.d_init[i_slice]
date_ishape = [] # list of tuples (date, ishape) date_ishape = [] # list of tuples (date, ishape)
reader = self.get_reader(i_slice, orientation, "extremum") reader = self.get_reader(i_slice, orientation, "extremum")
if reader is None: if reader is None:
print('Could not create "ishape_last.txt".') print('Could not create "ishape_last.txt".')
else: else:
for ishape, rec in enumerate(reader.iterRecords(fields = ["date"])): for ishape, rec in enumerate(reader.iterRecords(fields=["date"])):
if rec.date != last_date: if rec.date != last_date:
date_ishape.append((last_date, ishape - 1)) date_ishape.append((last_date, ishape - 1))
last_date = rec.date last_date = rec.date
...@@ -185,8 +207,8 @@ class SHPC_class: ...@@ -185,8 +207,8 @@ class SHPC_class:
date_ishape = np.array(date_ishape) date_ishape = np.array(date_ishape)
# Replace date by repetition of date factor: # Replace date by repetition of date factor:
date_ishape[:- 1, 0] = np.ediff1d(date_ishape[:, 0]) date_ishape[:-1, 0] = np.ediff1d(date_ishape[:, 0])
date_ishape[- 1, 0] = 1 date_ishape[-1, 0] = 1
ishape_last = [] ishape_last = []
...@@ -194,17 +216,19 @@ class SHPC_class: ...@@ -194,17 +216,19 @@ class SHPC_class:
ishape_last.extend(repet_factor * [ishape]) ishape_last.extend(repet_factor * [ishape])
self._ishape_last[i_slice][orientation] = np.array(ishape_last) self._ishape_last[i_slice][orientation] = np.array(ishape_last)
fname = path.join(self.dir, f"Slice_{i_slice}", orientation, fname = path.join(
"ishape_last.txt") self.dir, f"Slice_{i_slice}", orientation, "ishape_last.txt"
np.savetxt(fname, self._ishape_last[i_slice][orientation], )
fmt = "%.0d") np.savetxt(
fname, self._ishape_last[i_slice][orientation], fmt="%.0d"
)
def get_slice(self, date_index): def get_slice(self, date_index):
i_slice = bisect.bisect(self.d_init, date_index) - 1 i_slice = bisect.bisect(self.d_init, date_index) - 1
assert i_slice >= 0 assert i_slice >= 0
return i_slice return i_slice
def comp_ishape(self, date, eddy_index, orientation, i_slice = None): def comp_ishape(self, date, eddy_index, orientation, i_slice=None):
"""Compute the location in the shapefiles: return `(i_slice, ishape)`. """Compute the location in the shapefiles: return `(i_slice, ishape)`.
Compute i_slice if i_slice is None. Compute i_slice if i_slice is None.
...@@ -215,7 +239,9 @@ class SHPC_class: ...@@ -215,7 +239,9 @@ class SHPC_class:
""" """
if i_slice is None: i_slice = self.get_slice(date) if i_slice is None:
i_slice = self.get_slice(date)
ishape_last = self.get_ishape_last(i_slice, orientation) ishape_last = self.get_ishape_last(i_slice, orientation)
if ishape_last is None: if ishape_last is None:
...@@ -247,10 +273,15 @@ class SHPC_class: ...@@ -247,10 +273,15 @@ class SHPC_class:
if ishape_start is None: if ishape_start is None:
return [] return []
else: else:
ishape_stop = self._ishape_last[i_slice][orientation]\ ishape_stop = (
[date - self.d_init[i_slice]] + 1 self._ishape_last[i_slice][orientation][
date - self.d_init[i_slice]
]
+ 1
)
return range(ishape_start, ishape_stop) return range(ishape_start, ishape_stop)
def in_window(point, window): def in_window(point, window):
llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = window llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = window
longit = point[0] + math.ceil((llcrnrlon - point[0]) / 360) * 360 longit = point[0] + math.ceil((llcrnrlon - point[0]) / 360) * 360
...@@ -258,13 +289,15 @@ def in_window(point, window): ...@@ -258,13 +289,15 @@ def in_window(point, window):
return longit <= urcrnrlon and llcrnrlat <= point[1] <= urcrnrlat return longit <= urcrnrlon and llcrnrlat <= point[1] <= urcrnrlat
def other_orient(orientation): def other_orient(orientation):
if orientation == "Cyclones": if orientation == "Cyclones":
return "Anticyclones" return "Anticyclones"
else: else:
return "Cyclones" return "Cyclones"
def node_to_date_eddy(node_index, e_overestim, only_date = False):
def node_to_date_eddy(node_index, e_overestim, only_date=False):
"""Convert from node identification in graph to eddy identification in """Convert from node identification in graph to eddy identification in
shapefiles. Return date_index or (date_index, eddy_index). shapefiles. Return date_index or (date_index, eddy_index).
...@@ -276,6 +309,7 @@ def node_to_date_eddy(node_index, e_overestim, only_date = False): ...@@ -276,6 +309,7 @@ def node_to_date_eddy(node_index, e_overestim, only_date = False):
else: else:
return k, node_index - k * e_overestim return k, node_index - k * e_overestim
def date_eddy_to_node(date, eddy_index, e_overestim): def date_eddy_to_node(date, eddy_index, e_overestim):
"""Convert from eddy identification in shapefiles to node """Convert from eddy identification in shapefiles to node
identification in graph. identification in graph.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment