class BathyMap():
    """
    Create a bathymetry map figure and axis
    """
    def __init__(self, map_extent, bathy_file,
                 grid_x=None, grid_y=None, intens_file=None):
        """
        Set up a bathymetric map

        :param map_extent [list]: [minlon,maxlon,minlat,maxlat]
        :param bathy_file [str]: name of netcdf bathymetry file
        :param intens_file: name of netcdf intensity file
        :kind intens_file: str, optional
        :param grid_x: x grid spacing for the plot axis
        :param grid_y: y grid spacing for the plot axis
        """
        self.map_extent = map_extent
        # self.lon, self.lat, self.z = self._setup_bathy_vars()
        self.lon, self.lat, self.z = self._read_bathy_map(bathy_file)

        # set up figure and axes
        self.fig = plt.figure()
        self.ax = plt.axes(projection=ccrs.Mercator())
        self.ax.set_extent(self.map_extent)
        if grid_x:
            min_x = grid_x * np.floor(min(self.lon) / grid_x)
            grid_x = np.arange(min_x, max(self.lon), grid_x)
        if grid_y:
            min_y = grid_y * np.floor(min(self.lat) / grid_y)
            grid_y = np.arange(min_y, max(self.lat), grid_y)
        gl = self.ax.gridlines(xlocs=grid_x, ylocs=grid_y, draw_labels=True)
        gl.top_labels = False
        gl.right_labels = False

    def _read_bathy_map(self, fname):
        """
        Read netcdf grid file

        :param fname: filename
        """
        if not isinstance(fname, str):
            if fname is True:
                print('Sorry, ETOPO isn"t working')
            return self.map_extent[:2], self.map_extent[2:4], None

            # etopo = xr.open_dataset('etopo1.nc')
            # # Slice desired coordinates.
            # etopo = etopo.sel(x=slice(self.map_extent[0],self.map_extent[1]),
            #                   y=slice(self.map_extent[2],self.map_extent[3]))  
            # # Return sliced variables.
            # return etopo.x, etopo.y, etopo.z
        # Open ETOPO1 file and slice desired area.

        try:
            f = netcdf.netcdf_file(fname, 'r')
            if 'x' in f.variables:
                lon = f.variables['x'][:].copy()
                lat = f.variables['y'][:].copy()
                z = f.variables['z'][:].copy()
            elif 'x_range' in f.variables:
                x_spacing, y_spacing = f.variables['spacing'].data
                lon = np.arange(f.variables['x_range'].data[0],
                                f.variables['x_range'].data[1] + x_spacing / 2,
                                x_spacing)
                lat = np.arange(f.variables['y_range'][0],
                                f.variables['y_range'][1] + y_spacing / 2,
                                y_spacing)
                z = np.flipud(np.array(f.variables['z'].data).reshape(len(lat),
                                                                      len(lon)))
            f.close()
        except Exception as e:
            print(e)
            return None, None, None

        # Cut map down to desired range
        lon_step= lon[1] - lon[0]
        lat_step= lat[1] - lat[0]
        if self.map_extent[0] < min(lon):
            warnings.warn(f'Requested min lon < map minlon: increasing to {min(lon)}')
            self.map_extent[0] = min(lon) + lon_step/10
        if self.map_extent[1] > max(lon):
            warnings.warn(f'Requested max lon > map maxlon: decreasing to {max(lon)}')
            self.map_extent[1] = max(lon) - lon_step/10
        if self.map_extent[2] < min(lat):
            warnings.warn(f'Requested min lat < map minlat: increasing to {min(lat)}')
            self.map_extent[2] = min(lat) + lat_step/10
        if self.map_extent[3] > max(lat):
            warnings.warn(f'Requested max lat > map maxlat: decreasing to {max(lat)}')
            self.map_extent[3] = max(lat) - lat_step/10
        ixmin = np.flatnonzero((lon >= self.map_extent[0])).min()
        ixmax = np.flatnonzero((lon <= self.map_extent[1])).max()
        iymin = np.flatnonzero((lat >= self.map_extent[2])).min()
        iymax = np.flatnonzero((lat <= self.map_extent[3])).max()
        return lon[ixmin:ixmax], lat[iymin:iymax], z[iymin:iymax, ixmin:ixmax]

    def plot_image(self, cmap=cm.jet, pastel=False, shaded=False, intens_file=None, **kwargs):
        """
        Plot the bathymetric image

        :param cmap: colormap (e.g matplotlib.colors.LinearSegmentedColormap
              instance)
        :param pastel: lighten colors to make pastel
        :param shaded: [bool] apply shading
        :param intens_file: file containing intensities to plot (generally for shading)
        :param **kwargs: arguments to pass to hillshade (scale, azdeg, altdeg)
        """
        intens = None
        if intens_file is not None:
            lon, lat, intens = self._read_bathy_map(intens_file)
            if not (len(lon) == len(self.lon)) or not (len(lat) == len(self.lat)):
                warnings.warn('intensity lon/lat values different length than in bathy map')
            elif not all(abs(lon - self.lon) < ((lon[1] - lon[0]) / 10)) \
                    or not all(abs(lat - self.lat) < ((lat[1] - lat[0]) / 10)):
                warnings.warn('intensity lon/lat values do not match bathy map')
        if intens is not None:
            z_shade = _set_shade(self.z, intensity=intens, cmap=cmap, **kwargs)
        else:
            z_shade = _set_shade(self.z, cmap=cmap, shaded=shaded, **kwargs)
        if pastel:
            z_shade = 0.5 * z_shade + 0.5

        self.ax.imshow(z_shade, origin='lower', extent=self.map_extent,
                       transform=ccrs.PlateCarree())

    def plot_contours(self, levels=500, linewidth=1, color='k'):
        """
        Plot the bathymetric contours

        :param levels: list of contours, or contour interval (m)
        :param linewidth: contour linewidth (1)
        :param colors: contour line color ('k')
        """
        if not isinstance(levels, list):
            interval = levels
            min_level = interval * np.floor(np.nanmin(self.z)/interval)
            levels = np.arange(min_level, np.nanmax(self.z), interval)
        plt.contour(self.lon, self.lat, self.z,
                    levels, colors=color,
                    linestyles='solid', linewidths=linewidth,
                    transform=ccrs.PlateCarree())

    def show(self):
        """
        Show the plot on the screen
        """
        plt.show()

    def plot_coastlines(self, resolution):
        """
        Plot coastlines

        :param resolution: what resolution coastlines to include.  Must
            correspond to a NaturalEarth resolution ('10m', '50m', '110m')
            or a GSHSS resolution ('auto', 'low', 'high', 'full'...)
        """
        NE_coast_resolutions = ['10m', '50m', '110m']
        GSHSS_coast_resolutions = ['auto', 'coarse', 'low',
                                   'intermediate', 'high',
                                   'full']
        if resolution in GSHSS_coast_resolutions:
            coast = GSHHSFeature(scale=resolution)
            self.ax.add_feature(coast)
        elif resolution in NE_coast_resolutions:
            coast = NaturalEarthFeature(scale=resolution)
            self.ax.add_feature(coast)
            # self.ax.coastlines(resolution=resolution)
        else:
            print(f'Invalid coastline resolution: "{resolution}"')


    def plot(self, lon, lat, sym='o', color='deepskyblue',
                     mec='k', ms=8, **kwargs):
        """
        plot a point

        :param lon:  longitude
        :param lat:  latitude
        :param sym:  symbol
        :param color: marker face color
        :param kwargs: keyword arguments to pass on to pyplot.plot()
        """
        self.ax.plot(lon, lat, sym, color=color,
                     transform=ccrs.PlateCarree(),
                     **kwargs)

    def arrow(self, x, y, dx, dy, **kwargs):
        """
        plot a arrow

        :param x:  longitude
        :param y:  latitude
        :param dx:  arrow length
        :param dy:  arrow length
        :param kwargs: keyword arguments to pass on to pyplot.arrow()
        """
        self.ax.arrow(x, y, dx, dy, transform=ccrs.PlateCarree(), **kwargs)

    def text(self, lon, lat, name, text_box=True,
                   ha='left', va='center', color='black', **kwargs):
        """
        plot a station label

        :param m: map object
        :param lon: longitude
        :param lat: latitude
        :param name: station name
        :param text_box: Put a semi-opaque box under the text
        :param ha: horizontal text alignment ('left', 'center', 'right')
        :param va: vertical text alignment ('center', 'middle', 'bottom')
        :param color: text color
        :param kwargs: keyword arguments to pass on to pyplot.text()
        """
        if text_box:
            kwargs['bbox'] = dict(boxstyle="round", pad=0.1, ec='none', fc='white', alpha=0.5)
        self.ax.text(lon, lat, name, va=va, ha=ha,
                     color=color, transform=ccrs.PlateCarree(), **kwargs)

    def save_map(self, filename, title, fontsize=12):
        """
        Saves map to a file

        :param filename:
        :param title: plot title
        :param timedelta_total: total time spent (datetime.timedelta)
        :returns: figure object, base_name of output file
        """
        self.ax.set_title(title, fontsize=fontsize)
        self.fig.savefig(filename)

