Skip to content

Commit

Permalink
Merge pull request #903 from mraspaud/fix-hrit-tests
Browse files Browse the repository at this point in the history
Fix HRV area definition tests
  • Loading branch information
mraspaud authored Sep 16, 2019
2 parents 0f76aea + c9f7f2f commit 1389229
Show file tree
Hide file tree
Showing 3 changed files with 25 additions and 23 deletions.
28 changes: 15 additions & 13 deletions satpy/readers/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,9 @@ def get_geostationary_angle_extent(geos_area):
# TODO: take into account sweep_axis_angle parameter

# get some projection parameters
req = geos_area.proj_dict['a'] / 1000
rp = geos_area.proj_dict['b'] / 1000
h = geos_area.proj_dict['h'] / 1000 + req
req = float(geos_area.proj_dict['a']) / 1000
rp = float(geos_area.proj_dict['b']) / 1000
h = float(geos_area.proj_dict['h']) / 1000 + req

# compute some constants
aeq = 1 - req**2 / (h ** 2)
Expand All @@ -77,14 +77,15 @@ def get_geostationary_angle_extent(geos_area):


def get_geostationary_mask(area):
"""Compute a mask of the earth's shape as seen by a geostationary satellite
"""Compute a mask of the earth's shape as seen by a geostationary satellite.
Args:
area (pyresample.geometry.AreaDefinition) : Corresponding area
definition
Returns:
Boolean mask, True inside the earth's shape, False outside.
"""
# Compute projection coordinates at the earth's limb
h = area.proj_dict['h']
Expand All @@ -101,12 +102,12 @@ def get_geostationary_mask(area):

def _lonlat_from_geos_angle(x, y, geos_area):
"""Get lons and lats from x, y in projection coordinates."""
h = (geos_area.proj_dict['h'] + geos_area.proj_dict['a']) / 1000
b__ = (geos_area.proj_dict['a'] / geos_area.proj_dict['b']) ** 2
h = float(geos_area.proj_dict['h'] + geos_area.proj_dict['a']) / 1000
b__ = (geos_area.proj_dict['a'] / float(geos_area.proj_dict['b'])) ** 2

sd = np.sqrt((h * np.cos(x) * np.cos(y)) ** 2 -
(np.cos(y)**2 + b__ * np.sin(y)**2) *
(h**2 - (geos_area.proj_dict['a'] / 1000)**2))
(h**2 - (float(geos_area.proj_dict['a']) / 1000)**2))
# sd = 0

sn = (h * np.cos(x) * np.cos(y) - sd) / (np.cos(y)**2 + b__ * np.sin(y)**2)
Expand All @@ -126,6 +127,7 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):
Args:
nb_points: Number of points on the polygon
"""
xmax, ymax = get_geostationary_angle_extent(geos_area)

Expand All @@ -136,7 +138,7 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):

# clip the projection coordinates to fit the area extent of geos_area
ll_x, ll_y, ur_x, ur_y = (np.array(geos_area.area_extent) /
geos_area.proj_dict['h'])
float(geos_area.proj_dict['h']))

x = np.clip(np.concatenate([x, x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
y = np.clip(np.concatenate([y, -y]), min(ll_y, ur_y), max(ll_y, ur_y))
Expand All @@ -146,7 +148,6 @@ def get_geostationary_bounding_box(geos_area, nb_points=50):

def get_area_slices(data_area, area_to_cover):
"""Compute the slice to read from an *area* based on an *area_to_cover*."""

if data_area.proj_dict['proj'] != 'geos':
raise NotImplementedError('Only geos supported')

Expand Down Expand Up @@ -191,8 +192,7 @@ def get_sub_area(area, xslice, yslice):


def unzip_file(filename):
"""Unzip the file if file is bzipped = ending with 'bz2'"""

"""Unzip the file if file is bzipped = ending with 'bz2'."""
if filename.endswith('bz2'):
bz2file = bz2.BZ2File(filename)
fdn, tmpfilepath = tempfile.mkstemp()
Expand All @@ -212,12 +212,13 @@ def unzip_file(filename):


def bbox(img):
"""Find the bounding box around nonzero elements in the given array
"""Find the bounding box around nonzero elements in the given array.
Copied from https://stackoverflow.com/a/31402351/5703449 .
Returns:
rowmin, rowmax, colmin, colmax
"""
rows = np.any(img, axis=1)
cols = np.any(img, axis=0)
Expand All @@ -238,6 +239,7 @@ def get_earth_radius(lon, lat, a, b):
Returns:
Earth Radius (meters)
"""
geocent = pyproj.Proj(proj='geocent', a=a, b=b, units='m')
latlong = pyproj.Proj(proj='latlong', a=a, b=b, units='m')
Expand All @@ -246,7 +248,7 @@ def get_earth_radius(lon, lat, a, b):


def reduce_mda(mda, max_size=100):
"""Recursively remove arrays with more than `max_size` elements from the given metadata dictionary"""
"""Recursively remove arrays with more than `max_size` elements from the given metadata dictionary."""
reduced = {}
for key, val in mda.items():
if isinstance(val, dict):
Expand Down
13 changes: 7 additions & 6 deletions satpy/tests/reader_tests/test_seviri_l1b_hrit.py
Original file line number Diff line number Diff line change
Expand Up @@ -251,12 +251,13 @@ def test_get_area_def(self):
area = self.reader.get_area_def(DatasetID('HRV'))
self.assertEqual(area.area_extent,
(-45561979844414.07, -3720765401003.719, 45602912357076.38, 77771774058.38356))
self.assertEqual(area.proj_dict, {'a': 6378169.0,
'b': 6356583.8,
'h': 35785831.0,
'lon_0': 44,
'proj': 'geos',
'units': 'm'})
proj_dict = area.proj_dict
self.assertEqual(proj_dict['a'], 6378169.0)
self.assertEqual(proj_dict['b'], 6356583.8)
self.assertEqual(proj_dict['h'], 35785831.0)
self.assertEqual(proj_dict['lon_0'], 44.0)
self.assertEqual(proj_dict['proj'], 'geos')
self.assertEqual(proj_dict['units'], 'm')
self.reader.fill_hrv = False
area = self.reader.get_area_def(DatasetID('HRV'))
self.assertEqual(area.defs[0].area_extent,
Expand Down
7 changes: 3 additions & 4 deletions satpy/tests/test_config.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,7 @@
#
# You should have received a copy of the GNU General Public License along with
# satpy. If not, see <http://www.gnu.org/licenses/>.
"""Test objects and functions in the satpy.config module.
"""
"""Test objects and functions in the satpy.config module."""

import sys

Expand Down Expand Up @@ -111,7 +110,7 @@ def test_areas_rasterio(self):
# in a generic test so just skip this area
continue
proj_dict = area_obj.proj_dict
if proj_dict['proj'] in ('ob_tran', 'nsper') and \
if proj_dict.get('proj') in ('ob_tran', 'nsper') and \
'wktext' not in proj_dict:
# FIXME: rasterio doesn't understand ob_tran unless +wktext
# See: https://github.com/pyproj4/pyproj/issues/357
Expand All @@ -121,7 +120,7 @@ def test_areas_rasterio(self):


def suite():
"""The test suite for test_config."""
"""Test suite for test_config."""
loader = unittest.TestLoader()
my_suite = unittest.TestSuite()
my_suite.addTest(loader.loadTestsFromTestCase(TestCheckSatpy))
Expand Down

0 comments on commit 1389229

Please sign in to comment.