Skip to content

Commit 16bcf54

Browse files
knutfrodeknutfrode
authored
Merge pull request #1634 from knutfrode/dev
Replaced @abstractproperty with @Property and @AbstractMethod
2 parents dbb4afa + 320bb47 commit 16bcf54

2 files changed

Lines changed: 111 additions & 112 deletions

File tree

‎opendrift/legacy.py‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,3 +41,105 @@ def gls_tke(windstress, depth, sea_water_density,
4141

4242
return K
4343

44+
#@require_mode(mode=Mode.Ready)
45+
def seed_from_shapefile_old(self,
46+
shapefile,
47+
number,
48+
layername=None,
49+
featurenum=None,
50+
**kwargs):
51+
"""Seeds elements within contours read from a shapefile
52+
53+
Obsolete, as new method based on geopandas is simpler"""
54+
55+
try:
56+
from osgeo import ogr, osr
57+
except Exception as e:
58+
logger.warning(e)
59+
raise ValueError('OGR library is needed to read shapefiles.')
60+
61+
if 'timeformat' in kwargs:
62+
# Recondstructing time from filename, where 'timeformat'
63+
# is forwarded to datetime.strptime()
64+
kwargs['time'] = datetime.strptime(os.path.basename(shapefile),
65+
kwargs['timeformat'])
66+
del kwargs['timeformat']
67+
68+
num_seeded_before = self.num_elements_scheduled()
69+
70+
targetSRS = osr.SpatialReference()
71+
targetSRS.ImportFromEPSG(4326)
72+
try:
73+
s = ogr.Open(shapefile)
74+
except:
75+
s = shapefile
76+
77+
for layer in s:
78+
if layername is not None and layer.GetName() != layername:
79+
logger.info('Skipping layer: ' + layer.GetName())
80+
continue
81+
else:
82+
logger.info('Seeding for layer: %s (%s features)' %
83+
(layer.GetDescription(), layer.GetFeatureCount()))
84+
85+
coordTrans = osr.CoordinateTransformation(layer.GetSpatialRef(),
86+
targetSRS)
87+
88+
if featurenum is None:
89+
featurenum = range(1, layer.GetFeatureCount() + 1)
90+
else:
91+
featurenum = np.atleast_1d(featurenum)
92+
if max(featurenum) > layer.GetFeatureCount():
93+
raise ValueError('Only %s features in layer.' %
94+
layer.GetFeatureCount())
95+
96+
# Loop first through all features to determine total area
97+
layer.ResetReading()
98+
area_srs = osr.SpatialReference()
99+
area_srs.ImportFromEPSG(3857)
100+
areaTransform = osr.CoordinateTransformation(
101+
layer.GetSpatialRef(), area_srs)
102+
103+
areas = np.zeros(len(featurenum))
104+
for i, f in enumerate(featurenum):
105+
feature = layer.GetFeature(f - 1) # Note 1-indexing, not 0
106+
if feature is not None:
107+
gom = feature.GetGeometryRef().Clone()
108+
gom.Transform(areaTransform)
109+
areas[i] = gom.GetArea()
110+
111+
total_area = np.sum(areas)
112+
layer.ResetReading() # Rewind to first layer
113+
logger.info('Total area of all polygons: %s m2' % total_area)
114+
# Find number of points per polygon
115+
numbers = np.round(number * areas / total_area).astype(int)
116+
numbers[numbers.argmax()] += int(number - sum(numbers))
117+
118+
for i, f in enumerate(featurenum):
119+
feature = layer.GetFeature(f - 1)
120+
if feature is None:
121+
continue
122+
num_elements = numbers[i]
123+
geom = feature.GetGeometryRef()
124+
logger.info(f'\tSeeding {num_elements} elements within polygon number {featurenum[i]} of area {areas[i]} m3')
125+
try:
126+
geom.Transform(coordTrans)
127+
except Exception as e:
128+
logger.warning('Could not transform coordinates:')
129+
logger.warning(e)
130+
pass
131+
#b = geom.GetBoundary()
132+
#if b is not None:
133+
# points = b.GetPoints()
134+
# lons = [p[0] for p in points]
135+
# lats = [p[1] for p in points]
136+
#else:
137+
# Alternative if OGR is not built with GEOS support
138+
r = geom.GetGeometryRef(0)
139+
lons = [r.GetY(j) for j in range(r.GetPointCount())]
140+
lats = [r.GetX(j) for j in range(r.GetPointCount())]
141+
142+
self.seed_within_polygon(lons=lons,
143+
lats=lats,
144+
number=num_elements,
145+
**kwargs)

‎opendrift/models/basemodel/__init__.py‎

Lines changed: 9 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@
3333
from opendrift.readers import reader_global_landmask
3434

3535
from datetime import datetime, timedelta
36-
from abc import ABCMeta, abstractmethod, abstractproperty
36+
from abc import ABCMeta, abstractmethod
3737

3838
import geojson
3939
import xarray as xr
@@ -677,10 +677,8 @@ def interact_with_coastline(self, final=False):
677677
logger.debug('%s elements hit coastline, '
678678
'moving back to water' % len(on_land))
679679
on_land_ID = self.elements.ID[on_land]
680-
self.elements.lon[on_land] = \
681-
self._elements_previous.lon[on_land_ID]
682-
self.elements.lat[on_land] = \
683-
self._elements_previous.lat[on_land_ID]
680+
self.elements.lon[on_land] = self._elements_previous.lon[on_land_ID]
681+
self.elements.lat[on_land] = self._elements_previous.lat[on_land_ID]
684682
self.environment.land_binary_mask[on_land] = 0
685683

686684
def interact_with_seafloor(self):
@@ -717,8 +715,8 @@ def interact_with_seafloor(self):
717715
logger.warning('%s elements hit seafloor, '
718716
'moving back ' % len(below))
719717
below_ID = self.elements.ID[below]
720-
self.elements.lon[below] = self._elements_previous.lon[below_ID].data
721-
self.elements.lat[below] = self._elements_previous.lat[below_ID].data
718+
self.elements.lon[below] = self._elements_previous.lon[below_ID]
719+
self.elements.lat[below] = self._elements_previous.lat[below_ID]
722720

723721
@abstractmethod
724722
def update(self):
@@ -727,11 +725,13 @@ def update(self):
727725
update properties (including position) of its particles (self.elements)
728726
"""
729727

730-
@abstractproperty
728+
@property
729+
@abstractmethod
731730
def ElementType(self):
732731
"""Any trajectory model implementation must define an ElementType."""
733732

734-
@abstractproperty
733+
@property
734+
@abstractmethod
735735
def required_variables(self):
736736
"""Any trajectory model implementation must list needed variables."""
737737

@@ -1524,109 +1524,6 @@ def seed_from_geopandas(self,
15241524
time=time,
15251525
**kwargs)
15261526

1527-
# @require_mode(mode=Mode.Ready)
1528-
# def seed_from_shapefile_old(self,
1529-
# shapefile,
1530-
# number,
1531-
# layername=None,
1532-
# featurenum=None,
1533-
# **kwargs):
1534-
# """Seeds elements within contours read from a shapefile
1535-
1536-
# Obsolete, as new method based on geopandas is simpler"""
1537-
1538-
# try:
1539-
# from osgeo import ogr, osr
1540-
# except Exception as e:
1541-
# logger.warning(e)
1542-
# raise ValueError('OGR library is needed to read shapefiles.')
1543-
1544-
# if 'timeformat' in kwargs:
1545-
# # Recondstructing time from filename, where 'timeformat'
1546-
# # is forwarded to datetime.strptime()
1547-
# kwargs['time'] = datetime.strptime(os.path.basename(shapefile),
1548-
# kwargs['timeformat'])
1549-
# del kwargs['timeformat']
1550-
1551-
# num_seeded_before = self.num_elements_scheduled()
1552-
1553-
# targetSRS = osr.SpatialReference()
1554-
# targetSRS.ImportFromEPSG(4326)
1555-
# try:
1556-
# s = ogr.Open(shapefile)
1557-
# except:
1558-
# s = shapefile
1559-
1560-
# for layer in s:
1561-
# if layername is not None and layer.GetName() != layername:
1562-
# logger.info('Skipping layer: ' + layer.GetName())
1563-
# continue
1564-
# else:
1565-
# logger.info('Seeding for layer: %s (%s features)' %
1566-
# (layer.GetDescription(), layer.GetFeatureCount()))
1567-
1568-
# coordTrans = osr.CoordinateTransformation(layer.GetSpatialRef(),
1569-
# targetSRS)
1570-
1571-
# if featurenum is None:
1572-
# featurenum = range(1, layer.GetFeatureCount() + 1)
1573-
# else:
1574-
# featurenum = np.atleast_1d(featurenum)
1575-
# if max(featurenum) > layer.GetFeatureCount():
1576-
# raise ValueError('Only %s features in layer.' %
1577-
# layer.GetFeatureCount())
1578-
1579-
# # Loop first through all features to determine total area
1580-
# layer.ResetReading()
1581-
# area_srs = osr.SpatialReference()
1582-
# area_srs.ImportFromEPSG(3857)
1583-
# areaTransform = osr.CoordinateTransformation(
1584-
# layer.GetSpatialRef(), area_srs)
1585-
1586-
# areas = np.zeros(len(featurenum))
1587-
# for i, f in enumerate(featurenum):
1588-
# feature = layer.GetFeature(f - 1) # Note 1-indexing, not 0
1589-
# if feature is not None:
1590-
# gom = feature.GetGeometryRef().Clone()
1591-
# gom.Transform(areaTransform)
1592-
# areas[i] = gom.GetArea()
1593-
1594-
# total_area = np.sum(areas)
1595-
# layer.ResetReading() # Rewind to first layer
1596-
# logger.info('Total area of all polygons: %s m2' % total_area)
1597-
# # Find number of points per polygon
1598-
# numbers = np.round(number * areas / total_area).astype(int)
1599-
# numbers[numbers.argmax()] += int(number - sum(numbers))
1600-
1601-
# for i, f in enumerate(featurenum):
1602-
# feature = layer.GetFeature(f - 1)
1603-
# if feature is None:
1604-
# continue
1605-
# num_elements = numbers[i]
1606-
# geom = feature.GetGeometryRef()
1607-
# logger.info(f'\tSeeding {num_elements} elements within polygon number {featurenum[i]} of area {areas[i]} m3')
1608-
# try:
1609-
# geom.Transform(coordTrans)
1610-
# except Exception as e:
1611-
# logger.warning('Could not transform coordinates:')
1612-
# logger.warning(e)
1613-
# pass
1614-
# #b = geom.GetBoundary()
1615-
# #if b is not None:
1616-
# # points = b.GetPoints()
1617-
# # lons = [p[0] for p in points]
1618-
# # lats = [p[1] for p in points]
1619-
# #else:
1620-
# # Alternative if OGR is not built with GEOS support
1621-
# r = geom.GetGeometryRef(0)
1622-
# lons = [r.GetY(j) for j in range(r.GetPointCount())]
1623-
# lats = [r.GetX(j) for j in range(r.GetPointCount())]
1624-
1625-
# self.seed_within_polygon(lons=lons,
1626-
# lats=lats,
1627-
# number=num_elements,
1628-
# **kwargs)
1629-
16301527
@require_mode(mode=Mode.Ready)
16311528
def seed_letters(self, text, lon, lat, time, number, scale=1.2):
16321529
"""Seed elements within text polygons"""

0 commit comments

Comments
 (0)