aboutsummaryrefslogtreecommitdiffstats
path: root/list-observations
diff options
context:
space:
mode:
authorGuilhem Moulin <guilhem@fripost.org>2023-10-23 21:55:28 +0200
committerGuilhem Moulin <guilhem@fripost.org>2023-10-23 21:57:26 +0200
commit0070da1fc157bf8bf9f2703bea5505843c39bde6 (patch)
tree2beb4267f64b06214a27764eea372836ccc064a4 /list-observations
parent3b0f985835ec0a6b628da5cc4d7573ce1d36a467 (diff)
Add new script to list and filter observations from Artdatabanken.
Usage example: ./list-observations --redlisted --max-distance=250 --geometry=/tmp/nygruva.geojson /tmp/nygruva/fynd.gpkg
Diffstat (limited to 'list-observations')
-rwxr-xr-xlist-observations432
1 files changed, 432 insertions, 0 deletions
diff --git a/list-observations b/list-observations
new file mode 100755
index 0000000..be8961f
--- /dev/null
+++ b/list-observations
@@ -0,0 +1,432 @@
+#!/usr/bin/python3
+
+#----------------------------------------------------------------------
+# List and filter observations from Artdatabanken
+# Copyright © 2023 Guilhem Moulin <guilhem@fripost.org>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+#----------------------------------------------------------------------
+
+import argparse
+from datetime import datetime
+import time
+import locale
+from pathlib import Path
+from osgeo import gdal, ogr, osr
+gdal.UseExceptions()
+
+target_srs = osr.SpatialReference()
+target_srs.ImportFromEPSG(3006)
+target_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+
+locale.setlocale(locale.LC_TIME, 'sv_SE.UTF-8')
+
+programName = 'list-observations'
+
+def ePath(v):
+ return Path(v).expanduser()
+
+
+class geometryAction(argparse.Action):
+ def __init__(self, option_strings, dest, nargs=None, **kwargs):
+ super().__init__(option_strings, dest, **kwargs)
+ def __call__(self, parser, namespace, value, option_string=None):
+ items = getattr(namespace, self.dest, None)
+ if items is None:
+ items = [value]
+ else:
+ items = items.copy()
+ items.append(value)
+ setattr(namespace, self.dest, items)
+
+class redlistedAction(argparse.Action):
+ def __init__(self, option_strings, dest, nargs=None, **kwargs):
+ super().__init__(option_strings, dest, nargs=nargs, **kwargs)
+ def __call__(self, parser, namespace, values, option_string=None):
+ items = getattr(namespace, self.dest, None)
+ if values is None:
+ items = 'TaxonIsRedlisted'
+ else:
+ values = values.split(',')
+ for v in values:
+ if v not in ['EX', 'EW', 'RE', 'CR', 'EN', 'VU', 'NT', 'LC', 'DD', 'NE', 'NA']:
+ raise Exception(f'{v}: Invalid redlist category')
+ if items is None:
+ items = values
+ else:
+ items = items.copy()
+ items.extend(values)
+ setattr(namespace, self.dest, items)
+
+class taxonAction(argparse.Action):
+ def __init__(self, option_strings, dest, nargs=None, **kwargs):
+ super().__init__(option_strings, dest, nargs=nargs, **kwargs)
+ def __call__(self, parser, namespace, values, option_string=None):
+ items = getattr(namespace, self.dest, None)
+ values = values.lower().split(',')
+ if items is None:
+ items = values
+ else:
+ items = items.copy()
+ items.extend(values)
+ setattr(namespace, self.dest, items)
+
+parser = argparse.ArgumentParser(
+ description='List and filter observations from Artdatabanken.',
+ prog = programName,
+)
+parser.add_argument('--geometry', type=ePath, metavar='FILE', action=geometryAction,
+ help='Geometry file of interest')
+parser.add_argument('--max-uncertainty', metavar='N', type=int, default=250,
+ help='Maximum coordinate uncertainty in meters (default: %(default)s)')
+parser.add_argument('--max-distance', metavar='N', type=int,
+ help='Maximum distance in meters from target geometry')
+parser.add_argument('--since', metavar='YYYY-MM-DD',
+ help='Start date for observations in ISO 8601 format')
+parser.add_argument('--until', metavar='YYYY-MM-DD',
+ help='End date for observations in ISO 8601 format')
+parser.add_argument('--protected', action='store_true',
+ help='Only consider taxons that are protected by law')
+parser.add_argument('--redlisted', nargs='?', metavar='[NT,…]', action=redlistedAction,
+ help='Only consider redlisted taxons')
+parser.add_argument('--taxon', metavar='NAME', action=taxonAction,
+ help='Only consider these taxons')
+parser.add_argument('observations', type=ePath, metavar='OBSERVATIONS_FILE',
+ help='Observations file')
+
+args = parser.parse_args()
+
+for k in ['since', 'until']:
+ v = getattr(args, k, None)
+ if v is not None:
+ d = datetime.fromisoformat(v)
+ setattr(args, k, d)
+
+geometry = None
+for path in args.geometry:
+ path = path.as_posix()
+ ds = ogr.Open(path, update=0)
+
+ layer = None
+ for i in range(ds.GetLayerCount()):
+ lyr = ds.GetLayerByIndex(i)
+ if lyr.GetGeomType() in [ogr.wkbPolygon,
+ ogr.wkbMultiPolygon,
+ ogr.wkbPolygonM,
+ ogr.wkbMultiPolygonM,
+ ogr.wkbPolygonZM,
+ ogr.wkbMultiPolygonZM,
+ ogr.wkbPolygon25D,
+ ogr.wkbMultiPolygon25D]:
+ layer = lyr
+ break
+ if layer is None:
+ raise Exception(f'{path}: does not have any layer of geometry type Polygon')
+
+ source_srs = layer.GetSpatialRef()
+ transform_srs = osr.CoordinateTransformation(source_srs, target_srs)
+
+ feature = layer.GetNextFeature()
+ while feature is not None:
+ geom = feature.GetGeometryRef()
+ geom = geom.Clone()
+ geom.FlattenTo2D()
+ geom.Transform(transform_srs)
+
+ if geometry is None:
+ geometry = ogr.Geometry(ogr.wkbMultiPolygon)
+
+ geometryType = geom.GetGeometryType()
+ if geometryType == ogr.wkbPolygon:
+ geometry.AddGeometry(geom)
+ elif geometryType == ogr.wkbMultiPolygon:
+ for i in range(geom.GetGeometryCount()):
+ geom2 = geom.GetGeometryRef(i)
+ geometry.AddGeometry(geom2)
+ else:
+ raise Exception(f'Unsuported geometry {geometryType}')
+ feature = layer.GetNextFeature()
+ layer = None
+ ds = None
+
+if geometry is not None:
+ geometry = geometry.UnionCascaded()
+
+
+path = args.observations.as_posix()
+ds = ogr.Open(path, update=0)
+
+layer = None
+for i in range(ds.GetLayerCount()):
+ lyr = ds.GetLayerByIndex(i)
+ if lyr.GetGeomType() in [ogr.wkbPoint, ogr.wkbPointM, ogr.wkbMultiPointM, ogr.wkbPointZM, ogr.wkbMultiPointZM, ogr.wkbPoint25D]:
+ layer = lyr
+ break
+if layer is None:
+ raise Exception(f'{path}: does not have any layer of geometry type Point')
+
+source_srs = layer.GetSpatialRef()
+transform_srs = osr.CoordinateTransformation(source_srs, target_srs)
+
+def getField(feat, name):
+ idx = feature.GetFieldIndex(name)
+ if idx is not None and idx >= 0:
+ return feature.GetField(idx)
+ return None
+def getFieldAsDateTime(feat, name):
+ idx = feature.GetFieldIndex(name)
+ if idx is not None and idx >= 0:
+ v = feature.GetFieldAsDateTime(idx)
+ t = (v[0], v[1], v[2], v[3], v[4], int(v[5]), v[6], -1, -1)
+ return datetime.fromtimestamp(time.mktime(t))
+ return None
+
+sum_count = 0
+summary = {}
+feature = layer.GetNextFeature()
+while feature is not None:
+ v = feature.GetField('IsPositiveObservation')
+ if v is not None and not v:
+ feature = layer.GetNextFeature()
+ continue
+ v = feature.GetField('IsNaturalOccurrence')
+ if v is not None and not v:
+ feature = layer.GetNextFeature()
+ continue
+
+ uncertainty = feature.GetField('CoordinateUncertaintyInMeters')
+ if uncertainty is not None and uncertainty > args.max_uncertainty:
+ feature = layer.GetNextFeature()
+ continue
+
+ vernacularName = getField(feature, 'VernacularName')
+ if vernacularName is None:
+ raise Exception('Missing name')
+ if args.taxon is not None and vernacularName.lower() not in args.taxon:
+ feature = layer.GetNextFeature()
+ continue
+
+ protectedByLaw = feature.GetField('ProtectedByLaw')
+ if args.protected and (protectedByLaw is None or not protectedByLaw):
+ feature = layer.GetNextFeature()
+ continue
+
+ taxonIsRedlisted = feature.GetField('TaxonIsRedlisted')
+ redlistCategory = getField(feature, 'RedlistCategory')
+ if args.redlisted is not None:
+ if args.redlisted == 'TaxonIsRedlisted':
+ if taxonIsRedlisted is None or not taxonIsRedlisted:
+ feature = layer.GetNextFeature()
+ continue
+ elif type(args.redlisted) == list:
+ if redlistCategory is None or redlistCategory not in args.redlisted:
+ feature = layer.GetNextFeature()
+ continue
+
+ startDate = getFieldAsDateTime(feature, 'StartDate')
+ if args.until is not None and (startDate is None or args.until < startDate):
+ feature = layer.GetNextFeature()
+ continue
+
+ endDate = getFieldAsDateTime(feature, 'EndDate')
+ if args.since is not None and (endDate is None or endDate < args.since):
+ feature = layer.GetNextFeature()
+ continue
+
+ geom = feature.GetGeometryRef()
+ geom = geom.Clone()
+ geom.FlattenTo2D()
+ geom.Transform(transform_srs)
+
+ distance = None
+ if geometry is not None and not geometry.Contains(geom):
+ distance = geom.Distance(geometry)
+ if distance == -1:
+ raise Exception('distance')
+ max_distance = 0 if args.max_distance is None else args.max_distance
+ if uncertainty is not None:
+ max_distance += uncertainty
+ if max_distance < distance:
+ feature = layer.GetNextFeature()
+ continue
+
+ if sum_count > 0:
+ print('')
+
+ print(f'Artnamn: {vernacularName.title()}')
+
+ name = vernacularName
+ if taxonIsRedlisted is not None and taxonIsRedlisted and redlistCategory is not None:
+ name += f' ({redlistCategory})'
+ n = summary.get(name, 0)
+ summary[name] = n + 1
+
+ v = getField(feature, 'OrganismQuantityInt')
+ if v is None:
+ v = 'Noterad'
+ else:
+ v = str(v)
+ unit = getField(feature, 'OrganismQuantityUnit')
+ if unit is not None and unit != '':
+ v += ' ' + unit
+ print(f'Antal & enhet: {v}')
+
+ v = getField(feature, 'Sex')
+ if v is not None:
+ print(f'Kön: {v}')
+
+ v = getField(feature, 'LifeStage')
+ if v is not None:
+ print(f'Ålder/stadium: {v}')
+
+ v = getField(feature, 'Activity')
+ if v is not None:
+ print(f'Aktivitet: {v}')
+
+ v = getField(feature, 'DiscoveryMethod')
+ if v is not None:
+ print(f'Metod: {v}')
+
+ v = getField(feature, 'ReportedBy')
+ if ',' in v:
+ print(f'Rapportörer: {v}')
+ else:
+ print(f'Rapportör: {v}')
+ v = getField(feature, 'RecordedBy')
+ if ',' in v:
+ print(f'Observatörer: {v}')
+ else:
+ print(f'Observatör: {v}')
+
+ print(f'Startdatum: {startDate.strftime("%-d %B %Y %H:%M")}')
+ print(f'Slutdatum: {endDate.strftime("%-d %B %Y %H:%M")}')
+
+ v = getField(feature, 'Projects')
+ if v is not None:
+ print(f'Projects: {v}')
+
+ v = getField(feature, 'OccurrenceRemarks')
+ if v is not None:
+ print(f'Kommentar: {v}')
+
+ v = getField(feature, 'DatasetName')
+ if v is not None:
+ print(f'Källa: {v}')
+ v = getField(feature, 'Url')
+ if v is not None:
+ print(f'URL: {v}')
+
+ if geom.GetGeometryType() != ogr.wkbPoint:
+ raise Exception('Not a point')
+ x, y = geom.GetPoint_2D(i)
+ v = getField(feature, 'Locality')
+ if v is not None:
+ print(f'Lokalnamn: {v}')
+ coords = f'Koordinater: N{y:.0f}, Ö{x:.0f} (±{uncertainty:.0f}m) {target_srs.GetName()}'
+
+ if geometry is not None:
+ if distance is None:
+ dist = 'i'
+ else:
+ dist = f'{distance:.0f}m från'
+ coords += f'; {dist} målgeometrin'
+ print(coords)
+
+ if redlistCategory is None:
+ v = '—'
+ elif redlistCategory == 'EX':
+ v = 'Utdöd (EX)'
+ elif redlistCategory == 'EW':
+ v = 'Utdöd i vilt tillstånd (EW)'
+ elif redlistCategory == 'RE':
+ v = 'Nationellt utdöd (RE)'
+ elif redlistCategory == 'CR':
+ v = 'Akut hotad (CR)'
+ elif redlistCategory == 'EN':
+ v = 'Starkt hotad (EN)'
+ elif redlistCategory == 'VU':
+ v = 'Sårbar (VU)'
+ elif redlistCategory == 'NT':
+ v = 'Nära hotad (NT)'
+ elif redlistCategory == 'LC':
+ v = 'Livskraftig (LC)'
+ elif redlistCategory == 'DD':
+ v = 'Kunskapsbrist (DD)'
+ elif redlistCategory == 'NE':
+ v = 'Ej bedömd (NE)'
+ elif redlistCategory == 'NA':
+ v = 'Ej tillämplig (NA)'
+ else:
+ v = redlistCategory
+ print(f'Rödlistning: {v}')
+
+ v = []
+ if protectedByLaw is not None and protectedByLaw:
+ v.append('nationellt fridlyst')
+ v2 = getField(feature, 'TaxonIsInvasiveInSweden')
+ if v2 is not None and v2:
+ v.append('främmande i Sverige')
+ v2 = getField(feature, 'TaxonIsInvasiveEuRegulation')
+ if v2 is not None and v2:
+ v.append('invasiv enligt EU-förordning 1143/2014')
+ v2 = getField(feature, 'Natura2000HabitatsDirective')
+ if v2 is not None and v2:
+ v2 = []
+ v3 = getField(feature, 'Natura2000HabitatsDirectiveArticle2')
+ if v3 is not None and v3:
+ v3 = getField(feature, 'Natura2000HabitatsDirectiveArticle2PrioritySpecie')
+ if v3 is not None and v3:
+ v2.append('II (prioriterad art)')
+ else:
+ v2.append('II')
+ v3 = getField(feature, 'Natura2000HabitatsDirectiveArticle4')
+ if v3 is not None and v3:
+ v2.append('IV')
+ v3 = getField(feature, 'Natura2000HabitatsDirectiveArticle5')
+ if v3 is not None and v3:
+ v2.append('V')
+ if len(v2) == 0:
+ v.append('habitatdirektivet')
+ elif len(v2) == 1:
+ v.append('habitatdirektivets bilaga ' + v2[0])
+ else:
+ v.append('habitatdirektivets bilagor ' + ', '.join(v2[0:-1]) + ' & ' + v2[-1])
+ v2 = getField(feature, 'TaxonIsSignalSpecie')
+ if v2 is not None and v2:
+ v.append('signalart enligt Skogsstyrelsen')
+ v2 = getField(feature, 'BirdDirective')
+ if v2 is not None and v2:
+ v.append('fågeldirektivet')
+ if len(v) > 0:
+ v[0] = v[0][0].upper() + v[0][1:]
+ v = '; '.join(v)
+ else:
+ v = '—'
+ print(f'Naturvård: {v}')
+
+ sum_count += 1
+ feature = layer.GetNextFeature()
+
+if sum_count == 0:
+ print('No matching observations found!')
+ exit(0)
+
+max_count = max(summary.values())
+width = len(str(max_count))
+
+print('\n\nSummary\n-------')
+for name, n in sorted(summary.items(), key=lambda x: (-x[1], x[0])):
+ print('{0:{width}d}× {1}'.format(n, name, width=width))
+print('-' * (width+1))
+print(f'{sum_count} matching observations')