Source code for opendrift.readers.reader_shape

# This file is part of OpenDrift.
#
# OpenDrift 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, version 2
#
# OpenDrift 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 OpenDrift.  If not, see <https://www.gnu.org/licenses/>.
#
# Copyright 2020, Gaute Hope, MET Norway

from opendrift.readers.basereader import BaseReader, ContinuousReader

import pyproj
import shapely
import shapely.ops
import shapely.vectorized
import cartopy
import itertools
import logging
logger = logging.getLogger(__name__)

[docs] class Reader(BaseReader, ContinuousReader): """ The shape reader can be used to load generic shapes as the 'landmask' variable. Args: :param shapes: shapely geometries. :type shapes: iterable. :param proj4_str: Proj.4 string of shape file projection coordinates (default: '+proj=lonlat +ellps=WGS84'). :type proj4_str: string. """ name = 'shape' variables = ['land_binary_mask'] proj4 = None crs = None polys = None land = None always_valid = True
[docs] @staticmethod def from_shpfiles(shpfiles, proj4_str = '+proj=lonlat +ellps=WGS84', invert=False): """ Construct a shape-reader from shape-files (.shp) Args: :param shapes: shape-file or files (.shp) :type shapes: string or list of file names as strings. :param proj4_str: Proj.4 string of shape file projection coordinates (default: '+proj=lonlat +ellps=WGS84'). :type proj4_str: string. """ if isinstance(shpfiles, list): shpfiles = shpfiles else: shpfiles = [ shpfiles ] # reading shapefiles shp_iters = [] for shp in shpfiles: logger.debug("Reading shapefile: %s" % shp) from cartopy import io reader = io.shapereader.Reader(shp) shp_iters.append(reader.geometries()) return Reader(itertools.chain(*shp_iters), proj4_str, invert=invert)
def __init__(self, shapes, proj4_str = '+proj=lonlat +ellps=WGS84', invert=False): self.invert = invert # True if polygons are lakes and not land areas self.proj4 = proj4_str self.crs = pyproj.CRS(self.proj4) super (Reader, self).__init__ () # Depth self.z = None # loading and pre-processing shapes self.polys = list(shapes) # reads geometries if Shape Readers assert len(self.polys) > 0, "no geometries loaded" logger.info("Pre-processing %d geometries" % len(self.polys)) self.land = shapely.ops.unary_union(self.polys) self.xmin, self.ymin, self.xmax, self.ymax = self.land.bounds self.xmin, self.ymin = self.lonlat2xy(self.xmin, self.ymin) self.xmax, self.ymax = self.lonlat2xy(self.xmax, self.ymax) self.xmin = self.xmin[0] self.xmax = self.xmax[0] self.ymin = self.ymin[0] self.ymax = self.ymax[0]
[docs] def __on_land__(self, x, y): if self.invert is False: return shapely.vectorized.contains(self.land, x, y) else: # Inverse if polygons are lakes and not land areas return 1 - shapely.vectorized.contains(self.land, x, y)
[docs] def get_variables(self, requestedVariables, time = None, x = None, y = None, z = None): """ Get binary mask of whether elements are on land or not. Args: x (deg[]): longitude (decimal degrees) y (deg[]): latitude (decimal degrees) ... x, y is given in reader local projection. Returns: Binary mask of point x, y on land. """ self.check_arguments(requestedVariables, time, x, y, z) return { 'x' : x, 'y' : y, 'land_binary_mask': self.__on_land__(x,y) }