class MapCrossSection():
    """
    Create axes map with cross-sections to right (east) and below (south)

    the axis attributes to use are .ax_map, .ax_xz and .ax_zy
    .ax_map is a Mercator projection, when you plot to it you must
        specify transform=cartopy.crs.PlateCarree()
    """
    def __init__(self, figsize=None, fill_figure=True, pad_left=0.6,
                 pad_right=0.6, pad_top=0.4, pad_bottom=0.4,
                 map_extent=None, depth_extent=None, extent_buffer=0.2,
                 latitudes=None, longitudes=None, depths=None, gridlines=True,
                 coastlines=None):
        """
        :param figsize: width, height in inches
        :kind figsize: (float, float), optional
        :param fill_figure: Expand map bounds to fill the figure
        :param pad_left: inches to left of axes for labels/text/emptiness
        :param pad_right: inches to right of axes for labels/text/emptiness
        :param pad_top: inches above axes for labels/text/emptiness
        :param pad_bottom: inches below axes for labels/text/emptiness
        :param map_extent: bounds of map (x0, x1, y0, y1) or None
        :param depth_extent: depth bounds (z0, z1) or None
        :param latitudes: list of "event" latitudes, used to calculate map
            extents if map_extents == None
        :param longitudes: list of "event" longitudes, used to calculate map
            extents if map_extents == None
        :param longitudes: list of "event" depths, used to calculate depth
            extents if depth_extents == None
        :param extent_buffer: extents will exceed event limits by this
            fraction of the event range
        :param gridlines: Put gridlines on map view (boolean)
        :param coastlines: what resolution coastlines to include.  Must
            correspond to a NaturalEarth resolution ('10m', '50m', '110m')
            or a GSHSS resolution ('auto', 'low', 'high', 'full'...)
        """
        self.km_per_deg_lat = 1.852 * 60.
        self.map_extent = self._get_map_extent(map_extent, latitudes,
                                               longitudes, extent_buffer)
        self.depth_extent = self._get_depth_extent(depth_extent, depths,
                                                   extent_buffer)
        self.pad_left = pad_left
        self.pad_right = pad_right
        self.pad_top = pad_top
        self.pad_bottom = pad_bottom

        self.coastlines = None
        if coastlines is not None:
            if ((coastlines in NE_coast_resolutions)
                    or (coastlines in GSHSS_coast_resolutions)):
                self.coastlines = coastlines
            else:
                print(f'Invalid coastlines string: "{coastlines}"')

        self.fig = plt.figure(figsize)
        rect_map, rect_xz, rect_zy = self._make_rects(fill_figure)

        # Setup Lat-Z axis (lower left)
        self.ax_xz = plt.axes(rect_xz,
                              xlim=self.map_extent[:2],
                              ylim=self.depth_extent[::-1])
        self.ax_xz.tick_params(direction='in')

        # Setup Z-lon axis (upper right)
        self.ax_zy = plt.axes(rect_zy,
                              xlim=self.depth_extent,
                              ylim=self.map_extent[2:])
        self.ax_zy.tick_params(direction='in')
        self.ax_zy.xaxis.set_ticks_position('top')
        self.ax_zy.yaxis.set_ticks_position('right')

        # Setup map axis (upper left)
        self.ax_map = plt.axes(rect_map, projection=ccrs.Mercator())
        self.ax_map.set_extent(self.map_extent)
        if self.coastlines in GSHSS_coast_resolutions:
            coast = GSHHSFeature(scale=self.coastlines)
            self.ax_map.add_feature(coast)
            # feature = self.ax_map.add_feature(coast)
        elif self.coastlines in NE_coast_resolutions:
            self.ax_map.coastlines(resolution='10m')
        gl = self.ax_map.gridlines(draw_labels=True)
        gl.right_labels = False
        gl.bottom_labels = False
        if not gridlines:
            gl.xlines = False
            gl.ylines = False

    @property
    def lat_min(self):
        return self.map_extent[2]

    @property
    def lat_max(self):
        return self.map_extent[3]

    @property
    def lon_min(self):
        return self.map_extent[0]

    @property
    def lon_max(self):
        return self.map_extent[1]

    @property
    def depth_min(self):
        return self.depth_extent[0]

    @property
    def depth_max(self):
        return self.depth_extent[1]

    @property
    def km_per_deg_lon(self):
        return self.km_per_deg_lat\
            * np.cos(np.radians((self.lat_min + self.lat_max) / 2.))

    @property
    def lonrange(self):
        return self.lon_max - self.lon_min

    @property
    def latrange(self):
        return self.lat_max - self.lat_min

    @property
    def depthrange(self):
        return self.depth_extent[1] - self.depth_extent[0]

    @property
    def lonrange_km(self):
        return self.lonrange * self.km_per_deg_lon

    @property
    def latrange_km(self):
        return self.latrange * self.km_per_deg_lat

    @property
    def xrange_km(self):
        return self.lonrange_km + self.depthrange

    @property
    def yrange_km(self):
        return self.latrange_km + self.depthrange

    @property
    def plot_aspect(self):
        return self.xrange_km / self.yrange_km

    def _get_map_extent(self, map_extent, latitudes, longitudes,
                        extent_buffer):
        if map_extent:
            assert(len(map_extent) == 4),\
                f'len(map_extent) != 4 ({map_extent})'
            # print(map_extents)
            assert map_extent[1] > map_extent[0],\
                f'map_extent: x0({map_extent[0]:g}) >= x1({map_extent[1]:g})'
            assert map_extent[3] > map_extent[2], 'map_extent: y0 >= y1'
            lon_min, lon_max, lat_min, lat_max = map_extent
        else:
            assert latitudes is not None and longitudes is not None,\
                'need either extent or latitudes and longitudes'
            lon_min, lon_max = min(longitudes), max(longitudes)
            lat_min, lat_max = min(latitudes), max(latitudes)
            lon_buffer = (lon_max - lon_min) * extent_buffer
            lat_buffer = (lat_max - lat_min) * extent_buffer
            lon_min -= lon_buffer
            lon_max += lon_buffer
            lat_min -= lat_buffer
            lat_max += lat_buffer
            map_extent = (lon_min, lon_max, lat_min, lat_max)
        return map_extent

    def _get_depth_extent(self, depth_extent, depths, extent_buffer):
        if depth_extent:
            assert depth_extent[1] > depth_extent[0], 'depth_extent: z0 >= z1'
        else:
            assert depths is not None, 'need either depth_extents or depths'
            depth_min, depth_max = min(depths), max(depths)
            depth_buffer = (depth_max - depth_min) * extent_buffer
            depth_min -= depth_buffer
            depth_max += depth_buffer
            depth_extent = (depth_min, depth_max)
        return depth_extent

    def _adjustminmax(self, fig_aspect):
        if fig_aspect > self.plot_aspect:   # expand x_range
            add_kms = ((self.yrange_km * fig_aspect) - self.xrange_km) / 2.
            add_deg = add_kms / self.km_per_deg_lon
            self.map_extent[0] -= add_deg
            self.map_extent[1] += add_deg
        else:   # expand y_range
            add_kms = ((self.xrange_km / fig_aspect) - self.yrange_km) / 2.
            add_deg = add_kms / self.km_per_deg_lat
            self.map_extent[2] -= add_deg
            self.map_extent[3] += add_deg
        return

    def _make_figure(self, figsize):
        fig = plt.figure(figsize)
        # fig_inches = fig.get_size_inches()
        # fig_xrange = fig_inches[0] - self.pad_left - self.pad_right
        # fig_yrange = fig_inches[1] - self.pad_top - self.pad_bottom
        return fig

    def _make_rects(self, fill_figure):
        """
        Make rectangles for the three axes

        :param fill_figure: if true, extend map bounds to fill figure
        """
        fig_inches_x, fig_inches_y = self.fig.get_size_inches()
        fig_xrange = fig_inches_x - self.pad_left - self.pad_right
        fig_yrange = fig_inches_y - self.pad_top - self.pad_bottom
        fig_aspect = fig_xrange / fig_yrange

        if fill_figure:
            # Expand map extent to fill figure
            self._adjustminmax(fig_aspect)

        inch_per_km_x = fig_xrange / self.xrange_km
        inch_per_km_y = fig_yrange / self.yrange_km
        if self.plot_aspect > fig_aspect:
            inch_per_km_y *= fig_aspect / self.plot_aspect
        elif self.plot_aspect < fig_aspect:
            inch_per_km_x *= self.plot_aspect / fig_aspect

        # make axes rects
        left_width = self.lonrange_km * inch_per_km_x / fig_inches_x
        right_width = self.depthrange * inch_per_km_x / fig_inches_x
        lower_height = self.depthrange * inch_per_km_y / fig_inches_y
        upper_height = self.latrange_km * inch_per_km_y / fig_inches_y

        left_left = self.pad_left / fig_inches_x
        right_left = left_left + left_width
        lower_bottom = self.pad_bottom / fig_inches_y
        upper_bottom = lower_bottom + lower_height

        rect_map = [left_left, upper_bottom, left_width, upper_height]
        rect_xz = [left_left, lower_bottom, left_width, lower_height]
        rect_zy = [right_left, upper_bottom, right_width, upper_height]
        return rect_map, rect_xz, rect_zy

