From d7c946840788ea22363259c4209c725681c2847e Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Thu, 7 May 2020 07:18:17 -0600 Subject: [PATCH 1/2] new test, saving here so I can share online --- test/shake_test.py | 89 +++++++++++++++++++++++++++++++++------------- 1 file changed, 65 insertions(+), 24 deletions(-) diff --git a/test/shake_test.py b/test/shake_test.py index 8ad5869..23ea4e0 100755 --- a/test/shake_test.py +++ b/test/shake_test.py @@ -4,36 +4,21 @@ from __future__ import print_function # stdlib imports -from xml.dom import minidom from datetime import datetime from collections import OrderedDict -import re -import sys import tempfile import time import shutil - -if sys.version_info.major == 2: - import StringIO -else: - from io import StringIO import os.path -# hack the path so that I can debug these functions if I need to -homedir = os.path.dirname(os.path.abspath(__file__)) # where is this script? -mapiodir = os.path.abspath(os.path.join(homedir, '..')) -# put this at the front of the system path, ignoring any installed mapio stuff -sys.path.insert(0, mapiodir) - # third party from mapio.shake import ShakeGrid -from mapio.gridbase import Grid -from mapio.multiple import MultiGrid -from mapio.dataset import DataSetException -from mapio.grid2d import Grid2D from mapio.geodict import GeoDict import numpy as np +homedir = os.path.dirname(os.path.abspath(__file__)) # where is this script? +mapiodir = os.path.abspath(os.path.join(homedir, '..')) + def test_modify(): print('Testing ShakeGrid interpolate() method...') @@ -177,7 +162,7 @@ def test_save(): print('Passed save/read functionality for shakemap grids.') print('Testing getFileGeoDict method...') - fgeodict = ShakeGrid.getFileGeoDict(testfile) + _ = ShakeGrid.getFileGeoDict(testfile) print('Passed save/read functionality for shakemap grids.') print('Testing loading with bounds (no resampling or padding)...') @@ -186,7 +171,8 @@ def test_save(): 'dx': 1.0, 'dy': 1.0, 'ny': 5, 'nx': 5}) shake3 = ShakeGrid.load(testfile, samplegeodict=sampledict, - resample=False, doPadding=False, padValue=np.nan) + resample=False, doPadding=False, + padValue=np.nan) tdata = shake3.getLayer('pga').getData() np.testing.assert_almost_equal(tdata, layers['pga']) @@ -198,7 +184,8 @@ def test_save(): 'dx': 1.0, 'dy': 1.0, 'ny': 6, 'nx': 6}) shake4 = ShakeGrid.load(testfile, samplegeodict=newdict, - resample=False, doPadding=True, padValue=np.nan) + resample=False, doPadding=True, + padValue=np.nan) output = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan, np.nan], [np.nan, 0.0, 1.0, 2.0, 3.0, np.nan], [np.nan, 4.0, 5.0, 6.0, 7.0, np.nan], @@ -230,7 +217,8 @@ def test_save(): 'dx': 1.0, 'dy': 1.0, 'ny': 3, 'nx': 3}) shake5 = ShakeGrid.load(testfile, samplegeodict=littledict, - resample=True, doPadding=False, padValue=np.nan) + resample=True, doPadding=False, + padValue=np.nan) output = np.array([[10.5, 11.5, 12.5], [16.5, 17.5, 18.5], [22.5, 23.5, 24.5]]) @@ -257,7 +245,8 @@ def test_save(): 'dx': 1.0, 'dy': 1.0, 'ny': 5, 'nx': 5}) shake6 = ShakeGrid.load( - testfile, samplegeodict=bigdict, resample=True, doPadding=True, padValue=np.nan) + testfile, samplegeodict=bigdict, resample=True, doPadding=True, + padValue=np.nan) tdata = shake6.getLayer('pga').getData() output = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan], [np.nan, 2.5, 3.5, 4.5, np.nan], @@ -268,14 +257,66 @@ def test_save(): print('Passed resampling and padding...') except Exception as error: print('Failed to read grid.xml format file "%s". Error "%s".' % - (xmlfile, str(error))) + (testfile, str(error))) assert 0 == 1 finally: if os.path.isdir(tdir): shutil.rmtree(tdir) +def test_meridian(): + shakeDict = {'event_id': 'usabcd1234', + 'shakemap_id': 'usabcd1234', + 'shakemap_version': 1, + 'code_version': '4.0', + 'process_timestamp': datetime.utcnow(), + 'shakemap_originator': 'us', + 'map_status': 'RELEASED', + 'shakemap_event_type': 'ACTUAL'} + eventDict = {'event_id': 'usabcd1234', + 'magnitude': 7.6, + 'depth': 1.4, + 'lat': 2.0, + 'lon': 2.0, + 'event_timestamp': datetime.utcnow(), + 'event_network': 'us', + # line below tests escaping XML + 'event_description': 'sample event & stuff'} + uncDict = {'pga': (0.0, 0), + 'pgv': (0.0, 0), + 'mmi': (0.0, 0)} + mdict = {'digits': 4, + 'dx': 0.033333333333333333, + 'dy': 0.033333333333333333, + 'nx': 442, + 'ny': 319, + 'units': 'ln(g)', + 'xmax': 180.34999999999999, + 'xmin': 165.65000000000001, + 'ymax': -36.649999999999999, + 'ymin': -47.25} + gdict = GeoDict(mdict) + pga = np.ones((gdict.ny, gdict.nx)) + pgv = np.ones((gdict.ny, gdict.nx)) + mmi = np.ones((gdict.ny, gdict.nx)) + layers = OrderedDict() + layers['pga'] = pga + layers['pgv'] = pgv + layers['mmi'] = mmi + shake = ShakeGrid(layers, gdict, eventDict, shakeDict, uncDict) + tdir = tempfile.mkdtemp() + try: + fname = os.path.join(tdir, 'temp.xml') + shake.save(fname) + x = 1 + except Exception: + assert False + finally: + shutil.rmtree(tdir) + + if __name__ == '__main__': + test_meridian() test_modify() test_interpolate() test_read() From c17c0e3ea913733440b3163d9fe59130b7684b28 Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Fri, 8 May 2020 11:42:28 -0600 Subject: [PATCH 2/2] Now correctly writing out shake grid.xml files when crossing the 180 meridian --- mapio/gridbase.py | 123 ++++++++++++++++++++++++---------------------- 1 file changed, 64 insertions(+), 59 deletions(-) diff --git a/mapio/gridbase.py b/mapio/gridbase.py index f8f3e7a..b862c31 100755 --- a/mapio/gridbase.py +++ b/mapio/gridbase.py @@ -3,17 +3,18 @@ import abc import numpy as np -#third party imports +# third party imports from .dataset import DataSet from .geodict import GeoDict + class Grid(DataSet): """ An abstract class to represent lat/lon gridded datasets. Grids are assumed to be pixel-registered - that is, grid coordinates represent the value at the *center* of the cells. """ - @abc.abstractmethod #should be a classmethod when instantiated + @abc.abstractmethod # should be a classmethod when instantiated def getFileGeoDict(filename): """ Abstract method to return the bounding box, resolution, and shape of a file in whatever Grid format. @@ -24,8 +25,8 @@ def getFileGeoDict(filename): """ raise NotImplementedError - @abc.abstractmethod #should be a classmethod when instantiated - def getBoundsWithin(filename,geodict): + @abc.abstractmethod # should be a classmethod when instantiated + def getBoundsWithin(filename, geodict): """ Abstract method to return a geodict for this file that is guaranteed to be inside the input geodict defined, without resampling. :param filename: @@ -36,21 +37,23 @@ def getBoundsWithin(filename,geodict): Always in base class """ raise NotImplementedError - + @classmethod - def _getPadding(cls,geodict,paddict,padvalue): - #get pad left columns - go outside specified bounds if not exact edge - pxmin,pxmax,pymin,pymax = (paddict.xmin,paddict.xmax,paddict.ymin,paddict.ymax) - gxmin,gxmax,gymin,gymax = (geodict.xmin,geodict.xmax,geodict.ymin,geodict.ymax) - dx,dy = (geodict.dx,geodict.dy) - ny,nx = (geodict.ny,geodict.nx) - - padleftcols = int(np.ceil((gxmin - pxmin)/dx)) - padrightcols = int(np.ceil((pxmax - gxmax)/dx)) - padbottomrows = int(np.ceil((gymin - pymin)/dy)) - padtoprows = int(np.ceil((pymax - gymax)/dy)) - - #if any of these are negative, set them to zero + def _getPadding(cls, geodict, paddict, padvalue): + # get pad left columns - go outside specified bounds if not exact edge + pxmin, pxmax, pymin, pymax = ( + paddict.xmin, paddict.xmax, paddict.ymin, paddict.ymax) + gxmin, gxmax, gymin, gymax = ( + geodict.xmin, geodict.xmax, geodict.ymin, geodict.ymax) + dx, dy = (geodict.dx, geodict.dy) + ny, nx = (geodict.ny, geodict.nx) + + padleftcols = int(np.ceil((gxmin - pxmin) / dx)) + padrightcols = int(np.ceil((pxmax - gxmax) / dx)) + padbottomrows = int(np.ceil((gymin - pymin) / dy)) + padtoprows = int(np.ceil((pymax - gymax) / dy)) + + # if any of these are negative, set them to zero if padleftcols < 0: padleftcols = 0 if padrightcols < 0: @@ -60,36 +63,37 @@ def _getPadding(cls,geodict,paddict,padvalue): if padtoprows < 0: padtoprows = 0 - leftpad = np.ones((ny,padleftcols))*padvalue - rightpad = np.ones((ny,padrightcols))*padvalue + leftpad = np.ones((ny, padleftcols)) * padvalue + rightpad = np.ones((ny, padrightcols)) * padvalue nx += padrightcols + padleftcols - bottompad = np.ones((padbottomrows,nx))*padvalue - toppad = np.ones((padtoprows,nx))*padvalue + bottompad = np.ones((padbottomrows, nx)) * padvalue + toppad = np.ones((padtoprows, nx)) * padvalue - #now figure out what the new bounds are + # now figure out what the new bounds are outdict = {} outdict['nx'] = int(nx) outdict['ny'] = int(ny + bottompad.shape[0] + toppad.shape[0]) - - outdict['xmin'] = gxmin - (padleftcols)*dx - outdict['xmax'] = gxmax + (padrightcols)*dx - outdict['ymin'] = gymin - (padbottomrows)*dy - outdict['ymax'] = gymax + (padtoprows)*dy + + outdict['xmin'] = gxmin - (padleftcols) * dx + outdict['xmax'] = gxmax + (padrightcols) * dx + outdict['ymin'] = gymin - (padbottomrows) * dy + outdict['ymax'] = gymax + (padtoprows) * dy outdict['dx'] = dx outdict['dy'] = dy - + gd = GeoDict(outdict) - return (leftpad,rightpad,bottompad,toppad,gd) - - @classmethod - def checkGeoDict(cls,geodict): - reqfields = set(['xmin','xmax','ymin','ymax','dx','dy','ny','nx']) + return (leftpad, rightpad, bottompad, toppad, gd) + + @classmethod + def checkGeoDict(cls, geodict): + reqfields = set( + ['xmin', 'xmax', 'ymin', 'ymax', 'dx', 'dy', 'ny', 'nx']) if not reqfields.issubset(set(geodict.keys())): return False return True - + @abc.abstractmethod - def blockmean(self,geodict): + def blockmean(self, geodict): """ Abstract method to calculate average values for cells of larger size than the current grid. :param geodict: @@ -97,8 +101,8 @@ def blockmean(self,geodict): """ raise NotImplementedError - @abc.abstractmethod #should be a classmethod when instantiated - def loadFromCloud(cls,cloud,geodict): + @abc.abstractmethod # should be a classmethod when instantiated + def loadFromCloud(cls, cloud, geodict): """ Create a grid from a Cloud instance (scattered XY data). :param cloud: @@ -109,28 +113,34 @@ def loadFromCloud(cls,cloud,geodict): An instance of a Grid object. """ raise NotImplementedError - + @staticmethod def getLatLonMesh(geodict): - lons = np.linspace(geodict.xmin,geodict.xmax,num=geodict.nx) - lats = np.linspace(geodict.ymin,geodict.ymax,num=geodict.ny) - lon,lat = np.meshgrid(lons,lats) - return (lat,lon) - + if geodict.xmax < geodict.xmin: + lons = np.linspace(geodict.xmin, + geodict.xmax + 360, + num=geodict.nx) + else: + lons = np.linspace(geodict.xmin, geodict.xmax, num=geodict.nx) + lats = np.linspace(geodict.ymin, geodict.ymax, num=geodict.ny) + lon, lat = np.meshgrid(lons, lats) + return (lat, lon) + @abc.abstractmethod def getGeoDict(self): """ Return a reference to the geodict inside the Grid - + :returns: A reference to a dictionary (see constructor). """ - raise NotImplementedError('getGeoDict method not implemented in base class') + raise NotImplementedError( + 'getGeoDict method not implemented in base class') @abc.abstractmethod - def getLatLon(self,row,col): + def getLatLon(self, row, col): """Return geographic coordinates (lat/lon decimal degrees) for given data row and column. - + :param row: Row dimension index into internal data array. :param col: @@ -138,12 +148,13 @@ def getLatLon(self,row,col): :returns: Tuple of latitude and longitude. """ - raise NotImplementedError('getLatLon method not implemented in base class') + raise NotImplementedError( + 'getLatLon method not implemented in base class') @abc.abstractmethod - def getRowCol(self,lat,lon,returnFloat=False): + def getRowCol(self, lat, lon, returnFloat=False): """Return data row and column from given geographic coordinates (lat/lon decimal degrees). - + :param lat: Input latitude. :param lon: @@ -153,11 +164,5 @@ def getRowCol(self,lat,lon,returnFloat=False): :returns: Tuple of row and column. """ - raise NotImplementedError('getRowCol method not implemented in base class') - - - - - - - + raise NotImplementedError( + 'getRowCol method not implemented in base class')