-
Notifications
You must be signed in to change notification settings - Fork 95
/
cartopy.py
71 lines (62 loc) · 2.7 KB
/
cartopy.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (C) 2018-2021 Pyresample developers
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Classes for geometry operations.
See the related Cartopy pull request that implements this functionality
directly in Cartopy: https://github.com/SciTools/cartopy/pull/1888
"""
from __future__ import annotations
import numpy as np
import pyproj
try:
from xarray import DataArray
except ImportError:
DataArray = np.ndarray
import cartopy.crs as ccrs
class Projection(ccrs.Projection):
"""Flexible generic Projection with optionally specified bounds."""
def __init__(self,
crs: pyproj.CRS,
bounds: tuple[float, float, float, float] | None = None,
transform_bounds: bool = False):
"""Initialize CRS instance and compute bounds if possible."""
# NOTE: Skip the base cartopy Projection class __init__
super(ccrs.Projection, self).__init__(crs)
if bounds is None and self.area_of_use is not None:
bounds = (
self.area_of_use.west,
self.area_of_use.east,
self.area_of_use.south,
self.area_of_use.north,
)
transform_bounds = True
self.bounds = bounds
if bounds is not None and transform_bounds:
# Convert lat/lon bounds to projected bounds.
# Geographic area of the entire dataset referenced to WGS 84
# NB. We can't use a polygon transform at this stage because
# that relies on the existence of the map boundary... the very
# thing we're trying to work out! ;-)
x0, x1, y0, y1 = bounds
lons = np.array([x0, x0, x1, x1])
lats = np.array([y0, y1, y1, y0])
points = self.transform_points(self.as_geodetic(), lons, lats)
x = points[:, 0]
y = points[:, 1]
self.bounds = (x.min(), x.max(), y.min(), y.max())
if self.bounds is not None:
x0, x1, y0, y1 = self.bounds
self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.