Skip to content

Commit

Permalink
Merge pull request #141 from mhearne-usgs/shakefix
Browse files Browse the repository at this point in the history 
Shakefix
  • Loading branch information
@mhearne-usgs
mhearne-usgs authored May 8, 2020
2 parents ce270bb + c17c0e3 commit 2edf047
Show file tree
Hide file tree
Showing 2 changed files with 129 additions and 83 deletions.
123 changes: 64 additions & 59 deletions mapio/gridbase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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:
Expand All @@ -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:
Expand All @@ -60,45 +63,46 @@ 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:
Geodict that defines the new coarser grid.
"""
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:
Expand All @@ -109,41 +113,48 @@ 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:
Column dimension index into internal data array.
: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:
Expand All @@ -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')
89 changes: 65 additions & 24 deletions test/shake_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...')
Expand Down Expand Up @@ -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)...')
Expand All @@ -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'])

Expand All @@ -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],
Expand Down Expand Up @@ -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]])
Expand All @@ -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],
Expand All @@ -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()
Expand Down

0 comments on commit 2edf047

Please sign in to comment.