class Circle(ShapeSelection):
    """
    Selects events within a given distance of a given point
    """
    def __init__(self, longitude, latitude, radius):
        """
        :param latitude: latitude at the center of the circle
        :param longitude: longitude at the center of the circle
        :param radius: radius of the circle in degrees latitude
        :kind latitude, longitude, radius_km: num

        >>> Circle(-12.3, 45.6, 10)
        Circle(-12.3, 45.6, 10)
        """
        self.longitude = longitude
        self.latitude = latitude
        self.radius = radius

    def __repr__(self):
        return "Circle({:g}, {:g}, {:g})".format(
            self.longitude, self.latitude, self.radius)

    def contains(self, event_longitude, event_latitude):
        """
        Returns True if event is inside or on the Circle, False if not

        :param event_latitude: latitude of the event
        :param event_longitude: longitude of the event
        :type event_latitude, event_longitude: numeric

        >>> Circle(-12.3, 45.6, 10)._iswithin(-12.35, 45.65)
        True
        >>> Circle(-12.3, 45.6, 10)._iswithin(-12.4, 45.7)
        False
        """
        dist = gps2dist_azimuth(event_latitude, event_longitude,
                                self.latitude, self.longitude)[0]
        return dist/(1852*60) <= self.radius

    def plot(self, ax=None, **kwargs):
        if not ax:
            fig, ax = plt.subplots()
        degrees_long_per_lat = 1./np.cos(np.radians(self.latitude))
        thetas = np.radians(np.arange(0, 361))
        ax.plot(self.longitude
                + self.radius*np.cos(thetas) * degrees_long_per_lat,
                self.latitude + self.radius*np.sin(thetas))

    def add_patch(self, ax=None, **kwargs):
        patch = patches.Circle((self.longitude, self.latitude),
                               self.radius, **kwargs)
        if not ax:
            fig, ax = plt.subplots()
        ax.add_patch(patch)

class Polygon(ShapeSelection):
    """
    Selection of events within or outside of a polygon
    """
    def __init__(self, polygon):
        """
        :param polygon: list of [longitudes, latitudes]

        No need to close the polygon
        """
        self.polygon = polygon

    def contains(self, event_longitude, event_latitude):
        """
        Returns True if event is inside or on the Polygon, False if not

        :param event_longitude: longitude of the event
        :param event_latitude: latitude of the event
        :type event_latitude, event_longitude: numeric

        >>> Polygon([[12, 45], [12, 47], [14, 47], [14, 45]])._iswithin(13, 46)
        True
        >>> Polygon([[12, 45], [12, 47], [14, 47], [14, 45]])._iswithin(15, 46)
        False
        """
        path = paths.Path(self.polygon)
        return path.contains_point((event_longitude, event_latitude))

    def plot(self, ax=None, **kwargs):
        """Plots the shape's outline"""
        if not ax:
            fig, ax = plt.subplots()
        p = self.polygon
        if not p[0] == p[-1]:
            p.append(p[0])
        ax.plot([x[0] for x in p], [x[1] for x in p], **kwargs)

    def add_patch(self, ax=None, **kwargs):
        """Adds a patch corresponding to the shape"""
        path = paths.Path(self.polygon)
        if not ax:
            fig, ax = plt.subplots()
        patch = patches.PathPatch(path, lw=2, **kwargs)
        ax.add_patch(patch)

class ShapeSelection(object):
    """ Superclass for including or excluding events within the shape"""
    def include(self, catalog):
        """
        Return catalog with events within the shape

        :param catalog: earthquake catalog
        :kind catalog: ~class `obspy.core.event.catalog.Catalog`
        """
        new_catalog = Catalog()
        for event in catalog:
            origin = event.preferred_origin()
            if self.contains(origin.longitude, origin.latitude):
                new_catalog.append(event)
        return new_catalog

    def exclude(self, catalog):
        """
        Return catalog without any events in the shape

        :param catalog: earthquake catalog
        :kind catalog: ~class `obspy.core.event.catalog.Catalog`
        """
        new_catalog = Catalog()
        for event in catalog:
            origin = event.preferred_origin()
            if not self.contains(origin.longitude, origin.latitude):
                new_catalog.append(event)
        return new_catalog