diff --git a/mapio/basemapcity.py b/mapio/basemapcity.py index 505bd18..cfd24d5 100644 --- a/mapio/basemapcity.py +++ b/mapio/basemapcity.py @@ -1,28 +1,27 @@ -#stdlib imports -import os.path -import warnings - -#local imports +# local imports from .mapcity import MapCities -from .dataset import DataSetException,DataSetWarning +from .dataset import DataSetException + +# third party imports -#third party imports -import matplotlib.font_manager import matplotlib.patheffects as path_effects import matplotlib.pyplot as plt import numpy as np import pandas as pd -XOFFSET = 4 #how many pixels between the city dot and the city text +XOFFSET = 4 # how many pixels between the city dot and the city text + class BasemapCities(MapCities): """ A subclass of Cities that can remove cities whose labels on a map intersect with larger cities. - """ - - def limitByMapCollision(self,basemap,fontname='Bitstream Vera Sans',fontsize=10.0): + """ + + def limitByMapCollision( + self, basemap, fontname="Bitstream Vera Sans", fontsize=10.0 + ): """Create a smaller Cities dataset by removing smaller cities whose bounding boxes collide with larger cities. @@ -39,7 +38,6 @@ def limitByMapCollision(self,basemap,fontname='Bitstream Vera Sans',fontsize=10. resolution='l',area_thresh=1000.,projection='lcc',\ lat_1=50.,lon_0=-107.,ax=ax) ####################################### - :param fontsize: Desired font size for city labels. :param fontname: @@ -55,46 +53,41 @@ def limitByMapCollision(self,basemap,fontname='Bitstream Vera Sans',fontsize=10. New Cities instance where smaller colliding cities have been removed. """ if basemap.ax is None: - raise DataSetException('Basemap is missing a reference to an "axes" object.') + raise DataSetException( + 'Basemap is missing a reference to an "axes" object.' + ) ax = basemap.ax - #get the transformation from display units (pixels) to data units - #we'll use this to set an offset in pixels between the city dot and - #the city name. + # get the transformation from display units (pixels) to data units + # we'll use this to set an offset in pixels between the city dot and + # the city name. self._display_to_data_transform = ax.transData.inverted() - + self.project(basemap) - if 'x' not in self._dataframe.columns and 'y' not in self._dataframe.columns: - raise DataSetException('Cities object has not had project() called yet.') - + if "x" not in self._dataframe.columns and "y" not in self._dataframe.columns: + raise DataSetException("Cities object has not had project() called yet.") + if fontname not in self._fontlist: - raise DataSetException('Font %s not in supported list.' % fontname) - #TODO: - check placement column + raise DataSetException("Font %s not in supported list." % fontname) + # TODO: - check placement column newdf = self._dataframe.copy() - #older versions of pandas use a different sort function - if pd.__version__ < '0.17.0': - newdf = newdf.sort(columns='pop',ascending=False) + # older versions of pandas use a different sort function + if pd.__version__ < "0.17.0": + newdf = newdf.sort(columns="pop", ascending=False) else: - newdf = newdf.sort_values(by='pop',ascending=False) + newdf = newdf.sort_values(by="pop", ascending=False) - lefts,rights,bottoms,tops = self._getCityBoundingBoxes(newdf,fontname,fontsize,ax) + lefts, rights, bottoms, tops = self._getCityBoundingBoxes( + newdf, fontname, fontsize, ax + ) ikeep = ~np.isnan(lefts) newdf = newdf.iloc[ikeep] lefts = lefts[ikeep] rights = rights[ikeep] bottoms = bottoms[ikeep] tops = tops[ikeep] - ikeep = [0] #indices of rows to keep in dataframe - for i in range(1,len(tops)): - cname = newdf.iloc[i]['name'] - if cname.lower().find('pacific grove') > -1: - allnames = newdf['name'].tolist() - sidx = allnames.index('Salinas') - sleft = lefts[sidx] - sright = rights[sidx] - sbottom = bottoms[sidx] - stop = tops[sidx] - foo = 1 + ikeep = [0] # indices of rows to keep in dataframe + for i in range(1, len(tops)): if np.isnan(lefts[i]): continue left = lefts[i] @@ -104,23 +97,23 @@ def limitByMapCollision(self,basemap,fontname='Bitstream Vera Sans',fontsize=10. clrx = (left > rights[0:i]) | (right < lefts[0:i]) clry = (top < bottoms[0:i]) | (bottom > tops[0:i]) - allclr = (clrx | clry) + allclr = clrx | clry if all(allclr): ikeep.append(i) - else: - foo = 1 if len(newdf): try: newdf = newdf.iloc[ikeep] - newdf['top'] = tops[ikeep] - newdf['bottom'] = bottoms[ikeep] - newdf['left'] = lefts[ikeep] - newdf['right'] = rights[ikeep] + newdf["top"] = tops[ikeep] + newdf["bottom"] = bottoms[ikeep] + newdf["left"] = lefts[ikeep] + newdf["right"] = rights[ikeep] except Exception as e: - x = 1 + print(e) return type(self)(newdf) - def renderToMap(self,ax,fontname='Bitstream Vera Sans',fontsize=10.0,zorder=10,shadow=False): + def renderToMap( + self, ax, fontname="Bitstream Vera Sans", fontsize=10.0, zorder=10, shadow=False + ): """Render cities on Basemap axes. :param ax: @@ -130,23 +123,21 @@ def renderToMap(self,ax,fontname='Bitstream Vera Sans',fontsize=10.0,zorder=10,s :param fontsize: Font size in points. :param zorder: - Matplotlib plotting order - higher zorder is on top. + Matplotlib plotting order - higher zorder is on top. :param shadow: Boolean indicating whether "drop-shadow" effect should be used. """ - if 'x' not in self._dataframe.columns and 'y' not in self._dataframe.columns: - raise DataSetException('Cities object has not had project() called yet.') - #get the transformation from display units (pixels) to data units - #we'll use this to set an offset in pixels between the city dot and - #the city name. + if "x" not in self._dataframe.columns and "y" not in self._dataframe.columns: + raise DataSetException("Cities object has not had project() called yet.") + # get the transformation from display units (pixels) to data units + # we'll use this to set an offset in pixels between the city dot and + # the city name. self._display_to_data_transform = ax.transData.inverted() - - for index,row in self._dataframe.iterrows(): - th = self._renderRow(row,ax,fontname,fontsize,shadow=shadow, - zorder=zorder) - ax.plot(row['x'],row['y'],'k.') - - def _getCityBoundingBoxes(self,df,fontname,fontsize,ax): + + for index, row in self._dataframe.iterrows(): + ax.plot(row["x"], row["y"], "k.") + + def _getCityBoundingBoxes(self, df, fontname, fontsize, ax): """Get the axes coordinate system bounding boxes for each city. :param df: DataFrame containing information about cities. @@ -160,30 +151,34 @@ def _getCityBoundingBoxes(self,df,fontname,fontsize,ax): Numpy arrays of top,bottom,left and right edges of city bounding boxes. """ fig = ax.get_figure() - fwidth,fheight = fig.get_figwidth(),fig.get_figheight() + fwidth, fheight = fig.get_figwidth(), fig.get_figheight() plt.sca(ax) - axmin,axmax,aymin,aymax = plt.axis() + axmin, axmax, aymin, aymax = plt.axis() axbox = ax.get_position().bounds - newfig = plt.figure(figsize=(fwidth,fheight)) + newfig = plt.figure(figsize=(fwidth, fheight)) newax = newfig.add_axes(axbox) newfig.canvas.draw() plt.sca(newax) - plt.axis((axmin,axmax,aymin,aymax)) - #make arrays of the edges of all the bounding boxes - tops = np.ones(len(df))*np.nan - bottoms = np.ones(len(df))*np.nan - lefts = np.ones(len(df))*np.nan - rights = np.ones(len(df))*np.nan + plt.axis((axmin, axmax, aymin, aymax)) + # make arrays of the edges of all the bounding boxes + tops = np.ones(len(df)) * np.nan + bottoms = np.ones(len(df)) * np.nan + lefts = np.ones(len(df)) * np.nan + rights = np.ones(len(df)) * np.nan if len(df): - left,right,bottom,top = self._getCityEdges(df.iloc[0],newax,newfig,fontname,fontsize) + left, right, bottom, top = self._getCityEdges( + df.iloc[0], newax, newfig, fontname, fontsize + ) lefts[0] = left rights[0] = right bottoms[0] = bottom tops[0] = top - for i in range(1,len(df)): + for i in range(1, len(df)): row = df.iloc[i] - left,right,bottom,top = self._getCityEdges(row,newax,newfig,fontname,fontsize) - #remove cities that have any portion off the map + left, right, bottom, top = self._getCityEdges( + row, newax, newfig, fontname, fontsize + ) + # remove cities that have any portion off the map if left < axmin or right > axmax or bottom < aymin or top > aymax: continue lefts[i] = left @@ -191,11 +186,11 @@ def _getCityBoundingBoxes(self,df,fontname,fontsize,ax): bottoms[i] = bottom tops[i] = top - #get rid of new figure and its axes + # get rid of new figure and its axes plt.close(newfig) - return (lefts,rights,bottoms,tops) + return (lefts, rights, bottoms, tops) - def _renderRow(self,row,ax,fontname,fontsize,zorder=10,shadow=False): + def _renderRow(self, row, ax, fontname, fontsize, zorder=10, shadow=False): """Internal method to consistently render city names. :param row: pandas dataframe row. @@ -206,48 +201,70 @@ def _renderRow(self,row,ax,fontname,fontsize,zorder=10,shadow=False): :param fontsize: Font size in points. :param zorder: - Matplotlib plotting order - higher zorder is on top. + Matplotlib plotting order - higher zorder is on top. :param shadow: Boolean indicating whether "drop-shadow" effect should be used. :returns: Matplotlib Text instance. """ - ha = 'left' - va = 'center' - if 'placement' in row.index: - if row['placement'].find('E') > -1: - ha = 'left' - if row['placement'].find('W') > -1: - ha = 'right' + ha = "left" + va = "center" + if "placement" in row.index: + if row["placement"].find("E") > -1: + ha = "left" + if row["placement"].find("W") > -1: + ha = "right" else: - ha = 'center' - if row['placement'].find('N') > -1: - ha = 'top' - if row['placement'].find('S') > -1: - ha = 'bottom' + ha = "center" + if row["placement"].find("N") > -1: + ha = "top" + if row["placement"].find("S") > -1: + ha = "bottom" else: - ha = 'center' - - display1 = (1,1) - display2 = (1+XOFFSET,1) + ha = "center" + + display1 = (1, 1) + display2 = (1 + XOFFSET, 1) data1 = self._display_to_data_transform.transform((display1)) data2 = self._display_to_data_transform.transform((display2)) data_x_offset = data2[0] - data1[0] - tx = row['x'] + data_x_offset - ty = row['y'] - if shadow: - th = ax.text(tx,ty,row['name'],fontname=fontname,color='white', - fontsize=fontsize,ha=ha,va=va,zorder=zorder,clip_on=True) - th.set_path_effects([path_effects.Stroke(linewidth=1.5, foreground='black'), - path_effects.Normal()]) - else: - th = ax.text(tx,ty,row['name'],fontname=fontname, - fontsize=fontsize,ha=ha,va=va,zorder=zorder,clip_on=True) - + tx = row["x"] + data_x_offset + ty = row["y"] + if shadow: + th = ax.text( + tx, + ty, + row["name"], + fontname=fontname, + color="white", + fontsize=fontsize, + ha=ha, + va=va, + zorder=zorder, + clip_on=True, + ) + th.set_path_effects( + [ + path_effects.Stroke(linewidth=1.5, foreground="black"), + path_effects.Normal(), + ] + ) + else: + th = ax.text( + tx, + ty, + row["name"], + fontname=fontname, + fontsize=fontsize, + ha=ha, + va=va, + zorder=zorder, + clip_on=True, + ) + return th - - - def _getCityEdges(self,row,ax,fig,fontname,fontsize): + + def _getCityEdges(self, row, ax, fig, fontname, fontsize): """Return the edges of a city label on a given map axes. :param row: Row of a dataframe containing city information. @@ -260,11 +277,8 @@ def _getCityEdges(self,row,ax,fig,fontname,fontsize): :param fontsize: Font size in points. """ - th = self._renderRow(row,ax,fontname,fontsize) + th = self._renderRow(row, ax, fontname, fontsize) bbox = th.get_window_extent(fig.canvas.renderer) axbox = bbox.inverse_transformed(ax.transData) - left,bottom,right,top = axbox.extents - return (left,right,bottom,top) - - - + left, bottom, right, top = axbox.extents + return (left, right, bottom, top) diff --git a/mapio/city.py b/mapio/city.py index ea8ce62..9cb2651 100755 --- a/mapio/city.py +++ b/mapio/city.py @@ -1,9 +1,9 @@ #!/usr/bin/python -from xml.dom.minidom import parse import numpy as np import zipfile import tempfile import sys + if sys.version_info.major == 3: import urllib.request as request else: @@ -13,14 +13,16 @@ from impactutils.extern.openquake.geodetic import geodetic_distance import pandas as pd -from .dataset import DataSetException,DataSetWarning +from .dataset import DataSetException + +GEONAME_URL = "http://download.geonames.org/export/dump/cities1000.zip" -GEONAME_URL = 'http://download.geonames.org/export/dump/cities1000.zip' def _fetchGeoNames(): """ Internal method to retrieve a cities1000.txt file from GeoNames. - :returns: + + Returns: Path to local cities1000.txt file. """ fh = request.urlopen(GEONAME_URL) @@ -29,43 +31,51 @@ def _fetchGeoNames(): f = io.BytesIO(data) myzip = zipfile.ZipFile(f) fdir = tempfile.mkdtemp() - myzip.extract('cities1000.txt',fdir) + myzip.extract("cities1000.txt", fdir) myzip.close() - return os.path.join(fdir,'cities1000.txt') + return os.path.join(fdir, "cities1000.txt") + class Cities(object): """ Handles loading and searching for cities. """ - REQFIELDS = ['name','lat','lon'] #class variable - def __init__(self,dataframe): + + REQFIELDS = ["name", "lat", "lon"] # class variable + + def __init__(self, dataframe): """Construct a Cities object from a pandas DataFrame. - :param dataframe: - pandas DataFrame, where each row represents a city. - Columns include: - - name Name of the city (required). - - lat Latitude of city (required). - - lon Longitude of city (required). - - pop Population of city (optional). - - iscap Boolean indicating capital status (optional). - - placement String indicating where city label - should be placed relative to city coordinates, - one of: E,W,N,S,NE,SE,SW,NW (optional). - -xoff Longitude offset for label relative to city coordinates (optional). - -yoff Latitude offset for label relative to city coordinates (optional). - :raises DataSetException: - When any of required columns are missing. - :returns: - Cities instance. + + Returns: + dataframe (pandas dataframe): + where each row represents a city. + Columns include: + - name Name of the city (required). + - lat Latitude of city (required). + - lon Longitude of city (required). + - pop Population of city (optional). + - iscap Boolean indicating capital status (optional). + - placement String indicating where city label + should be placed relative to city coordinates, + one of: E,W,N,S,NE,SE,SW,NW (optional). + -xoff Longitude offset for label relative to city coordinates + (optional). + -yoff Latitude offset for label relative to city coordinates (optional). + Raises: + DataSetException: + When any of required columns are missing. + Returns: + Cities instance. """ if len(set(dataframe.columns).intersection(set(self.REQFIELDS))) < 3: - raise DataSetException('Missing some of required keys: %s' % self.REQFIELDS) + raise DataSetException("Missing some of required keys: %s" % self.REQFIELDS) self._dataframe = dataframe.copy() - #"magic" methods + # "magic" methods def __len__(self): """Return the number of cities in the Cities object. - :returns: + + Returns: Number of cities in the Cities object. """ return len(self._dataframe) @@ -73,201 +83,229 @@ def __len__(self): def __repr__(self): """ Return the string to represent the Cities instance. - :returns: + + Returns: String representing Cities instance. """ return str(self._dataframe) - + @classmethod - def loadFromGeoNames(cls,cityfile=None): - """Load a cities data set from a GeoNames cities1000.txt file or by downloading + def loadFromGeoNames(cls, cityfile=None): + """Load a cities data set from a GeoNames cities1000.txt file or by downloading it from GeoNames, then loading it. - :param cityfile: - Path to cities1000.txt file from GeoNames, or None (file will be downloaded). - :returns: - Cities instance. + Args: + cityfile (str): + Path to cities1000.txt file from GeoNames, or None + (file will be downloaded). + Returns: + Cities instance. """ - CAPFLAG1 = 'PPLC' - CAPFLAG2 = 'PPLA' + CAPFLAG1 = "PPLC" + CAPFLAG2 = "PPLA" delete_folder = False if cityfile is None: cityfile = _fetchGeoNames() delete_folder = True - mydict = {'name':[], - 'ccode':[], - 'lat':[], - 'lon':[], - 'iscap':[], - 'pop':[]} - f = open(cityfile,'rb') + mydict = {"name": [], "ccode": [], "lat": [], "lon": [], "iscap": [], "pop": []} + f = open(cityfile, "rb") for line in f.readlines(): - line = line.decode('utf-8') - parts = line.split('\t') + line = line.decode("utf-8") + parts = line.split("\t") tname = parts[2].strip() if not tname: continue myvals = np.array([ord(c) for c in tname]) if len((myvals > 127).nonzero()[0]): continue - mydict['name'].append(tname) - mydict['ccode'].append(parts[8].strip()) - mydict['lat'].append(float(parts[4].strip())) - mydict['lon'].append(float(parts[5].strip())) + mydict["name"].append(tname) + mydict["ccode"].append(parts[8].strip()) + mydict["lat"].append(float(parts[4].strip())) + mydict["lon"].append(float(parts[5].strip())) capfield = parts[7].strip() iscap = (capfield == CAPFLAG1) or (capfield == CAPFLAG2) - mydict['iscap'].append(iscap) - mydict['pop'].append(int(parts[14].strip())) + mydict["iscap"].append(iscap) + mydict["pop"].append(int(parts[14].strip())) f.close() if delete_folder: - fdir,bname = os.path.split(cityfile) + fdir, bname = os.path.split(cityfile) os.remove(cityfile) os.rmdir(fdir) df = pd.DataFrame.from_dict(mydict) return cls(df) @classmethod - def loadFromCSV(cls,csvfile): + def loadFromCSV(cls, csvfile): """Load data from a csv file - :param csvfile: - CSV file containing city data, must have at least columns name,lat,lon. - :returns: - Cities instance + + Args: + csvfile (str): + CSV file containing city data, must have at least columns name,lat,lon. + Returns: + Cities instance """ df = pd.read_csv(csvfile) return cls(df) - - def save(self,filename): + + def save(self, filename): """Save City internal dataframe to CSV file. - - :param filename: - Output filename to save data to. - :returns: - None + + Args: + filename (str): + Output filename to save data to. + Returns: + None """ self._dataframe.to_csv(filename) - def sortByColumns(self,columns,ascending=True): + def sortByColumns(self, columns, ascending=True): """Sort list of cities by input column names. - :param columns: - String name or list of names of any of the columns that are in the internal - dataframe. See getColumns(). Only the required set of columns (see __init__ method) - are guaranteed to be present, but subclasses of Cities may add more. - :param ascending: - Boolean indicating which direction values should be sorted. - :raises DataSetException: - When column(s) are not in the list of dataframe columns. + Args: + columns (str or :tuple:'list'): + String name or list of names of any of the columns that are in the + internal dataframe. See getColumns(). + Only the required set of columns (see __init__ method) + are guaranteed to be present, but subclasses of Cities may add more. + ascending (bool): + indicating which direction values should be sorted. + Raises: + DataSetException: + When column(s) are not in the list of dataframe columns. """ - if isinstance(columns,str): + if isinstance(columns, str): columns = [columns] for column in columns: if column not in self._dataframe.columns: - raise DataSetException('Column not in list of DataFrame columns') - if pd.__version__ < '0.17.0': - self._dataframe = self._dataframe.sort(columns=columns,ascending=ascending) + raise DataSetException("Column not in list of DataFrame columns") + if pd.__version__ < "0.17.0": + self._dataframe = self._dataframe.sort(columns=columns, ascending=ascending) else: - self._dataframe = self._dataframe.sort_values(by=columns,ascending=ascending) + self._dataframe = self._dataframe.sort_values( + by=columns, ascending=ascending + ) def getColumns(self): """Return list of column names in internal data frame. - :returns: - List of column names in internal data frame. + Returns: + List of column names in internal data frame. """ return list(self._dataframe.columns) - def limitByBounds(self,bounds): + def limitByBounds(self, bounds): """Search for cities within a bounding box (xmin,xmax,ymin,ymax). - :param bounds: - Sequence containing xmin,xmax,ymin,ymax (decimal degrees). - :returns: - New Cities instance containing smaller cities data set. + + Args: + bounds (sequence): + Sequence containing xmin,xmax,ymin,ymax (decimal degrees). + Returns: + New Cities instance containing smaller cities data set. """ - #TODO - figure out what to do with a meridian crossing? + # TODO - figure out what to do with a meridian crossing? newdf = self._dataframe.copy() - xmin,xmax,ymin,ymax = bounds - newdf = newdf.loc[(newdf['lat'] > ymin) & (newdf['lat'] <= ymax) & - (newdf['lon'] > xmin) & (newdf['lon'] <= xmax)] + xmin, xmax, ymin, ymax = bounds + newdf = newdf.loc[ + (newdf["lat"] > ymin) + & (newdf["lat"] <= ymax) + & (newdf["lon"] > xmin) + & (newdf["lon"] <= xmax) + ] return type(self)(newdf) - def limitByRadius(self,lat,lon,radius): + def limitByRadius(self, lat, lon, radius): """Search for cities within a radius (km) around a central point. - :param lat: - Central latitude coordinate (dd). - :param lon: - Central longitude coordinate (dd). - :param radius: - Radius (km) around which cities will be searched. - :returns: - New Cities instance containing smaller cities data set. + + Args: + lat (float): + Central latitude coordinate (dd). + lon (float): + Central longitude coordinate (dd). + radius (foat): + Radius (km) around which cities will be searched. + Returns: + New Cities instance containing smaller cities data set. """ - #TODO - figure out what to do with a meridian crossing? + # TODO - figure out what to do with a meridian crossing? newdf = self._dataframe.copy() - dist = geodetic_distance(lon,lat,self._dataframe['lon'],self._dataframe['lat']) - newdf['dist'] = dist - newdf = newdf[newdf['dist'] <= radius] - del newdf['dist'] + dist = geodetic_distance( + lon, lat, self._dataframe["lon"], self._dataframe["lat"] + ) + newdf["dist"] = dist + newdf = newdf[newdf["dist"] <= radius] + del newdf["dist"] return type(self)(newdf) - def limitByPopulation(self,pop,minpop=0): + def limitByPopulation(self, pop, minpop=0): """Search for cities above a certain population threshold. - :param pop: - Population threshold. - :param minpop: - Population above which cities should be included. - :raises DataSetException: - When Cities instance does not contain population data. - When minpop >= pop. - :returns: - New Cities instance containing cities where population > pop. + + Args: + pop (float): + Population threshold. + minpop (:obj:`float`, optional): + Population above which cities should be included. + Raises: + DataSetException: + When Cities instance does not contain population data. + When minpop >= pop. + Returns: + New Cities instance containing cities where population > pop. """ - if 'pop' not in self._dataframe.columns: - raise DataSetException('Cities instance does not contain population information') + if "pop" not in self._dataframe.columns: + raise DataSetException( + "Cities instance does not contain population information" + ) if minpop >= pop: - raise DataSetException('Minimum population must be less than population threshold.') + raise DataSetException( + "Minimum population must be less than population threshold." + ) newdf = self._dataframe.copy() - newdf = newdf[newdf['pop'] >= pop] + newdf = newdf[newdf["pop"] >= pop] return type(self)(newdf) - def limitByGrid(self,nx=2,ny=2,cities_per_grid=20): + def limitByGrid(self, nx=2, ny=2, cities_per_grid=20): """Create a smaller Cities dataset by gridding cities, then limiting cities in each grid by population. - :param nx: - Desired number of columns for grid. - :param ny: - Desired number of rows for grid. - :param cities_per_cell: - Maximum number of cities allowed per grid cell. - :raises DataSetException: - When Cities instance does not contain population data. - :returns: - New Cities instance containing cities limited by number in each grid cell. + + Args: + nx (int): + Desired number of columns for grid. + ny (int): + Desired number of rows for grid. + cities_per_cell (int): + Maximum number of cities allowed per grid cell. + Raises: + DataSetException: + When Cities instance does not contain population data. + Returns: + New Cities instance containing cities limited by number in each grid cell. """ - if 'pop' not in self._dataframe.columns: - raise DataSetException('Cities instance does not contain population information') - xmin = self._dataframe['lon'].min() - xmax = self._dataframe['lon'].max() - ymin = self._dataframe['lat'].min() - ymax = self._dataframe['lat'].max() - dx = (xmax-xmin)/nx - dy = (ymax-ymin)/ny + if "pop" not in self._dataframe.columns: + raise DataSetException( + "Cities instance does not contain population information" + ) + xmin = self._dataframe["lon"].min() + xmax = self._dataframe["lon"].max() + ymin = self._dataframe["lat"].min() + ymax = self._dataframe["lat"].max() + dx = (xmax - xmin) / nx + dy = (ymax - ymin) / ny newdf = None - #start from the bottom left of our grid, and trim our cities. - - for i in range(0,ny): - cellymin = ymin + (i*dy) + # start from the bottom left of our grid, and trim our cities. + + for i in range(0, ny): + cellymin = ymin + (i * dy) cellymax = cellymin + dy - for j in range(0,nx): - cellxmin = xmin + (j*dx) + for j in range(0, nx): + cellxmin = xmin + (j * dx) cellxmax = cellxmin + dx - tcities = self.limitByBounds((cellxmin,cellxmax,cellymin,cellymax)) - #older versions of pandas use a different sort function - if pd.__version__ < '0.17.0': - tdf = tcities._dataframe.sort(columns='pop',ascending=False) + tcities = self.limitByBounds((cellxmin, cellxmax, cellymin, cellymax)) + # older versions of pandas use a different sort function + if pd.__version__ < "0.17.0": + tdf = tcities._dataframe.sort(columns="pop", ascending=False) else: - tdf = tcities._dataframe.sort_values(by='pop',ascending=False) + tdf = tcities._dataframe.sort_values(by="pop", ascending=False) tdf = tdf[0:cities_per_grid] if newdf is None: newdf = tdf.copy() @@ -275,55 +313,58 @@ def limitByGrid(self,nx=2,ny=2,cities_per_grid=20): newdf = newdf.append(tdf) return type(self)(newdf) - def limitByName(self,cityname): + def limitByName(self, cityname): """Find all cities that match a given cityname (or regular expression). - :param cityname: - Input city name (i.e., "Los Angeles"). - :returns: - Cities instance containing cities with names that match the input name/regular expression. + + Args: + cityname (str): + Input city name (i.e., "Los Angeles"). + Returns: + Cities instance containing cities + with names that match the input name/regular expression. """ newdf = self._dataframe[self._dataframe.name.str.contains(cityname)] return type(self)(newdf) def getDataFrame(self): """Return a copy of the internal pandas DataFrame containing city data. - :returns: - pandas DataFrame at least containing columns 'name','lat','lon' + + Returns: + pandas DataFrame at least containing columns 'name','lat','lon' """ return self._dataframe.copy() - + def getBounds(self): """Return the bounds of the Cities dataset. - :returns: - Tuple containing (xmin,xmax,ymin,ymax). + + Returns: + Tuple containing (xmin,xmax,ymin,ymax). """ - #TODO - figure out meridian crossing?? - xmin = self._dataframe['lon'].min() - xmax = self._dataframe['lon'].max() - ymin = self._dataframe['lat'].min() - ymax = self._dataframe['lat'].max() - return (xmin,xmax,ymin,ymax) + # TODO - figure out meridian crossing?? + xmin = self._dataframe["lon"].min() + xmax = self._dataframe["lon"].max() + ymin = self._dataframe["lat"].min() + ymax = self._dataframe["lat"].max() + return (xmin, xmax, ymin, ymax) def getCities(self): """Return arrays of lat,lon,names from Cities DataFrame. - :returns: - tuple of (lat,lon,names) where each is a numpy array. + + Returns: + tuple of (lat,lon,names) where each is a numpy array. """ - lat = self._dataframe['lat'].values - lon = self._dataframe['lon'].values - names = self._dataframe['name'].values - return (lat,lon,names) + lat = self._dataframe["lat"].values + lon = self._dataframe["lon"].values + names = self._dataframe["name"].values + return (lat, lon, names) - def project(self,mbasemap): + def project(self, mbasemap): """Use a Basemap instance to project city lat/lon to projected x/y. - :param mbasemap: - Basemap instance. + + Args: + mbasemap (basemap): + Basemap instance. """ - x,y = mbasemap(self._dataframe['lon'].values,self._dataframe['lat'].values) - self._dataframe['x'] = x - self._dataframe['y'] = y - - - - - + x, y = mbasemap(self._dataframe["lon"].values, self._dataframe["lat"].values) + self._dataframe["x"] = x + self._dataframe["y"] = y diff --git a/mapio/dataset.py b/mapio/dataset.py index 712f4d8..849455b 100755 --- a/mapio/dataset.py +++ b/mapio/dataset.py @@ -1,145 +1,174 @@ #!/usr/bin/env python -#python 3 compatibility +# python 3 compatibility from __future__ import print_function -#stdlib imports +# stdlib imports import abc - class DataSetException(Exception): """ Class to represent errors in the DataSet class. """ + def __init__(self, value): self.value = value + def __str__(self): return repr(self.value) + class DataSetWarning(Warning): """ Class to represent warnings in the DataSet class. """ + def __init__(self, value): self.value = value + def __str__(self): return repr(self.value) + class DataSet(object): - #This should be a @classmethod in subclasses + # This should be a @classmethod in subclasses @abc.abstractmethod - def load(filename,bounds=None,resample=False,padValue=None): + def load(filename, bounds=None, resample=False, padValue=None): """ Load data into a Grid subclass. Parameters below are suggested for subclasses. - :param filename: - File where data is stored - :param bounds: - Optional tuple of (lonmin,lonmax,latmin,latmax) used to subset data from file. - :param resample: - If subsetting data, True indicates that *exact* bounds are desired and data should be resampled to fit. - :param padValue: - If asking for data outside bounds of grid, any value not None means fill in those cells with padValue. - None means don't pad the grid at all. - :raises NotImplementedError: - Always for this base class. + Args: + filename: + File where data is stored + bounds: + Optional tuple of (lonmin,lonmax,latmin,latmax) to subset data from + file. + resample: + If subsetting data, True indicates that *exact* bounds are desired and + data should be resampled to fit. + padValue: + If asking for data outside bounds of grid, any value not None means + fill in those cells with padValue. None means don't pad the grid at all. + Raises: + NotImplementedError: + Always for this base class. """ - raise NotImplementedError('load method not implemented in base class') + raise NotImplementedError("load method not implemented in base class") - #TODO: Figure out how format-specific attributes will be handled (ShakeMap, for example) - #This should be a @classmethod in subclasses + # TODO: Figure out how format-specific attributes will be handled + # (ShakeMap, for example) + # This should be a @classmethod in subclasses @abc.abstractmethod - def save(self,filename): #would we ever want to save a subset of the data? + def save(self, filename): # would we ever want to save a subset of the data? """ - Save the data contained in the grid to a format specific file. Other attributes may be required for - format specific files. + Save the data contained in the grid to a format specific file. + Other attributes may be required for format specific files. - :param filename: - Where file containing data should be written. + Args: + filename: + Where file containing data should be written. """ - raise NotImplementedError('Save method not implemented in base class') + raise NotImplementedError("Save method not implemented in base class") @abc.abstractmethod - def getData(self,getCopy=False): + def getData(self, getCopy=False): """ Return a reference to or copy of the data inside the Grid - :param getCopy: - True indicates that the user wants a copy of the data, not a reference to it. - :returns: - A reference to or copy of a numpy array of data. + Args: + getCopy: + True indicates that the user wants a copy of the data, + not a reference to it. + Returns: + A reference to or copy of a numpy array of data. """ - raise NotImplementedError('getData method not implemented in base class') + raise NotImplementedError("getData method not implemented in base class") @abc.abstractmethod - def setData(self,data): + def setData(self, data): """ Modify the data inside the Grid. - :param data: - numpy array of desired data. + Args: + data: + numpy array of desired data. """ - raise NotImplementedError('setData method not implemented in base class') + raise NotImplementedError("setData method not implemented in base class") @abc.abstractmethod def getBounds(self): """ Return the lon/lat range of the data. - - :returns: + + Returns: Tuple of (lonmin,lonmax,latmin,latmax) """ - raise NotImplementedError('getBounds method not implemented in base class') + raise NotImplementedError("getBounds method not implemented in base class") @abc.abstractmethod - def trim(self,geodict,resample=False,method='linear'): + def trim(self, geodict, resample=False, method="linear"): """ - Trim data to a smaller set of bounds, resampling if requested. If not resampling, - data will be trimmed to smallest grid boundary possible. - - :param geodict: - GeoDict object used to specify subset bounds and resolution (if resample is selected) - :param resample: - Boolean indicating whether the data should be resampled to *exactly* match input bounds. - :param method: - If resampling, method used, one of ('linear','nearest','cubic','quintic') + Trim data to a smaller set of bounds, resampling if requested. + If not resampling, data will be trimmed to smallest grid boundary possible. + + Args: + geodict: + GeoDict object used to specify subset bounds and resolution + (if resample is selected) + resample: + Boolean indicating whether the data should be resampled to *exactly* + match input bounds. + method: + If resampling, method used, one of + ('linear','nearest','cubic','quintic') """ - raise NotImplementedError('trim method not implemented in base class') + raise NotImplementedError("trim method not implemented in base class") @abc.abstractmethod - def getValue(self,lat,lon,method='nearest',default=None): #return nearest neighbor value + def getValue( + self, lat, lon, method="nearest", default=None + ): # return nearest neighbor value """Return numpy array at given latitude and longitude (using given resampling method). - - :param lat: - Latitude (in decimal degrees) of desired data value. - :param lon: - Longitude (in decimal degrees) of desired data value. - :param method: - Interpolation method, one of ('nearest','linear','cubic','quintic') - :param default: - Default value to return when lat/lon is outside of grid bounds. - :return: + + Arg: + lat: + Latitude (in decimal degrees) of desired data value. + Arg: + lon: + Longitude (in decimal degrees) of desired data value. + method: + Interpolation method, one of ('nearest','linear','cubic','quintic') + default: + Default value to return when lat/lon is outside of grid bounds. + Returns: Value at input latitude,longitude position. """ - raise NotImplementedError('getValue method not implemented in base class') + raise NotImplementedError("getValue method not implemented in base class") @abc.abstractmethod - def interpolateToGrid(self,geodict,method='linear'): + def interpolateToGrid(self, geodict, method="linear"): """ - Given a geodict specifying a grid extent and resolution, resample current data set to match. - - :param geodict: - geodict object from a grid whose extents are inside the extent of this grid. - :param method: - Optional interpolation method - ['linear', 'cubic','quintic','nearest'] - :returns: - Interpolated grid. - :raises DataSetException: - If the Grid object upon which this function is being called is not completely contained by the grid to which this Grid is being resampled. - :raises DataSetException: - If the resulting interpolated grid shape does not match input geodict. + Given a geodict specifying a grid extent and resolution, resample current data + set to match. + + Args: + geodict: + geodict object from a grid whose extents are inside the extent of this + grid. + method: + Optional interpolation method - ['linear', 'cubic','quintic','nearest'] + Returns: + Interpolated grid. + Raises: + DataSetException: + If the Grid object upon which this function is being called is not + completely contained by the grid to which this Grid is being resampled. + DataSetException: + If the resulting interpolated grid shape does not match input geodict. This function modifies the internal griddata and geodict object variables. """ - raise NotImplementedError('interpolateToGrid method not implemented in base class') + raise NotImplementedError( + "interpolateToGrid method not implemented in base class" + ) diff --git a/mapio/gdal.py b/mapio/gdal.py index b5fdfe6..b7e284c 100644 --- a/mapio/gdal.py +++ b/mapio/gdal.py @@ -32,26 +32,35 @@ def get_affine(src): class GDALGrid(Grid2D): def __init__(self, data, geodict): """Construct a GMTGrid object. - :param data: - 2D numpy data array (must match geodict spec) - :param geodict: - GeoDict Object specifying the spatial extent,resolution and shape - of the data. - :returns: - A GMTGrid object. - :raises DataSetException: - When data and geodict dimensions do not match. + + Args: + data (np array): + 2D numpy data array (must match geodict spec) + geodict (GeoDict): + GeoDict Object specifying the spatial extent,resolution and shape + of the data. + Returns: + A GMTGrid object. + Raises: + DataSetException: + When data and geodict dimensions do not match. """ m, n = data.shape if m != geodict.ny or n != geodict.nx: - raise DataSetException( - 'Input geodict does not match shape of input data.') + raise DataSetException("Input geodict does not match shape of input data.") self._data = data self._geodict = geodict @classmethod - def load(cls, filename, samplegeodict=None, resample=False, - method='linear', doPadding=False, padValue=np.nan): + def load( + cls, + filename, + samplegeodict=None, + resample=False, + method="linear", + doPadding=False, + padValue=np.nan, + ): """ This method should do the following: 1) If resampling, buffer bounds outwards. @@ -78,20 +87,21 @@ def load(cls, filename, samplegeodict=None, resample=False, # filegeodict.isAligned(samplegeodict): # resample = False - if samplegeodict is not None and \ - not filegeodict.intersects(samplegeodict): + if samplegeodict is not None and not filegeodict.intersects(samplegeodict): if not doPadding: - raise DataSetException('Your sampling grid is not contained ' - 'by the file. To load anyway, ' - 'use doPadding=True.') + raise DataSetException( + "Your sampling grid is not contained " + "by the file. To load anyway, " + "use doPadding=True." + ) # buffer out the sample geodict (if resampling) enough to allow # interpolation. if samplegeodict is not None: # parent static method - sampledict = cls.bufferBounds(samplegeodict, filegeodict, - resample=resample, - doPadding=doPadding) + sampledict = cls.bufferBounds( + samplegeodict, filegeodict, resample=resample, doPadding=doPadding + ) else: sampledict = filegeodict @@ -105,8 +115,10 @@ def load(cls, filename, samplegeodict=None, resample=False, except DataSetException: if doPadding: if not filegeodict.contains(sampledict): - data = np.ones((sampledict.ny, sampledict.nx), - dtype=np.float32)*padValue + data = ( + np.ones((sampledict.ny, sampledict.nx), dtype=np.float32) + * padValue + ) return cls(data=data, geodict=sampledict) sampledict = filegeodict.getIntersection(sampledict) @@ -114,14 +126,13 @@ def load(cls, filename, samplegeodict=None, resample=False, # sampledict.ymin, sampledict.ymax) data_range = cls.getDataRange( - filegeodict, sampledict, - first_column_duplicated=first_column_duplicated) + filegeodict, sampledict, first_column_duplicated=first_column_duplicated + ) data, geodict = cls.readFile(filename, data_range) - if (doPadding or resample) and 'float' not in str(data.dtype): + if (doPadding or resample) and "float" not in str(data.dtype): data = data.astype(np.float32) # parent static method - pad_dict = cls.getPadding(filegeodict, samplegeodict, - doPadding=doPadding) + pad_dict = cls.getPadding(filegeodict, samplegeodict, doPadding=doPadding) data, geodict = cls.padGrid(data, geodict, pad_dict) if doPadding: data[np.isinf(data)] = padValue @@ -135,65 +146,70 @@ def load(cls, filename, samplegeodict=None, resample=False, def getFileGeoDict(cls, filename): """Get the spatial extent, resolution, and shape of grid inside ESRI grid file. - :param filename: - File name of ESRI grid file. - :returns: - - GeoDict object specifying spatial extent, resolution, and - shape of grid inside ESRI grid file. - :raises DataSetException: - When the file contains a grid with more than one band. - When the file geodict is internally inconsistent. + + Args: + filename (str): + File name of ESRI grid file. + Returns: + GeoDict object specifying spatial extent, resolution, and + shape of grid inside ESRI grid file. + Raises: + DataSetException: + When the file contains a grid with more than one band. + When the file geodict is internally inconsistent. """ geodict = {} with rasterio.open(filename) as src: aff = get_affine(src) if aff is None: - fmt = 'Could not find .transform attribute from GDAL dataset.' + fmt = "Could not find .transform attribute from GDAL dataset." raise AttributeError(fmt) - geodict['dx'] = aff.a - geodict['dy'] = -1*aff.e - geodict['xmin'] = aff.xoff + geodict['dx']/2.0 - geodict['ymax'] = aff.yoff - geodict['dy']/2.0 + geodict["dx"] = aff.a + geodict["dy"] = -1 * aff.e + geodict["xmin"] = aff.xoff + geodict["dx"] / 2.0 + geodict["ymax"] = aff.yoff - geodict["dy"] / 2.0 shp = src.shape if len(shp) > 2: - raise DataSetException( - 'Cannot support grids with more than one band') - geodict['ny'] = src.height - geodict['nx'] = src.width - geodict['xmax'] = geodict['xmin'] + (geodict['nx']-1)*geodict['dx'] - if geodict['xmax'] == geodict['xmin']: + raise DataSetException("Cannot support grids with more than one band") + geodict["ny"] = src.height + geodict["nx"] = src.width + geodict["xmax"] = geodict["xmin"] + (geodict["nx"] - 1) * geodict["dx"] + if geodict["xmax"] == geodict["xmin"]: pass - geodict['ymin'] = geodict['ymax'] - (geodict['ny']-1)*geodict['dy'] + geodict["ymin"] = geodict["ymax"] - (geodict["ny"] - 1) * geodict["dy"] gd = GeoDict(geodict) - newgeodict, first_column_duplicated = cls.checkFirstColumnDuplicated( - gd) + newgeodict, first_column_duplicated = cls.checkFirstColumnDuplicated(gd) return (newgeodict, first_column_duplicated) @classmethod def _subsetRegions(self, src, sampledict, fgeodict, firstColumnDuplicated): """Internal method used to do subsampling of data for all three GMT formats. - :param zvar: - A numpy array-like thing (CDF/HDF variable, or actual numpy array) - :param sampledict: - GeoDict object with bounds and row/col information. - :param fgeodict: - GeoDict object with the file information. - :param firstColumnDuplicated: - Boolean - is this a file where the last column of data is the same - as the first (for grids that span entire globe). - :returns: - Tuple of (data,geodict) (subsetted data and geodict describing that - data). + + Args: + sampledict (GeoDict object): + with bounds and row/col information. + fgeodict (GeoDict object): + with the file information. + firstColumnDuplicated (boolean): + is this a file where the last column of data is the same + as the first (for grids that span entire globe). + Returns: + Tuple of (data,geodict) (subsetted data and geodict describing that + data). """ txmin, txmax, tymin, tymax = ( - sampledict.xmin, sampledict.xmax, sampledict.ymin, sampledict.ymax) + sampledict.xmin, + sampledict.xmax, + sampledict.ymin, + sampledict.ymax, + ) # trows, tcols = (sampledict.ny, sampledict.nx) if fgeodict.xmin > fgeodict.xmax: fxmax = fgeodict.xmax + 360 @@ -223,8 +239,8 @@ def _subsetRegions(self, src, sampledict, fgeodict, firstColumnDuplicated): data = np.squeeze(data) if firstColumnDuplicated: data = data[:, 0:-1] - tfdict['xmax'] -= geodict.dx - tfdict['nx'] -= 1 + tfdict["xmax"] -= geodict.dx + tfdict["nx"] -= 1 geodict = GeoDict(tfdict) else: # what are the nearest grid coordinates to our desired bounds? @@ -236,10 +252,10 @@ def _subsetRegions(self, src, sampledict, fgeodict, firstColumnDuplicated): # txmin = 1.5 and # txmax = 4.5. if not fgeodict.isAligned(sampledict): - txmin2 = gxmin + dx*np.floor((txmin - gxmin)/dx) - txmax2 = gxmin + dx*np.ceil((txmax - gxmin)/dx) - tymin2 = gymin + dy*np.floor((tymin - gymin)/dy) - tymax2 = gymin + dy*np.ceil((tymax - gymin)/dy) + txmin2 = gxmin + dx * np.floor((txmin - gxmin) / dx) + txmax2 = gxmin + dx * np.ceil((txmax - gxmin) / dx) + tymin2 = gymin + dy * np.floor((tymin - gymin) / dy) + tymax2 = gymin + dy * np.ceil((tymax - gymin) / dy) else: txmin2 = xmin txmax2 = xmax @@ -263,21 +279,21 @@ def _subsetRegions(self, src, sampledict, fgeodict, firstColumnDuplicated): # tny = (ilry1 - iuly1)+1 # tnx = (ilrx1 - iulx1)+1 + (ilrx2 - iulx2)+1 - window1 = ((iuly1, ilry1+1), (iulx1, ilrx1+1)) - window2 = ((iuly2, ilry2+1), (iulx2, ilrx2+1)) + window1 = ((iuly1, ilry1 + 1), (iulx1, ilrx1 + 1)) + window2 = ((iuly2, ilry2 + 1), (iulx2, ilrx2 + 1)) section1 = src.read(1, window=window1) section2 = src.read(1, window=window2) data = np.hstack((section1, section2)) tfdict = {} newymax, newxmin = fgeodict.getLatLon(iuly1, iulx1) newymin, newxmax = fgeodict.getLatLon(ilry2, ilrx2) - tfdict['xmin'] = newxmin - tfdict['xmax'] = newxmax - tfdict['ymin'] = newymin - tfdict['ymax'] = newymax - tfdict['dx'] = dx - tfdict['dy'] = dy - tfdict['ny'], tfdict['nx'] = data.shape + tfdict["xmin"] = newxmin + tfdict["xmax"] = newxmax + tfdict["ymin"] = newymin + tfdict["ymax"] = newymax + tfdict["dx"] = dx + tfdict["dy"] = dy + tfdict["ny"], tfdict["nx"] = data.shape geodict = GeoDict(tfdict) else: iuly, iulx = fgeodict.getRowCol(tymax2, txmin2) @@ -285,20 +301,20 @@ def _subsetRegions(self, src, sampledict, fgeodict, firstColumnDuplicated): # tny = (ilry - iuly)+1 # tnx = (ilrx - iulx)+1 - window = ((iuly, ilry+1), (iulx, ilrx+1)) + window = ((iuly, ilry + 1), (iulx, ilrx + 1)) tfdict = {} newymax, newxmin = fgeodict.getLatLon(iuly, iulx) newymin, newxmax = fgeodict.getLatLon(ilry, ilrx) - tfdict['xmin'] = newxmin - tfdict['xmax'] = newxmax - tfdict['ymin'] = newymin - tfdict['ymax'] = newymax - tfdict['dx'] = dx - tfdict['dy'] = dy + tfdict["xmin"] = newxmin + tfdict["xmax"] = newxmax + tfdict["ymin"] = newymin + tfdict["ymax"] = newymax + tfdict["dx"] = dx + tfdict["dy"] = dy # window = ((iymin,iymax+1),(ixmin,ixmax+1)) data = src.read(1, window=window) data = np.squeeze(data) - tfdict['ny'], tfdict['nx'] = data.shape + tfdict["ny"], tfdict["nx"] = data.shape geodict = GeoDict(tfdict) return (data, geodict) @@ -308,33 +324,35 @@ def readFile(cls, filename, data_range): """ Read an ESRI flt/bip/bil/bsq formatted file using rasterIO (GDAL Python wrapper). - :param filename: - Input ESRI formatted grid file. - :param data_range: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. - :returns: - A tuple of (data,geodict) where data is a 2D numpy array of all - data found inside bounds, and - geodict gives the geo-referencing information for the data. + + Args: + filename (str): + Input ESRI formatted grid file. + data_range (dict): + containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. + Returns: + A tuple of (data,geodict) where data is a 2D numpy array of all + data found inside bounds, and + geodict gives the geo-referencing information for the data. """ - iulx1 = data_range['iulx1'] - ilrx1 = data_range['ilrx1'] - iuly1 = data_range['iuly1'] - ilry1 = data_range['ilry1'] - if 'iulx2' in data_range: - iulx2 = data_range['iulx2'] - ilrx2 = data_range['ilrx2'] - iuly2 = data_range['iuly2'] - ilry2 = data_range['ilry2'] + iulx1 = data_range["iulx1"] + ilrx1 = data_range["ilrx1"] + iuly1 = data_range["iuly1"] + ilry1 = data_range["ilry1"] + if "iulx2" in data_range: + iulx2 = data_range["iulx2"] + ilrx2 = data_range["ilrx2"] + iuly2 = data_range["iuly2"] + ilry2 = data_range["ilry2"] else: iulx2 = None ilrx2 = None @@ -343,12 +361,10 @@ def readFile(cls, filename, data_range): data = None with rasterio.open(filename) as src: - window1 = ((iuly1, ilry1), - (iulx1, ilrx1)) + window1 = ((iuly1, ilry1), (iulx1, ilrx1)) section1 = src.read(1, window=window1) - if 'iulx2' in data_range: - window2 = ((iuly2, ilry2), - (iulx2, ilrx2)) + if "iulx2" in data_range: + window2 = ((iuly2, ilry2), (iulx2, ilrx2)) section2 = src.read(1, window=window2) data = np.hstack((section1, section2)) else: @@ -364,75 +380,81 @@ def readFile(cls, filename, data_range): ny, nx = data.shape filegeodict, first_column_duplicated = cls.getFileGeoDict(filename) ymax1, xmin1 = filegeodict.getLatLon(iuly1, iulx1) - ymin1, xmax1 = filegeodict.getLatLon(ilry1-1, ilrx1-1) + ymin1, xmax1 = filegeodict.getLatLon(ilry1 - 1, ilrx1 - 1) xmin = xmin1 ymax = ymax1 ymin = ymin1 xmax = xmax1 if iulx2 is not None: - ymin2, xmax2 = filegeodict.getLatLon(ilry2-1, ilrx2-1) + ymin2, xmax2 = filegeodict.getLatLon(ilry2 - 1, ilrx2 - 1) xmax = xmax2 dx = filegeodict.dx dy = filegeodict.dy - geodict = GeoDict({'xmin': xmin, - 'xmax': xmax, - 'ymin': ymin, - 'ymax': ymax, - 'nx': nx, - 'ny': ny, - 'dx': dx, - 'dy': dy}, adjust='res') + geodict = GeoDict( + { + "xmin": xmin, + "xmax": xmax, + "ymin": ymin, + "ymax": ymax, + "nx": nx, + "ny": ny, + "dx": dx, + "dy": dy, + }, + adjust="res", + ) return (data, geodict) def _getHeader(self): hdr = {} - if sys.byteorder == 'little': - hdr['BYTEORDER'] = 'LSBFIRST' + if sys.byteorder == "little": + hdr["BYTEORDER"] = "LSBFIRST" else: - hdr['BYTEORDER'] = 'MSBFIRST' - hdr['LAYOUT'] = 'BIL' - hdr['NROWS'], hdr['NCOLS'] = self._data.shape - hdr['NBANDS'] = 1 + hdr["BYTEORDER"] = "MSBFIRST" + hdr["LAYOUT"] = "BIL" + hdr["NROWS"], hdr["NCOLS"] = self._data.shape + hdr["NBANDS"] = 1 if self._data.dtype == np.uint8: - hdr['NBITS'] = 8 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 8 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif self._data.dtype == np.int8: - hdr['NBITS'] = 8 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 8 + hdr["PIXELTYPE"] = "SIGNEDINT" elif self._data.dtype == np.uint16: - hdr['NBITS'] = 16 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 16 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif self._data.dtype == np.int16: - hdr['NBITS'] = 16 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 16 + hdr["PIXELTYPE"] = "SIGNEDINT" elif self._data.dtype == np.uint32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif self._data.dtype == np.int32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "SIGNEDINT" elif self._data.dtype == np.float32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'FLOAT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "FLOAT" elif self._data.dtype == np.float64: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'FLOAT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "FLOAT" else: raise DataSetException( - 'Data type "%s" not supported.' % str(self._data.dtype)) - hdr['BANDROWBYTES'] = hdr['NCOLS']*(hdr['NBITS']/8) - hdr['TOTALROWBYTES'] = hdr['NCOLS']*(hdr['NBITS']/8) - hdr['ULXMAP'] = self._geodict.xmin - hdr['ULYMAP'] = self._geodict.ymax - hdr['XDIM'] = self._geodict.dx - hdr['YDIM'] = self._geodict.dy + 'Data type "%s" not supported.' % str(self._data.dtype) + ) + hdr["BANDROWBYTES"] = hdr["NCOLS"] * (hdr["NBITS"] / 8) + hdr["TOTALROWBYTES"] = hdr["NCOLS"] * (hdr["NBITS"] / 8) + hdr["ULXMAP"] = self._geodict.xmin + hdr["ULYMAP"] = self._geodict.ymax + hdr["XDIM"] = self._geodict.dx + hdr["YDIM"] = self._geodict.dy # try to have a nice readable NODATA value in the header file zmin = np.nanmin(self._data) zmax = np.nanmax(self._data) if self._data.dtype in [np.int8, np.int16, np.int32]: - nodata = np.array([-1*int('9'*i) for i in range(3, 20)]) + nodata = np.array([-1 * int("9" * i) for i in range(3, 20)]) if zmin > nodata[-1]: NODATA = nodata[np.where(nodata < zmin)[0][0]] else: @@ -440,42 +462,58 @@ def _getHeader(self): # smallest NODATA = zmin - 1 else: - nodata = np.array([int('9'*i) for i in range(3, 20)]) + nodata = np.array([int("9" * i) for i in range(3, 20)]) if zmin < nodata[-1]: NODATA = nodata[np.where(nodata > zmin)[0][0]] else: # otherwise just pick an arbitrary value smaller than our # smallest NODATA = zmax + 1 - hdr['NODATA'] = NODATA - keys = ['BYTEORDER', 'LAYOUT', 'NROWS', 'NCOLS', 'NBANDS', 'NBITS', - 'BANDROWBYTES', 'TOTALROWBYTES', 'PIXELTYPE', - 'ULXMAP', 'ULYMAP', 'XDIM', 'YDIM', 'NODATA'] + hdr["NODATA"] = NODATA + keys = [ + "BYTEORDER", + "LAYOUT", + "NROWS", + "NCOLS", + "NBANDS", + "NBITS", + "BANDROWBYTES", + "TOTALROWBYTES", + "PIXELTYPE", + "ULXMAP", + "ULYMAP", + "XDIM", + "YDIM", + "NODATA", + ] hdr2 = OrderedDict() for key in keys: hdr2[key] = hdr[key] return hdr2 - def save(self, filename, format='EHdr'): + def save(self, filename, format="EHdr"): """ Save the data contained in this grid to a float or integer ESRI grid file. Described here: http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=BIL,_BIP,_and_BSQ_raster_files http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/conversion_tools/float_to_raster_conversion_.htm. - :param filename: - String representing file to which data should be saved. - :param format: - Currently this code only supports the GDAL format 'EHdr' (see - formats above.) As rasterIO write support is expanded, this - code should add functionality accordingly. - :raises DataSetException: - When format is not 'EHdr'. + Args: + filename (str): + String representing file to which data should be saved. + format (:obj:`str`:'EHdr', optional): + Currently this code only supports the GDAL format 'EHdr' (see + formats above.) As rasterIO write support is expanded, this + code should add functionality accordingly. + Raises: + DataSetException: + When format is not 'EHdr'. """ - supported = ['EHdr'] + supported = ["EHdr"] if format not in supported: raise DataSetException( - 'Only "%s" file formats supported for saving' % str(supported)) + 'Only "%s" file formats supported for saving' % str(supported) + ) hdr = self._getHeader() # create a reference to the data - this may be overridden by a # downcasted version for doubles @@ -483,19 +521,23 @@ def save(self, filename, format='EHdr'): if self._data.dtype == np.float32: # so we can find/reset nan values without screwing up original data data = self._data.astype(np.float32) - data[np.isnan(data)] = hdr['NODATA'] + data[np.isnan(data)] = hdr["NODATA"] elif self._data.dtype == np.float64: data = self._data.astype(np.float32) - data[np.isnan(data)] = hdr['NODATA'] - warnings.warn(DataSetWarning('Down-casting double precision ' - 'floating point to single precision')) + data[np.isnan(data)] = hdr["NODATA"] + warnings.warn( + DataSetWarning( + "Down-casting double precision " + "floating point to single precision" + ) + ) data.tofile(filename) # write out the header file basefile, ext = os.path.splitext(filename) - hdrfile = basefile+'.hdr' - f = open(hdrfile, 'wt') + hdrfile = basefile + ".hdr" + f = open(hdrfile, "wt") for (key, value) in hdr.items(): value = hdr[key] - f.write('%s %s\n' % (key, str(value))) + f.write("%s %s\n" % (key, str(value))) f.close() diff --git a/mapio/geodict.py b/mapio/geodict.py index dae4131..a52fe41 100755 --- a/mapio/geodict.py +++ b/mapio/geodict.py @@ -29,26 +29,33 @@ def __init__(self, geodict, adjust="bounds"): An object which represents the spatial information for a grid and is guaranteed to be self-consistent. - :param geodict: - A dictionary containing at least some of the following fields: - - xmin Longitude minimum (decimal degrees) (Center of upper left cell) - - xmax Longitude maximum (decimal degrees) (Center of upper right cell) - - ymin Longitude minimum (decimal degrees) (Center of lower left cell) - - ymax Longitude maximum (decimal degrees) (Center of lower right cell) - - dx Cell width (decimal degrees) - - dy Cell height (decimal degrees) - - ny Number of rows of input data (must match input data dimensions) - - nx Number of columns of input data (must match input data dimensions). - - projection proj4 string defining projection. Optional. If not specified, - assumed to be - "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" (WGS-84 geographic). - - :param adjust: - String (one of 'bounds','res') - 'bounds': dx/dy, nx/ny, xmin/ymax are assumed to be correct, xmax/ymin - will be recalculated. - 'res': nx/ny, xmin/ymax, xmax/ymin and assumed to be correct, dx/dy will - be recalculated. + Args: + geodict (GeoDict): + contains at least some of the following fields: + - xmin Longitude minimum (decimal degrees) (Center of upper left + cell) + - xmax Longitude maximum (decimal degrees) (Center of upper right + cell) + - ymin Longitude minimum (decimal degrees) (Center of lower left + cell) + - ymax Longitude maximum (decimal degrees) (Center of lower right + cell) + - dx Cell width (decimal degrees) + - dy Cell height (decimal degrees) + - ny Number of rows of input data (must match input data dimensions) + - nx Number of columns of input data (must match input data + dimensions). + - projection proj4 string defining projection. Optional. If not + specified, + assumed to be + "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" + (WGS-84 geographic). + adjust (:obj:`str`, optional): + String (one of 'bounds','res') + 'bounds': dx/dy, nx/ny, xmin/ymax are assumed to be correct, + xmax/ymin will be recalculated. + 'res': nx/ny, xmin/ymax, xmax/ymin and assumed to be correct, + dx/dy will be recalculated. """ for key in self.REQ_KEYS: if key not in geodict.keys(): @@ -82,16 +89,23 @@ def __init__(self, geodict, adjust="bounds"): def createDictFromBox(cls, xmin, xmax, ymin, ymax, dx, dy, inside=False): """Create GeoDict from two corner points and an x/y resolution. - :param xmin: X coordinate of center of upper left pixel. - :param xmax: X coordinate of center of lower right pixel. - :param ymin: Y coordinate of center of lower right pixel. - :param ymax: Y coordinate of center of upper left pixel. - :param dx: Width of pixels. - :param dy: Height of pixels. - :param inside: - Boolean, indicating whether to ensure that resulting GeoDict - falls inside or outside the input bounds. - :returns: + Args: + xmin (float): + X coordinate of center of upper left pixel. + xmax (float): + X coordinate of center of lower right pixel. + ymin(float): + Y coordinate of center of lower right pixel. + ymax (float): + Y coordinate of center of upper left pixel. + dx (float): + Width of pixels. + dy (float): + Height of pixels. + inside (bool): + indicating whether to ensure that resulting GeoDict + falls inside or outside the input bounds. + Returns: GeoDict object. """ if xmin > xmax: @@ -123,14 +137,21 @@ def createDictFromBox(cls, xmin, xmax, ymin, ymax, dx, dy, inside=False): def createDictFromCenter(cls, cx, cy, dx, dy, xspan, yspan): """Create GeoDict from a center point, dx/dy and a width and height. - :param cx: X coordinate of center point. - :param cy: Y coordinate of center point. - :param dx: Width of pixels. - :param dy: Height of pixels. - :param xspan: Width of desired box. - :param yspan: Height of desired box. - :returns: - GeoDict object. + Args: + cx (float): + X coordinate of center point. + cy (float): + Y coordinate of center point. + dx (float): + Width of pixels. + dy (float): + Height of pixels. + xspan (float): + Width of desired box. + yspan (float): + Height of desired box. + Returns: + GeoDict object. """ xmin = cx - xspan / 2.0 xmax = cx + xspan / 2.0 @@ -141,30 +162,32 @@ def createDictFromCenter(cls, cx, cy, dx, dy, xspan, yspan): def setProjection(self, projection): """Set a new projection for the GeoDict. - :param projection: - Valid proj4 string. - :raises DataSetException: - When input is not valid proj4. + Args: + projection (str): + Valid proj4 string. + Raises: + DataSetException: + When input is not valid proj4. """ try: pyproj.Proj(projection) except BaseException: - raise DataSetException( - "%s is not a valid proj4 string." % geodict["projection"] - ) + raise DataSetException("%s is not a valid proj4 string." % projection) self._projection = projection def getAligned(self, geodict, inside=False): """Return a geodict that is close to the input geodict bounds but aligned with this GeoDict. - :param geodict: - Input GeoDict object, whose bounds will be used, but dx/dy and nx/ny ignored. - :param inside: - Boolean indicating whether the aligned geodict should be inside or outside - input geodict. - :returns: - GeoDict object which is guaranteed to be grid-aligned with this GeoDict. + Args: + geodict (GeoDict): + Input GeoDict object, whose bounds will be used, but dx/dy and nx/ny + ignored. + inside (bool): + indicating whether the aligned geodict should be inside or + outside input geodict. + Returns: + GeoDict object which is guaranteed to be grid-aligned with this GeoDict. """ dx = self.dx dy = self.dy @@ -225,13 +248,15 @@ def getIntersection(self, geodict): """Return a geodict defining intersected area, retaining resolution of the input geodict. - :param geodict: - Input GeoDict object, which should intersect with this GeoDict. - :returns: - GeoDict which represents the intersected area, and is aligned with the input - geodict. - :raises DataSetException: - When input geodict does not intersect at all with this GeoDict. + Args: + geodict (GeoDict): + Input GeoDict object, which should intersect with this GeoDict. + Returns: + GeoDict which represents the intersected area, and is aligned with the input + geodict. + Raises: + DataSetException: + When input geodict does not intersect at all with this GeoDict. """ # return modified input geodict where bounds have been adjusted to be inside @@ -293,12 +318,14 @@ def getBoundsWithin(self, geodict): box aligned with enclosing GeoDict that is guaranteed to be inside input GeoDict. - :param geodict: - GeoDict object that output GeoDict will be contained by. - :raises DataSetException: - When input geodict is not fully contained by this GeoDict, or - if the output GeoDict cannot be aligned with this GeoDict (this - shouldn't happen). + Args: + geodict (GeoDict): + GeoDict object that output GeoDict will be contained by. + Raises: + DataSetException: + When input geodict is not fully contained by this GeoDict, or + if the output GeoDict cannot be aligned with this GeoDict (this + shouldn't happen). """ if not self.contains(geodict): raise DataSetException( @@ -425,10 +452,11 @@ def getBoundsWithin(self, geodict): def isAligned(self, geodict): """Determines if input geodict is grid-aligned with this GeoDict. - :param geodict: - Input geodict whose cells must all be grid-aligned with this GeoDict. - :returns: - True when geodict is grid-aligned, and False if not. + Args: + geodict (GeoDict): + Input geodict whose cells must all be grid-aligned with this GeoDict. + Returns: + True when geodict is grid-aligned, and False if not. """ dx1 = self._dx dx2 = geodict.dx @@ -461,11 +489,12 @@ def isAligned(self, geodict): def doesNotContain(self, geodict): """Determine if input geodict is completely outside this GeoDict. - :param geodict: - Input GeoDict object. - :returns: - True if input geodict is completely outside this GeoDict, - False if not. + Args: + geodict (GeoDict): + Input GeoDict object. + Returns: + True if input geodict is completely outside this GeoDict, + False if not. """ a, b = (self.xmin, self.ymax) e, f = (geodict.xmin, geodict.ymax) @@ -482,11 +511,12 @@ def doesNotContain(self, geodict): def intersects(self, geodict): """Determine if input geodict intersects this GeoDict. - :param geodict: - Input GeoDict object. - :returns: - True if input geodict intersects with this GeoDict, - False if not. + Args: + geodict (GeoDict): + Input GeoDict object. + Returns: + True if input geodict intersects with this GeoDict, + False if not. """ try: xmin = self.xmin @@ -504,37 +534,16 @@ def intersects(self, geodict): input_inside_self = fymin <= ymin and fymax >= ymin if self_inside_input or input_inside_self: return True - # if self.xmin > self.xmax: - # c, d = (self.xmax + 360, self.ymin) - # else: - # c, d = (self.xmax, self.ymin) - - # a, b = (self.xmin, self.ymax) - # e, f = (geodict.xmin, geodict.ymax) - - # if geodict.xmin == 180 and geodict.xmax < 0: - # e, f = (-geodict.xmin, geodict.ymax) - # else: - # e, f = (geodict.xmin, geodict.ymax) - - # g, h = (geodict.xmax, geodict.ymin) - - # inside_x = (e >= a and e < c) or (a >= e and a < c) - # inside_y = (h >= d and h < b) or (d >= h and d < f) - # # inside_x = geodict.xmin >= self._xmin and geodict.xmax <= self._xmax - # # inside_y = geodict.ymin >= self._ymin and geodict.ymax <= self._ymax - # if inside_x and inside_y: - # return True - # return False def contains(self, geodict): """Determine if input geodict is completely outside this GeoDict. - :param geodict: - Input GeoDict object. - :returns: - True if input geodict is completely outside this GeoDict, - False if not. + Args: + geodict (GeoDict): + Input GeoDict object. + Returns: + True if input geodict is completely outside this GeoDict, + False if not. """ inside_x = False if np.abs((self.xmax + self.dx) - (self.xmin + 360)) < self.dx * 0.01: @@ -559,8 +568,9 @@ def contains(self, geodict): def asDict(self): """Return GeoDict object in dictionary representation. - :returns: - Dictionary containing the same fields as found in constructor. + + Returns: + Dictionary containing the same fields as found in constructor. """ tdict = {} tdict["xmin"] = self._xmin @@ -592,9 +602,9 @@ def __repr__(self): def copy(self): """Return an object that is a complete copy of this GeoDict. - :returns: - A GeoDict object whose elements (xmin, xmax, ...) are an exact copy of - this object's elements. + Returns: + A GeoDict object whose elements (xmin, xmax, ...) are an exact copy of + this object's elements. """ if not hasattr(self.__class__, "projection"): self._projection = self.DEFAULT_PROJ4 @@ -614,11 +624,12 @@ def copy(self): def __eq__(self, other): """Check for equality between one GeoDict object and another. - :param other: - Another GeoDict object. - :returns: - True when all GeoDict parameters are no more different than 1e-12, False - otherwise. + Args: + other (GeoDict): + Another GeoDict object. + Returns: + True when all GeoDict parameters are no more different than 1e-12, False + otherwise. """ if np.abs(self._xmin - other._xmin) > self.EPS: return False @@ -642,12 +653,13 @@ 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. + Args: + row (int): + Row dimension index into internal data array. + col (int): + Column dimension index into internal data array. + Returns: + Tuple of latitude and longitude. """ scalar_types = (int, float) if type(row) != type(col): @@ -677,19 +689,20 @@ def getRowCol(self, lat, lon, returnFloat=False, intMethod="round"): """Return data row and column from given geographic coordinates (lat/lon decimal degrees). - :param lat: - Input latitude. - :param lon: - Input longitude. - :param returnFloat: - Boolean indicating whether floating point row/col coordinates should be - returned. - :param intMethod: - String indicating the desired method by which floating point row/col values - should - be converted to integers. Choices are 'round' (default), 'floor', or 'ceil'. - :returns: - Tuple of row and column. + Args: + lat (int,float): + Input latitude. + lon (int,float): + Input longitude. + returnFloat (bool): + indicating whether floating point row/col coordinates should be + returned. + intMethod (:obj:`str`, optional): + String indicating the desired method by which floating point row/col + values should be converted to integers. + Choices are 'round' (default), 'floor', or 'ceil'. + Returns: + Tuple of row and column. """ scalar_types = (int, float) if type(lat) != type(lon): @@ -735,90 +748,103 @@ def getRowCol(self, lat, lon, returnFloat=False, intMethod="round"): @property def xmin(self): """Get xmin value. - :returns: - xmin value. + + Returns: + xmin value. """ return self._xmin @property def xmax(self): """Get xmin value. - :returns: - xmin value. + + Returns: + xmin value. """ return self._xmax @property def ymin(self): """Get xmax value. - :returns: - xmax value. + + Returns: + xmax value. """ return self._ymin @property def ymax(self): """Get xmax value. - :returns: - xmax value. + + Returns: + xmax value. """ return self._ymax @property def dx(self): """Get dx value. - :returns: - dx value. + + Returns: + dx value. """ return self._dx @property def dy(self): """Get dy value. - :returns: - dy value. + + Returns: + dy value. """ return self._dy @property def ny(self): """Get ny value. - :returns: - ny value. + + Returns: + ny value. """ return self._ny @property def nx(self): """Get nx value. - :returns: - nx value. + + Returns: + nx value. """ return self._nx @property def projection(self): """Get projection Proj4 string. - :returns: - Valid Proj4 string. + + Returns: + Valid Proj4 string. """ return self._projection @property def nodata(self): """Get projection Proj4 string. - :returns: - Valid nodata value. + + Returns: + Valid nodata value. """ return self._nodata @xmin.setter def xmin(self, value): """Set xmin value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._xmin = value self.validate() @@ -826,10 +852,13 @@ def xmin(self, value): @xmax.setter def xmax(self, value): """Set xmax value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._xmax = value self.validate() @@ -837,10 +866,13 @@ def xmax(self, value): @ymin.setter def ymin(self, value): """Set ymin value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._ymin = value self.validate() @@ -848,10 +880,13 @@ def ymin(self, value): @ymax.setter def ymax(self, value): """Set ymax value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._ymax = value self.validate() @@ -859,10 +894,13 @@ def ymax(self, value): @dx.setter def dx(self, value): """Set dx value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._dx = value self.validate() @@ -870,10 +908,13 @@ def dx(self, value): @dy.setter def dy(self, value): """Set dy value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._dy = value self.validate() @@ -881,10 +922,13 @@ def dy(self, value): @ny.setter def ny(self, value): """Set ny value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._ny = value self.validate() @@ -892,10 +936,13 @@ def ny(self, value): @nx.setter def nx(self, value): """Set nx value, re-validate object. - :param value: - Value to set. - :raises DataSetException: - When validation fails. + + Args: + value: + Value to set. + Raises: + DataSetException: + When validation fails. """ self._nx = value self.validate() @@ -903,8 +950,9 @@ def nx(self, value): @nodata.setter def nodata(self, value): """Set nodata value. - :param value: - Value to set. + Args: + value: + Value to set. """ self._nodata = value @@ -972,7 +1020,8 @@ def geodict_from_src(src): """Return a geodict from a dataset object. Args: - src (DatasetReader): Open rasterio DatasetReader object. + src (DatasetReader): + Open rasterio DatasetReader object. Returns: GeoDict: GeoDict describing the DatasetReader object. """ @@ -983,6 +1032,19 @@ def geodict_from_src(src): def geodict_from_affine(affine, ny, nx): + """********** + + Args: + affine (): + ********** + ny(): + Number of rows. + nx(): + Number of columns. + Returns: + gd: + ************ + """ geodict = {} geodict["dx"] = affine.a geodict["dy"] = -1 * affine.e @@ -998,6 +1060,12 @@ def geodict_from_affine(affine, ny, nx): def affine_from_geodict(geodict): + """************** + + Args: + geodict (): + ***** + """ xoff = geodict.xmin - geodict.dx / 2.0 yoff = geodict.ymax + geodict.dy / 2.0 xres = geodict.dx @@ -1007,6 +1075,18 @@ def affine_from_geodict(geodict): def get_longitude_intersection(xmin, xmax, fxmin, fxmax): + """************** + + Args: + xmin(): + ***** + xmax(): + ***** + fxmin(): + ***** + fxmax(): + ***** + """ # if we're straddling the 180 meridian, adjust the xmax value so that our # new range of longitudes is 0-360. We'll adjust back at the end if fxmin > fxmax: diff --git a/mapio/gmt.py b/mapio/gmt.py index aaf9263..c6b3c71 100644 --- a/mapio/gmt.py +++ b/mapio/gmt.py @@ -14,7 +14,7 @@ import h5py -'''Grid2D subclass for reading, writing, and manipulating GMT format grids. +"""Grid2D subclass for reading, writing, and manipulating GMT format grids. Usage: @@ -26,76 +26,112 @@ This class supports reading and writing of all three GMT formats: NetCDF, HDF, and the GMT "native" format. -''' +""" -NETCDF_TYPES = {'B': np.uint8, - 'b': np.int8, - 'h': np.int16, - 'i': np.int32, - 'f': np.float32, - 'd': np.float64} +NETCDF_TYPES = { + "B": np.uint8, + "b": np.int8, + "h": np.int16, + "i": np.int32, + "f": np.float32, + "d": np.float64, +} -INVERSE_NETCDF_TYPES = {'uint8': 'B', - 'int8': 'b', - 'int16': 'h', - 'int32': 'i', - 'float32': 'f', - 'float64': 'd'} +INVERSE_NETCDF_TYPES = { + "uint8": "B", + "int8": "b", + "int16": "h", + "int32": "i", + "float32": "f", + "float64": "d", +} def subsetArray(data, data_range, fgeodict): - iulx1 = data_range['iulx1'] - iuly1 = data_range['iuly1'] - ilrx1 = data_range['ilrx1'] - ilry1 = data_range['ilry1'] + """ + Function description + + Args: + data(): + data_range(): + fgeodict(): + Return: + data(): + geodict(): + """ + iulx1 = data_range["iulx1"] + iuly1 = data_range["iuly1"] + ilrx1 = data_range["ilrx1"] + ilry1 = data_range["ilry1"] data1 = data[iuly1:ilry1, iulx1:ilrx1] ymax1, xmin1 = fgeodict.getLatLon(iuly1, iulx1) - ymin1, xmax1 = fgeodict.getLatLon(ilry1-1, ilrx1-1) - if 'iulx2' in data_range: - iulx2 = data_range['iulx2'] - iuly2 = data_range['iuly2'] - ilrx2 = data_range['ilrx2'] - ilry2 = data_range['ilry2'] + ymin1, xmax1 = fgeodict.getLatLon(ilry1 - 1, ilrx1 - 1) + if "iulx2" in data_range: + iulx2 = data_range["iulx2"] + iuly2 = data_range["iuly2"] + ilrx2 = data_range["ilrx2"] + ilry2 = data_range["ilry2"] data2 = data[iuly2:ilry2, iulx2:ilrx2] data = np.hstack((data1, data2)).copy() ymax2, xmin2 = fgeodict.getLatLon(iuly2, iulx2) - ymin2, xmax2 = fgeodict.getLatLon(ilry2-1, ilrx2-1) + ymin2, xmax2 = fgeodict.getLatLon(ilry2 - 1, ilrx2 - 1) ny, nx = data.shape - ymin2, xmax2 = fgeodict.getLatLon(ilry2-1, ilrx2-1) + ymin2, xmax2 = fgeodict.getLatLon(ilry2 - 1, ilrx2 - 1) xmax = xmax2 else: data = data1.copy() xmax = xmax1 ny, nx = data.shape - geodict = GeoDict({'xmin': xmin1, - 'xmax': xmax, - 'ymin': ymin1, - 'ymax': ymax1, - 'dx': fgeodict.dx, - 'dy': fgeodict.dy, - 'nx': nx, - 'ny': ny}) + geodict = GeoDict( + { + "xmin": xmin1, + "xmax": xmax, + "ymin": ymin1, + "ymax": ymax1, + "dx": fgeodict.dx, + "dy": fgeodict.dy, + "nx": nx, + "ny": ny, + } + ) return (data, geodict) def sub2ind(shape, subtpl): """ Convert 2D subscripts into 1D index. - @param shape: Tuple indicating size of 2D array. - @param subtpl: Tuple of (possibly) numpy arrays of row,col values. - @return: 1D array of indices. + + Args: + shape(Tuple): indicating size of 2D array. + subtpl(Tuple): of (possibly) numpy arrays of row,col values. + Returns: + 1D array of indices. """ if len(shape) != 2 or len(shape) != len(subtpl): - raise IndexError("Input size and subscripts must have length 2 and " - "be equal in length") + raise IndexError( + "Input size and subscripts must have length 2 and " "be equal in length" + ) row, col = subtpl ny, nx = shape - ind = nx*row + col + ind = nx * row + col return ind def indexArray(array, shp, i1, i2, j1, j2): + """ + Function description + + Args: + array(): + shp(): + i1(): + i2(): + j1(): + j2(): + Return: + data(): + """ if not isinstance(j1, int) and len(j1): j1 = j1[0] if not isinstance(j2, int) and len(j2): @@ -105,9 +141,9 @@ def indexArray(array, shp, i1, i2, j1, j2): if not isinstance(i2, int) and len(i2): i2 = i2[0] if len(array.shape) == 1: - ny = i2-i1 - nx = j2-j1 - if hasattr(array, 'dtype'): + ny = i2 - i1 + nx = j2 - j1 + if hasattr(array, "dtype"): data = np.zeros((ny, nx), dtype=array.dtype) else: typecode = array.typecode() @@ -122,41 +158,54 @@ def indexArray(array, shp, i1, i2, j1, j2): i += 1 else: ny, nx = array.shape - i1r = ny-i1 - i2r = ny-i2 + i1r = ny - i1 + i2r = ny - i2 data = array[i2r:i1r, j1:j2].copy() return data def createSampleXRange(M, N, filename, bounds=None, dx=None, dy=None): + """ + Function description + + Args: + M(): + N(): + filename(): + bounds(): + dx(): + dy(): + Return: + data(): + """ if dx is None: dx = 1.0 if dy is None: dy = 1.0 if bounds is None: xmin = 0.5 - xmax = xmin + (N-1)*dx + xmax = xmin + (N - 1) * dx ymin = 0.5 - ymax = ymin + (M-1)*dy + ymax = ymin + (M - 1) * dy else: xmin, xmax, ymin, ymax = bounds - data = np.arange(0, M*N).reshape(M, N).astype(np.int32) - cdf = netcdf.netcdf_file(filename, 'w') - cdf.createDimension('side', 2) - cdf.createDimension('xysize', M*N) - dim = cdf.createVariable('dimension', 'i', ('side',)) + data = np.arange(0, M * N).reshape(M, N).astype(np.int32) + cdf = netcdf.netcdf_file(filename, "w") + cdf.createDimension("side", 2) + cdf.createDimension("xysize", M * N) + dim = cdf.createVariable("dimension", "i", ("side",)) dim[:] = np.array([N, M]) - spacing = cdf.createVariable('spacing', 'i', ('side',)) + spacing = cdf.createVariable("spacing", "i", ("side",)) spacing[:] = np.array([dx, dy]) zrange = cdf.createVariable( - 'z_range', INVERSE_NETCDF_TYPES[str(data.dtype)], ('side',)) + "z_range", INVERSE_NETCDF_TYPES[str(data.dtype)], ("side",) + ) zrange[:] = np.array([data.min(), data.max()]) - x_range = cdf.createVariable('x_range', 'd', ('side',)) + x_range = cdf.createVariable("x_range", "d", ("side",)) x_range[:] = np.array([xmin, xmax]) - y_range = cdf.createVariable('y_range', 'd', ('side',)) + y_range = cdf.createVariable("y_range", "d", ("side",)) y_range[:] = np.array([ymin, ymax]) - z = cdf.createVariable( - 'z', INVERSE_NETCDF_TYPES[str(data.dtype)], ('xysize',)) + z = cdf.createVariable("z", INVERSE_NETCDF_TYPES[str(data.dtype)], ("xysize",)) z[:] = data.flatten() cdf.close() return data @@ -166,34 +215,38 @@ def createSampleGrid(M, N): """Used for internal testing - create an NxN grid with lower left corner at 0.5,0.5, dx/dy = 1.0. - :param M: - Number of rows in output grid - :param N: - Number of columns in output grid - :returns: - GMTGrid object where data values are an NxN array of values from 0 - to N-squared minus 1, and geodict lower left corner is at 0.5/0.5 - and cell dimensions are 1.0. + Args: + M (int): + Number of rows in output grid + N (int): + Number of columns in output grid + Returns: + GMTGrid(class): + object where data values are an NxN array of values from 0 + to N-squared minus 1, and geodict lower left corner is at 0.5/0.5 + and cell dimensions are 1.0. """ - data = np.arange(0, M*N).reshape(M, N) + data = np.arange(0, M * N).reshape(M, N) # arange gives int64 by default, not supported by netcdf3 data = data.astype(np.int32) - xvar = np.arange(0.5, 0.5+N, 1.0) - yvar = np.arange(0.5, 0.5+M, 1.0) - geodict = {'ny': N, - 'nx': N, - 'xmin': 0.5, - 'xmax': xvar[-1], - 'ymin': 0.5, - 'ymax': yvar[-1], - 'dx': 1.0, - 'dy': 1.0} + xvar = np.arange(0.5, 0.5 + N, 1.0) + yvar = np.arange(0.5, 0.5 + M, 1.0) + geodict = { + "ny": N, + "nx": N, + "xmin": 0.5, + "xmax": xvar[-1], + "ymin": 0.5, + "ymax": yvar[-1], + "dx": 1.0, + "dy": 1.0, + } gmtgrid = GMTGrid(data, geodict) return gmtgrid class GMTGrid(Grid2D): - '''Grid2D subclass for reading,writing, and manipulating GMT format grids. + """Grid2D subclass for reading,writing, and manipulating GMT format grids. Usage: @@ -205,66 +258,70 @@ class GMTGrid(Grid2D): This class supports reading and writing of all three GMT formats: NetCDF, HDF, and the GMT "native" format. - ''' + """ def __init__(self, data, geodict): """Construct a GMTGrid object. - :param data: - 2D numpy data array (must match geodict spec) - :param geodict: - Dictionary specifying the spatial extent,resolution and shape of + + + Args: + data(): + 2D numpy data array (must match geodict spec) + geodict(): + Dictionary specifying the spatial extent,resolution and shape of the data. - :returns: - A GMTGrid object. - :raises DataSetException: - When data and geodict dimensions do not match. + Returns: + A GMTGrid object. + Raises: + DataSetException: + When data and geodict dimensions do not match. """ m, n = data.shape if m != geodict.ny or n != geodict.nx: - raise DataSetException( - 'Input geodict does not match shape of input data.') + raise DataSetException("Input geodict does not match shape of input data.") self._data = data self._geodict = geodict @classmethod def getFileType(cls, grdfile): """Get the GMT sub-format (netcdf, hdf, or GMT binary). - :param grdfile: - File name of suspected GMT grid file. - :returns: - One of 'netcdf' (NetCDF version 3), 'hdf' (NetCDF version 4), - 'native' (so-called GMT native format), or 'unknown'. + Args: + grdfile(): + File name of suspected GMT grid file. + Returns: + One of 'netcdf' (NetCDF version 3), 'hdf' (NetCDF version 4), + 'native' (so-called GMT native format), or 'unknown'. """ isThree = True if sys.version_info.major == 2: isThree = False - f = open(grdfile, 'rb') + f = open(grdfile, "rb") # check to see if it's HDF or CDF - ftype = 'unknown' + ftype = "unknown" try: f.seek(1, 0) if isThree: - hdfsig = f.read(3).decode('utf-8') + hdfsig = f.read(3).decode("utf-8") else: - hdfsig = ''.join(struct.unpack('ccc', f.read(3))) + hdfsig = "".join(struct.unpack("ccc", f.read(3))) - if hdfsig == 'HDF': - ftype = 'hdf' + if hdfsig == "HDF": + ftype = "hdf" else: f.seek(0, 0) if isThree: - cdfsig = f.read(3).decode('utf-8') + cdfsig = f.read(3).decode("utf-8") else: - cdfsig = ''.join(struct.unpack('ccc', f.read(3))) - if cdfsig == 'CDF': - ftype = 'netcdf' + cdfsig = "".join(struct.unpack("ccc", f.read(3))) + if cdfsig == "CDF": + ftype = "netcdf" else: f.seek(0, 0) - nx = struct.unpack('I', f.read(4))[0] + nx = struct.unpack("I", f.read(4))[0] f.seek(8, 0) - offset = struct.unpack('I', f.read(4))[0] + offset = struct.unpack("I", f.read(4))[0] if (offset == 0 or offset == 1) and nx > 0: - ftype = 'native' + ftype = "native" except Exception: pass f.close() @@ -275,42 +332,42 @@ def getFileGeoDict(cls, filename): """Get the spatial extent, resolution, and shape of grid inside NetCDF file. - :param filename: - File name of NetCDF file. - :returns: - GeoDict specifying spatial extent, resolution, and shape of grid - inside NetCDF file. - :raises DataSetException: - When the file is not detectable as one of the three flavors of - GMT grids. + Args: + filename(): + File name of NetCDF file. + Returns: + GeoDict specifying spatial extent, resolution, and shape of grid + inside NetCDF file. + Raises: + DataSetException: + When the file is not detectable as one of the three flavors of + GMT grids. """ ftype = cls.getFileType(filename) - if ftype == 'native': - geodict, xvar, yvar, fmt, zscale, zoffset = cls.getNativeHeader( - filename) - elif ftype == 'netcdf': + if ftype == "native": + geodict, xvar, yvar, fmt, zscale, zoffset = cls.getNativeHeader(filename) + elif ftype == "netcdf": geodict, xvar, yvar = cls.getNetCDFHeader(filename) - elif ftype == 'hdf': + elif ftype == "hdf": geodict, xvar, yvar = cls.getHDFHeader(filename) else: - raise DataSetException( - 'Unknown file type for file "%s".' % filename) - newgeodict, first_column_duplicated = cls.checkFirstColumnDuplicated( - geodict) + raise DataSetException('Unknown file type for file "%s".' % filename) + newgeodict, first_column_duplicated = cls.checkFirstColumnDuplicated(geodict) return (newgeodict, first_column_duplicated) @classmethod def getNativeHeader(cls, fname, fmt=None): """Get the header information from a GMT native grid file. - :param fname: - File name of GMT native grid - :param fmt: - One of: - - 'h' for 16 bit signed integer - - 'i' for 32 bit signed integer - - 'f' for 32 bit float - - 'd' for 64 bit float - :returns: + Args: + fname: + File name of GMT native grid + fmt: + One of: + - 'h' for 16 bit signed integer + - 'i' for 32 bit signed integer + - 'f' for 32 bit float + - 'd' for 64 bit float + Returns: - GeoDict specifying spatial extent, resolution, and shape of grid inside NetCDF file. - xvar array specifying X coordinates of data columns @@ -320,7 +377,7 @@ def getNativeHeader(cls, fname, fmt=None): - zscale Data multiplier - zoffset Value to be added to data """ - f = open(fname, 'rb') + f = open(fname, "rb") # Given that we can't automatically distinguish between 32 bit ints # and 32 bit floats, we'll use this value as a cutoff for "detecting" # when a value read from disk as a float has an "unreasonably" high @@ -328,61 +385,61 @@ def getNativeHeader(cls, fname, fmt=None): # keyword if you want to be sure. MAX_FLOAT_EXP = 30 fsize = os.path.getsize(fname) - datalen = fsize-892 + datalen = fsize - 892 f.seek(0, 0) geodict = {} - geodict['nx'] = struct.unpack('I', f.read(4))[0] - geodict['ny'] = struct.unpack('I', f.read(4))[0] - offset = struct.unpack('I', f.read(4))[0] - geodict['xmin'] = struct.unpack('d', f.read(8))[0] - geodict['xmax'] = struct.unpack('d', f.read(8))[0] + geodict["nx"] = struct.unpack("I", f.read(4))[0] + geodict["ny"] = struct.unpack("I", f.read(4))[0] + offset = struct.unpack("I", f.read(4))[0] + geodict["xmin"] = struct.unpack("d", f.read(8))[0] + geodict["xmax"] = struct.unpack("d", f.read(8))[0] # GMT grids can have longitudes from 0-360, we're correcting for # that here - if geodict['xmax'] > 180: - geodict['xmax'] = geodict['xmax'] - 360 - if geodict['xmax'] > 180: - geodict['xmin'] = geodict['xmin'] - 360 - - geodict['ymin'] = struct.unpack('d', f.read(8))[0] - geodict['ymax'] = struct.unpack('d', f.read(8))[0] - zmin = struct.unpack('d', f.read(8))[0] # noqa - zmax = struct.unpack('d', f.read(8))[0] # noqa - geodict['dx'] = struct.unpack('d', f.read(8))[0] - geodict['dy'] = struct.unpack('d', f.read(8))[0] - zscale = struct.unpack('d', f.read(8))[0] - zoffset = struct.unpack('d', f.read(8))[0] + if geodict["xmax"] > 180: + geodict["xmax"] = geodict["xmax"] - 360 + if geodict["xmax"] > 180: + geodict["xmin"] = geodict["xmin"] - 360 + + geodict["ymin"] = struct.unpack("d", f.read(8))[0] + geodict["ymax"] = struct.unpack("d", f.read(8))[0] + zmin = struct.unpack("d", f.read(8))[0] # noqa + zmax = struct.unpack("d", f.read(8))[0] # noqa + geodict["dx"] = struct.unpack("d", f.read(8))[0] + geodict["dy"] = struct.unpack("d", f.read(8))[0] + zscale = struct.unpack("d", f.read(8))[0] + zoffset = struct.unpack("d", f.read(8))[0] xunits = f.read(80).strip() # noqa yunits = f.read(80).strip() # noqa zunits = f.read(80).strip() # noqa title = f.read(80).strip() # noqa command = f.read(320).strip() # noqa remark = f.read(160).strip() # noqa - npixels = geodict['ny']*geodict['nx'] - lenshort = npixels*2 - lenfloat = npixels*4 - lendouble = npixels*8 + npixels = geodict["ny"] * geodict["nx"] + lenshort = npixels * 2 + lenfloat = npixels * 4 + lendouble = npixels * 8 if fmt is None: if datalen == lenshort: - fmt = 'h' + fmt = "h" if datalen == lendouble: - fmt = 'd' + fmt = "d" # let's try to guess whether this is float or int if datalen == lenfloat: fpos = f.tell() # read 1 byte, check to see if it's nan or 0 - if it is, then # we definitely have a float - dbytes = struct.unpack('f', f.read(4))[0] + dbytes = struct.unpack("f", f.read(4))[0] while dbytes == 0.0: - dbytes = struct.unpack('f', f.read(4))[0] + dbytes = struct.unpack("f", f.read(4))[0] f.seek(fpos) # go back to where we were if np.isnan(dbytes): - fmt = 'f' + fmt = "f" # does this have a crazy large exponent? elif int(np.abs(np.log10(dbytes))) > MAX_FLOAT_EXP: - fmt = 'i' + fmt = "i" else: - fmt = 'f' + fmt = "f" f.close() # We are going to represent all grids internally as grid-line @@ -390,14 +447,16 @@ def getNativeHeader(cls, fname, fmt=None): # grids is depicted well here: # http://gmt.soest.hawaii.edu/doc/5.1.0/GMT_Docs.html#grid-registration-the-r-option if offset == 1: - geodict['xmin'] += geodict['dx']/2.0 - geodict['xmax'] -= geodict['dx']/2.0 - geodict['ymin'] += geodict['dy']/2.0 - geodict['ymax'] -= geodict['dy']/2.0 + geodict["xmin"] += geodict["dx"] / 2.0 + geodict["xmax"] -= geodict["dx"] / 2.0 + geodict["ymin"] += geodict["dy"] / 2.0 + geodict["ymax"] -= geodict["dy"] / 2.0 xvar, dx2 = np.linspace( - geodict['xmin'], geodict['xmax'], num=geodict['nx'], retstep=True) + geodict["xmin"], geodict["xmax"], num=geodict["nx"], retstep=True + ) yvar, dy2 = np.linspace( - geodict['ymin'], geodict['ymax'], num=geodict['ny'], retstep=True) + geodict["ymin"], geodict["ymax"], num=geodict["ny"], retstep=True + ) gd = GeoDict(geodict) return (gd, xvar, yvar, fmt, zscale, zoffset) @@ -423,11 +482,11 @@ def readFile(cls, filename, data_range): geodict gives the geo-referencing information for the data. """ ftype = cls.getFileType(filename) - if ftype == 'netcdf': + if ftype == "netcdf": data, geodict = cls.readNetCDF(filename, data_range) - elif ftype == 'hdf': + elif ftype == "hdf": data, geodict = cls.readHDF(filename, data_range) - elif ftype == 'native': + elif ftype == "native": data, geodict = cls.readGMTNative(filename, data_range) return (data, geodict) @@ -435,22 +494,23 @@ def readFile(cls, filename, data_range): def readNetCDF(cls, filename, data_range): """Read a NetCDF formatted GMT file. - :param filename: - Input GMT formatted grid file. - :param data_range: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. - :returns: - data is a 2D numpy array of all data found inside bounds, and - geodict gives the geo-referencing information for the data. + Args: + filename(): + Input GMT formatted grid file. + data_range(): + Dictionary containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. + Returns: + data is a 2D numpy array of all data found inside bounds, and + geodict gives the geo-referencing information for the data. """ fgeodict, first_column_duplicated = cls.getFileGeoDict(filename) # dx = fgeodict.dx @@ -462,10 +522,10 @@ def readNetCDF(cls, filename, data_range): # ilrx1 = data_range['ilrx1'] # ilry1 = data_range['ilry1'] cdf = netcdf.netcdf_file(filename) - if 'z' in cdf.variables: - zvar = cdf.variables['z'] + if "z" in cdf.variables: + zvar = cdf.variables["z"] else: - zvar = cdf.variables['elevation'] + zvar = cdf.variables["elevation"] isScanLine = len(zvar.shape) == 1 if isScanLine: data = indexArray(zvar, (gny, gnx), 0, gny, 0, gnx) @@ -488,50 +548,49 @@ def readGMTNative(cls, fname, data_range, fmt=None): http://gmt.soest.hawaii.edu/doc/5.1.2/GMT_Docs.html#native-binary-grid-files - :param fname: - File name of GMT native grid - :param data_range: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. - :param fmt: - Data width, one of: - - 'i' (16 bit signed integer) - - 'l' (32 bit signed integer) - - 'f' (32 bit float) - - 'd' (64 bit float) - - Strictly speaking, this is only necessary when the data file - is 32 bit float or 32 bit integer, as there is no *sure* way - to tell from the header or data which data type is contained - in the file. If fmt is None, then the code will try to - guess as best it can from the data whether it is integer or - floating point data. Caveat emptor! - - :returns: - Tuple of data (2D numpy array of data, possibly subsetted from - file) and geodict (see above). + Args: + fname (str): + File name of GMT native grid + data_range (dict): + containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. + fmt (str): + Data width, one of: + - 'i' (16 bit signed integer) + - 'l' (32 bit signed integer) + - 'f' (32 bit float) + - 'd' (64 bit float) + + Strictly speaking, this is only necessary when the data file + is 32 bit float or 32 bit integer, as there is no *sure* way + to tell from the header or data which data type is contained + in the file. If fmt is None, then the code will try to + guess as best it can from the data whether it is integer or + floating point data. Caveat emptor! + Returns: + Tuple of data (2D numpy array of data, possibly subsetted from + file) and geodict (see above). """ HDRLEN = 892 - fgeodict, xvar, yvar, fmt, zscale, zoffset = cls.getNativeHeader( - fname, fmt) + fgeodict, xvar, yvar, fmt, zscale, zoffset = cls.getNativeHeader(fname, fmt) t, first_column_duplicated = cls.getFileGeoDict(fname) - sfmt = '%i%s' % (fgeodict.nx*fgeodict.ny, fmt) - dwidths = {'h': 2, 'i': 4, 'f': 4, 'd': 8} + sfmt = "%i%s" % (fgeodict.nx * fgeodict.ny, fmt) + dwidths = {"h": 2, "i": 4, "f": 4, "d": 8} dwidth = dwidths[fmt] - f = open(fname, 'rb') + f = open(fname, "rb") f.seek(HDRLEN) # TODO - do something clever with memory mapping some day so that # we don't always have to read the whole file just to subset it. # Currently not an issue. - dbytes = f.read(fgeodict.nx*fgeodict.ny*dwidth) + dbytes = f.read(fgeodict.nx * fgeodict.ny * dwidth) f.close() data = np.flipud(np.array(struct.unpack(sfmt, dbytes))) @@ -548,26 +607,30 @@ def readGMTNative(cls, fname, data_range, fmt=None): @classmethod def getNetCDFHeader(cls, filename): """Get the header information from a GMT NetCDF3 file. - :param fname: - File name of GMT NetCDF3 grid - :returns: - - GeoDict specifying spatial extent, resolution, and shape of - grid inside NetCDF file. - - xvar array specifying X coordinates of data columns - - xvar array specifying Y coordinates of data rows + Args: + fname(str): + File name of GMT NetCDF3 grid + Returns: + GeoDict(geodict): + GeoDict specifying spatial extent, resolution, and shape of + grid inside NetCDF file. + xvar(array): + specifying X coordinates of data columns + yvar(array): + specifying Y coordinates of data rows """ cdf = netcdf.netcdf_file(filename) geodict = {} xvarname = None - registration = 'gridline' - if hasattr(cdf, 'node_offset') and getattr(cdf, 'node_offset') == 1: - registration = 'pixel' - if 'x' in cdf.variables.keys(): - xvarname = 'x' - yvarname = 'y' - elif 'lon' in cdf.variables.keys(): - xvarname = 'lon' - yvarname = 'lat' + registration = "gridline" + if hasattr(cdf, "node_offset") and getattr(cdf, "node_offset") == 1: + registration = "pixel" + if "x" in cdf.variables.keys(): + xvarname = "x" + yvarname = "y" + elif "lon" in cdf.variables.keys(): + xvarname = "lon" + yvarname = "lat" if xvarname is not None: xvar = cdf.variables[xvarname].data.copy() @@ -580,76 +643,80 @@ def getNetCDFHeader(cls, filename): xvar = xvar - 360 yvar = cdf.variables[yvarname].data.copy() - geodict['nx'] = len(xvar) - geodict['ny'] = len(yvar) - geodict['xmin'] = xvar.min() - geodict['xmax'] = xvar.max() - geodict['ymin'] = yvar.min() - geodict['ymax'] = yvar.max() - newx = np.linspace( - geodict['xmin'], geodict['xmax'], num=geodict['nx']) - newy = np.linspace( - geodict['ymin'], geodict['ymax'], num=geodict['ny']) - geodict['dx'] = newx[1]-newx[0] - geodict['dy'] = newy[1]-newy[0] - elif 'x_range' in cdf.variables.keys(): - geodict['xmin'] = cdf.variables['x_range'].data[0] - geodict['xmax'] = cdf.variables['x_range'].data[1] - geodict['ymin'] = cdf.variables['y_range'].data[0] - geodict['ymax'] = cdf.variables['y_range'].data[1] - geodict['nx'], geodict['ny'] = cdf.variables['dimension'].data + geodict["nx"] = len(xvar) + geodict["ny"] = len(yvar) + geodict["xmin"] = xvar.min() + geodict["xmax"] = xvar.max() + geodict["ymin"] = yvar.min() + geodict["ymax"] = yvar.max() + newx = np.linspace(geodict["xmin"], geodict["xmax"], num=geodict["nx"]) + newy = np.linspace(geodict["ymin"], geodict["ymax"], num=geodict["ny"]) + geodict["dx"] = newx[1] - newx[0] + geodict["dy"] = newy[1] - newy[0] + elif "x_range" in cdf.variables.keys(): + geodict["xmin"] = cdf.variables["x_range"].data[0] + geodict["xmax"] = cdf.variables["x_range"].data[1] + geodict["ymin"] = cdf.variables["y_range"].data[0] + geodict["ymax"] = cdf.variables["y_range"].data[1] + geodict["nx"], geodict["ny"] = cdf.variables["dimension"].data # geodict['dx'],geodict['dy'] = cdf.variables['spacing'].data - xvar = np.linspace( - geodict['xmin'], geodict['xmax'], num=geodict['nx']) - yvar = np.linspace( - geodict['ymin'], geodict['ymax'], num=geodict['ny']) - geodict['dx'] = xvar[1] - xvar[0] - geodict['dy'] = yvar[1] - yvar[0] + xvar = np.linspace(geodict["xmin"], geodict["xmax"], num=geodict["nx"]) + yvar = np.linspace(geodict["ymin"], geodict["ymax"], num=geodict["ny"]) + geodict["dx"] = xvar[1] - xvar[0] + geodict["dy"] = yvar[1] - yvar[0] else: - raise DataSetException('No support for CDF data file with ' - 'variables: %s' % str(cdf.variables.keys())) + raise DataSetException( + "No support for CDF data file with " + "variables: %s" % str(cdf.variables.keys()) + ) # We are going to represent all grids internally as grid-line # registered. The difference between pixel and gridline-registered # grids is depicted well here: # http://gmt.soest.hawaii.edu/doc/5.1.0/GMT_Docs.html#grid-registration-the-r-option - if registration == 'pixel': - geodict['xmin'] += geodict['dx']/2.0 - geodict['xmax'] -= geodict['dx']/2.0 - geodict['ymin'] += geodict['dy']/2.0 - geodict['ymax'] -= geodict['dy']/2.0 + if registration == "pixel": + geodict["xmin"] += geodict["dx"] / 2.0 + geodict["xmax"] -= geodict["dx"] / 2.0 + geodict["ymin"] += geodict["dy"] / 2.0 + geodict["ymax"] -= geodict["dy"] / 2.0 # because dx/dy are not explicitly defined in netcdf headers, here # we'll assume that those values are adjustable, and we'll preserve # the shape and extent. - gd = GeoDict(geodict, adjust='res') + gd = GeoDict(geodict, adjust="res") return (gd, xvar, yvar) @classmethod - def _subsetRegions(self, zvar, sampledict, fgeodict, xvar, yvar, - firstColumnDuplicated): + def _subsetRegions( + self, zvar, sampledict, fgeodict, xvar, yvar, firstColumnDuplicated + ): """Internal method used to do subsampling of data for all three GMT formats. - :param zvar: - A numpy array-like thing (CDF/HDF variable, or actual numpy array) - :param sampledict: - GeoDict indicating the bounds where data should be sampled. - :param fgeodict: - Geo dictionary with the file information. - :param xvar: - Numpy array specifying X coordinates of data columns - :param yvar: - Numpy array specifying Y coordinates of data rows - :param firstColumnDuplicated: - Boolean - is this a file where the last column of data is the same - as the first (for grids that span entire globe). - :returns: - Tuple of (data,geodict) (subsetted data and geodict describing that - data). + Args: + zvar (array): + A numpy array-like thing (CDF/HDF variable, or actual numpy array) + sampledict (GeDict): + GeoDict indicating the bounds where data should be sampled. + fgeodict (GeoDict): + with file information. + xvar (numpy array): + specifying X coordinates of data columns + yvar (numpy array): + specifying Y coordinates of data rows + firstColumnDuplicated (bool): + is this a file where the last column of data is the same + as the first (for grids that span entire globe). + Returns: + Tuple of (data,geodict) (subsetted data and geodict describing that + data). """ isScanLine = len(zvar.shape) == 1 txmin, txmax, tymin, tymax = ( - sampledict.xmin, sampledict.xmax, sampledict.ymin, sampledict.ymax) + sampledict.xmin, + sampledict.xmax, + sampledict.ymin, + sampledict.ymax, + ) # we're not doing anything fancy with the data here, just cutting out # what we need if fgeodict.xmin > fgeodict.xmax: @@ -679,17 +746,17 @@ def _subsetRegions(self, zvar, sampledict, fgeodict, xvar, yvar, data = np.flipud(indexArray(zvar, (gny, gnx), 0, gny, 0, gnx)) if firstColumnDuplicated: data = data[:, 0:-1] - geodict['xmax'] -= geodict['dx'] + geodict["xmax"] -= geodict["dx"] else: # what are the nearest grid coordinates to our desired bounds? # for example, if the grid starts at xmin = 0.5 and xmax = 6.5 # with dx=1.0, and the user wants txmin = 2.0 and txmax of 4.0, # the pixel coordinates (erring on the side of including more # data), would be at txmin = 1.5 and txmax = 4.5. - txmin2 = gxmin + dx*np.floor((txmin - gxmin)/dx) - txmax2 = gxmin + dx*np.ceil((txmax - gxmin)/dx) - tymin2 = gymin + dy*np.floor((tymin - gymin)/dy) - tymax2 = gymin + dy*np.ceil((tymax - gymin)/dy) + txmin2 = gxmin + dx * np.floor((txmin - gxmin) / dx) + txmax2 = gxmin + dx * np.ceil((txmax - gxmin) / dx) + tymin2 = gymin + dy * np.floor((tymin - gymin) / dy) + tymax2 = gymin + dy * np.ceil((tymax - gymin) / dy) if xmin > xmax: # cut user's request into two regions - one from the minimum # to the meridian, then another from the meridian to the @@ -710,34 +777,35 @@ def _subsetRegions(self, zvar, sampledict, fgeodict, xvar, yvar, # in Python 3 the "long" type was integrated into the int type. if sys.version_info.major == 2: - outcols1 = long(ilrx1-iulx1) # noqa - outcols2 = long(ilrx2-iulx2) # noqa - outcols = long(outcols1+outcols2) # noqa - outrows = long(ilry1-iuly1) # noqa + outcols1 = long(ilrx1 - iulx1) # noqa + outcols2 = long(ilrx2 - iulx2) # noqa + outcols = long(outcols1 + outcols2) # noqa + outrows = long(ilry1 - iuly1) # noqa else: - outcols1 = int(ilrx1-iulx1) - outcols2 = int(ilrx2-iulx2) - outcols = int(outcols1+outcols2) - outrows = int(ilry1-iuly1) + outcols1 = int(ilrx1 - iulx1) + outcols2 = int(ilrx2 - iulx2) + outcols = int(outcols1 + outcols2) + outrows = int(ilry1 - iuly1) section1 = indexArray( - zvar, (gny, gnx), iuly1, ilry1+1, iulx1, ilrx1+1) + zvar, (gny, gnx), iuly1, ilry1 + 1, iulx1, ilrx1 + 1 + ) # section1 = zvar[iuly1:ilry1,iulx1:ilrx1].copy() # section2 = zvar[iuly2:ilry2,iulx2:ilrx2].copy() section2 = indexArray( - zvar, (gny, gnx), iuly2, ilry2+1, iulx2, ilrx2+1) + zvar, (gny, gnx), iuly2, ilry2 + 1, iulx2, ilrx2 + 1 + ) if isScanLine: data = np.concatenate((section1, section2), axis=1) else: - data = np.flipud(np.concatenate( - (section1, section2), axis=1)) + data = np.flipud(np.concatenate((section1, section2), axis=1)) outrows, outcols = data.shape newymax, newxmin = fgeodict.getLatLon(iuly1, iulx1) newymin, newxmax = fgeodict.getLatLon(ilry2, ilrx2) - geodict['xmin'] = newxmin - geodict['xmax'] = newxmax - geodict['ymin'] = newymin - geodict['ymax'] = newymax - geodict['ny'], geodict['nx'] = data.shape + geodict["xmin"] = newxmin + geodict["xmax"] = newxmax + geodict["ymin"] = newymin + geodict["ymax"] = newymax + geodict["ny"], geodict["nx"] = data.shape else: iuly, iulx = fgeodict.getRowCol(tymax2, txmin2) ilry, ilrx = fgeodict.getRowCol(tymin2, txmax2) @@ -745,19 +813,19 @@ def _subsetRegions(self, zvar, sampledict, fgeodict, xvar, yvar, # tnx = (ilrx - iulx)+1 if isScanLine: - data = indexArray(zvar, (gny, gnx), iuly, - ilry+1, iulx, ilrx+1) + data = indexArray(zvar, (gny, gnx), iuly, ilry + 1, iulx, ilrx + 1) else: - data = np.flipud(indexArray( - zvar, (gny, gnx), iuly, ilry+1, iulx, ilrx+1)) + data = np.flipud( + indexArray(zvar, (gny, gnx), iuly, ilry + 1, iulx, ilrx + 1) + ) newymax, newxmin = fgeodict.getLatLon(iuly, iulx) newymin, newxmax = fgeodict.getLatLon(ilry, ilrx) - geodict['xmin'] = newxmin - geodict['xmax'] = newxmax - geodict['ymin'] = newymin - geodict['ymax'] = newymax - geodict['ny'], geodict['nx'] = data.shape + geodict["xmin"] = newxmin + geodict["xmax"] = newxmax + geodict["ymin"] = newymin + geodict["ymax"] = newymax + geodict["ny"], geodict["nx"] = data.shape gd = GeoDict(geodict) return (data, gd) @@ -765,23 +833,24 @@ def _subsetRegions(self, zvar, sampledict, fgeodict, xvar, yvar, @classmethod def getHDFHeader(cls, hdffile): """Get the header information from a GMT NetCDF4 (HDF) file. - :param fname: - File name of GMT NetCDF4 grid - :returns: - GeoDict specifying spatial extent, resolution, and shape of grid - inside NetCDF file. + Args: + fname (str): + File name of GMT NetCDF4 grid + Returns: + GeoDict specifying spatial extent, resolution, and shape of grid + inside NetCDF file. """ geodict = {} - f = h5py.File(hdffile, 'r') - registration = 'gridline' - if f.get('node_offset') is not None and f.attrs['node_offset'][0] == 1: - registration = 'pixel' - if 'x' in f.keys(): - xvarname = 'x' - yvarname = 'y' - elif 'lon' in f.keys(): - xvarname = 'lon' - yvarname = 'lat' + f = h5py.File(hdffile, "r") + registration = "gridline" + if f.get("node_offset") is not None and f.attrs["node_offset"][0] == 1: + registration = "pixel" + if "x" in f.keys(): + xvarname = "x" + yvarname = "y" + elif "lon" in f.keys(): + xvarname = "lon" + yvarname = "lat" if xvarname is not None: xvar = f[xvarname][:] yvar = f[yvarname][:] @@ -794,73 +863,69 @@ def getHDFHeader(cls, hdffile): if (xvar > 180.01).any(): xvar = xvar - 360 - geodict['ny'] = len(yvar) - geodict['nx'] = len(xvar) - geodict['xmin'] = xvar[0] - geodict['xmax'] = xvar[-1] - geodict['ymin'] = yvar[0] - geodict['ymax'] = yvar[-1] - newx = np.linspace( - geodict['xmin'], geodict['xmax'], num=geodict['nx']) - newy = np.linspace( - geodict['ymin'], geodict['ymax'], num=geodict['ny']) - geodict['dx'] = newx[1]-newx[0] - geodict['dy'] = newy[1]-newy[0] + geodict["ny"] = len(yvar) + geodict["nx"] = len(xvar) + geodict["xmin"] = xvar[0] + geodict["xmax"] = xvar[-1] + geodict["ymin"] = yvar[0] + geodict["ymax"] = yvar[-1] + newx = np.linspace(geodict["xmin"], geodict["xmax"], num=geodict["nx"]) + newy = np.linspace(geodict["ymin"], geodict["ymax"], num=geodict["ny"]) + geodict["dx"] = newx[1] - newx[0] + geodict["dy"] = newy[1] - newy[0] else: - geodict['xmin'] = f['x_range'][0] - geodict['xmax'] = f['x_range'][1] - geodict['ymin'] = f['y_range'][0] - geodict['ymax'] = f['y_range'][1] - geodict['nx'], geodict['ny'] = ( - f['dimension'][0], f['dimension'][1]) - xvar = np.linspace(geodict['xmin'], geodict['xmax'], - num=geodict['nx']) - yvar = np.linspace(geodict['ymin'], geodict['ymax'], - num=geodict['ny']) - geodict['dx'] = xvar[1] - xvar[0] - geodict['dy'] = yvar[1] - yvar[0] + geodict["xmin"] = f["x_range"][0] + geodict["xmax"] = f["x_range"][1] + geodict["ymin"] = f["y_range"][0] + geodict["ymax"] = f["y_range"][1] + geodict["nx"], geodict["ny"] = (f["dimension"][0], f["dimension"][1]) + xvar = np.linspace(geodict["xmin"], geodict["xmax"], num=geodict["nx"]) + yvar = np.linspace(geodict["ymin"], geodict["ymax"], num=geodict["ny"]) + geodict["dx"] = xvar[1] - xvar[0] + geodict["dy"] = yvar[1] - yvar[0] # We are going to represent all grids internally as grid-line # registered. The difference between pixel and gridline-registered # grids is depicted well here: # http://gmt.soest.hawaii.edu/doc/5.1.0/GMT_Docs.html#grid-registration-the-r-option - if registration == 'pixel': - geodict['xmin'] += geodict['dx']/2.0 - geodict['xmax'] -= geodict['dx']/2.0 - geodict['ymin'] += geodict['dy']/2.0 - geodict['ymax'] -= geodict['dy']/2.0 + if registration == "pixel": + geodict["xmin"] += geodict["dx"] / 2.0 + geodict["xmax"] -= geodict["dx"] / 2.0 + geodict["ymin"] += geodict["dy"] / 2.0 + geodict["ymax"] -= geodict["dy"] / 2.0 f.close() # because dx/dy are not explicitly defined in hdf headers, here we'll # assume that those values are adjustable, and we'll preserve the # shape and extent. - gd = GeoDict(geodict, adjust='res') + gd = GeoDict(geodict, adjust="res") return (gd, xvar, yvar) @classmethod def readHDF(cls, hdffile, data_range): """Read a NetCDF formatted GMT file. - :param filename: - Input GMT formatted grid file. - :param data_range: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. - :returns: - data is a 2D numpy array of all data found inside bounds, and - geodict gives the geo-referencing information for the data. + Args: + hdffile?() + Args: + data_range (dict): + Dictionary containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. + Returns: + data is a 2D numpy array of all data found inside bounds, and + geodict gives the geo-referencing information for the data. """ # fgeodict,xvar,yvar = cls.getHDFHeader(hdffile) fgeodict, first_column_duplicated = cls.getFileGeoDict(hdffile) - f = h5py.File(hdffile, 'r') - zvar = f['z'] + f = h5py.File(hdffile, "r") + zvar = f["z"] data = np.flipud(zvar) if first_column_duplicated: data = data[:, 0:-1] @@ -868,102 +933,109 @@ def readHDF(cls, hdffile, data_range): f.close() return (data, geodict) - def save(self, filename, format='netcdf'): + def save(self, filename, format="netcdf"): """Save a GMTGrid object to a file. - :param filename: - Name of desired output file. - :param format: - One of 'netcdf','hdf' or 'native'. - :raises DataSetException: - When format not one of ('netcdf,'hdf','native') + + Args: + filename (str): + Name of desired output file. + format (:obj:`str`, optional): + One of 'netcdf','hdf' or 'native'. + Raises: + DataSetException: + When format not one of ('netcdf,'hdf','native') """ - if format not in ['netcdf', 'hdf', 'native']: - raise DataSetException('Only NetCDF3, HDF (NetCDF4), and GMT ' - 'native output are supported.') - if format == 'netcdf': - f = netcdf.NetCDFFile(filename, 'w') + if format not in ["netcdf", "hdf", "native"]: + raise DataSetException( + "Only NetCDF3, HDF (NetCDF4), and GMT " "native output are supported." + ) + if format == "netcdf": + f = netcdf.NetCDFFile(filename, "w") m, n = self._data.shape - dx = f.createDimension('x', n) # noqa - dy = f.createDimension('y', m) # noqa - x = f.createVariable('x', np.float64, ('x')) - y = f.createVariable('y', np.float64, ('y')) - x[:] = np.linspace(self._geodict.xmin, - self._geodict.xmax, self._geodict.nx) - y[:] = np.linspace(self._geodict.ymin, - self._geodict.ymax, self._geodict.ny) - z = f.createVariable('z', self._data.dtype, ('y', 'x')) + dx = f.createDimension("x", n) # noqa + dy = f.createDimension("y", m) # noqa + x = f.createVariable("x", np.float64, ("x")) + y = f.createVariable("y", np.float64, ("y")) + x[:] = np.linspace(self._geodict.xmin, self._geodict.xmax, self._geodict.nx) + y[:] = np.linspace(self._geodict.ymin, self._geodict.ymax, self._geodict.ny) + z = f.createVariable("z", self._data.dtype, ("y", "x")) z[:] = np.flipud(self._data) - z.actual_range = np.array( - (np.nanmin(self._data), np.nanmax(self._data))) + z.actual_range = np.array((np.nanmin(self._data), np.nanmax(self._data))) f.close() - elif format == 'hdf': + elif format == "hdf": # Create the file and the top-level attributes GMT wants to see - f = h5py.File(filename, 'w') - f.attrs['Conventions'] = 'COARDS, CF-1.5' - f.attrs['title'] = 'filename' - f.attrs['history'] = 'Created with python ' \ - 'GMTGrid.save(%s,format="hdf")' % filename - f.attrs['GMT_version'] = 'NA' + f = h5py.File(filename, "w") + f.attrs["Conventions"] = "COARDS, CF-1.5" + f.attrs["title"] = "filename" + f.attrs["history"] = ( + "Created with python " 'GMTGrid.save(%s,format="hdf")' % filename + ) + f.attrs["GMT_version"] = "NA" # Create the x array and the attributes of that GMT wants to see - xvar = np.linspace(self._geodict.xmin, - self._geodict.xmax, self._geodict.nx) + xvar = np.linspace(self._geodict.xmin, self._geodict.xmax, self._geodict.nx) x = f.create_dataset( - 'x', data=xvar, shape=xvar.shape, dtype=str(xvar.dtype)) - x.attrs['CLASS'] = 'DIMENSION_SCALE' - x.attrs['NAME'] = 'x' - x.attrs['_Netcdf4Dimid'] = 0 # no idea what this is - x.attrs['long_name'] = 'x' - x.attrs['actual_range'] = np.array((xvar[0], xvar[-1])) + "x", data=xvar, shape=xvar.shape, dtype=str(xvar.dtype) + ) + x.attrs["CLASS"] = "DIMENSION_SCALE" + x.attrs["NAME"] = "x" + x.attrs["_Netcdf4Dimid"] = 0 # no idea what this is + x.attrs["long_name"] = "x" + x.attrs["actual_range"] = np.array((xvar[0], xvar[-1])) # Create the x array and the attributes of that GMT wants to see - yvar = np.linspace(self._geodict.ymin, - self._geodict.ymax, self._geodict.ny) + yvar = np.linspace(self._geodict.ymin, self._geodict.ymax, self._geodict.ny) y = f.create_dataset( - 'y', data=yvar, shape=yvar.shape, dtype=str(yvar.dtype)) - y.attrs['CLASS'] = 'DIMENSION_SCALE' - y.attrs['NAME'] = 'y' - y.attrs['_Netcdf4Dimid'] = 1 # no idea what this is - y.attrs['long_name'] = 'y' - y.attrs['actual_range'] = np.array((yvar[0], yvar[-1])) + "y", data=yvar, shape=yvar.shape, dtype=str(yvar.dtype) + ) + y.attrs["CLASS"] = "DIMENSION_SCALE" + y.attrs["NAME"] = "y" + y.attrs["_Netcdf4Dimid"] = 1 # no idea what this is + y.attrs["long_name"] = "y" + y.attrs["actual_range"] = np.array((yvar[0], yvar[-1])) # create the z data set - z = f.create_dataset('z', data=np.flipud(self._data), - shape=self._data.shape, - dtype=str(self._data.dtype)) - z.attrs['long_name'] = 'z' + z = f.create_dataset( + "z", + data=np.flipud(self._data), + shape=self._data.shape, + dtype=str(self._data.dtype), + ) + z.attrs["long_name"] = "z" # zvar.attrs['_FillValue'] = array([ nan], dtype=float32) - z.attrs['actual_range'] = np.array( - (np.nanmin(self._data), np.nanmax(self._data))) + z.attrs["actual_range"] = np.array( + (np.nanmin(self._data), np.nanmax(self._data)) + ) # close the hdf file f.close() - elif format == 'native': - f = open(filename, 'wb') - f.write(struct.pack('I', self._geodict.nx)) - f.write(struct.pack('I', self._geodict.ny)) - f.write(struct.pack('I', 0)) # gridline registration - f.write(struct.pack('d', self._geodict.xmin)) - f.write(struct.pack('d', self._geodict.xmax)) - f.write(struct.pack('d', self._geodict.ymin)) - f.write(struct.pack('d', self._geodict.ymax)) - f.write(struct.pack('d', self._data.min())) - f.write(struct.pack('d', self._data.max())) - f.write(struct.pack('d', self._geodict.dx)) - f.write(struct.pack('d', self._geodict.dy)) - f.write(struct.pack('d', 1.0)) # scale factor to multiply data by - f.write(struct.pack('d', 0.0)) # offfset to add to data - f.write(struct.pack('80s', b'X units (probably degrees)')) - f.write(struct.pack('80s', b'Y units (probably degrees)')) - f.write(struct.pack('80s', b'Z units unknown')) - f.write(struct.pack('80s', b'')) # title - f.write(struct.pack('320s', b'Created with GMTGrid() class, a ' - b'product of the NEIC.')) # command - f.write(struct.pack('160s', b'')) # remark - if self._data.dtype not in [np.int16, np.int32, np.float32, - np.float64]: - msg = ('Data type of "%s" is not supported by the GMT native ' - 'format.') + elif format == "native": + f = open(filename, "wb") + f.write(struct.pack("I", self._geodict.nx)) + f.write(struct.pack("I", self._geodict.ny)) + f.write(struct.pack("I", 0)) # gridline registration + f.write(struct.pack("d", self._geodict.xmin)) + f.write(struct.pack("d", self._geodict.xmax)) + f.write(struct.pack("d", self._geodict.ymin)) + f.write(struct.pack("d", self._geodict.ymax)) + f.write(struct.pack("d", self._data.min())) + f.write(struct.pack("d", self._data.max())) + f.write(struct.pack("d", self._geodict.dx)) + f.write(struct.pack("d", self._geodict.dy)) + f.write(struct.pack("d", 1.0)) # scale factor to multiply data by + f.write(struct.pack("d", 0.0)) # offfset to add to data + f.write(struct.pack("80s", b"X units (probably degrees)")) + f.write(struct.pack("80s", b"Y units (probably degrees)")) + f.write(struct.pack("80s", b"Z units unknown")) + f.write(struct.pack("80s", b"")) # title + f.write( + struct.pack( + "320s", b"Created with GMTGrid() class, a " b"product of the NEIC." + ) + ) # command + f.write(struct.pack("160s", b"")) # remark + if self._data.dtype not in [np.int16, np.int32, np.float32, np.float64]: + msg = 'Data type of "%s" is not supported by the GMT native ' "format." raise DataSetException(msg % str(self._data.dtype)) # fpos1 = f.tell() # the left-right flip is necessary because of the way tofile() @@ -975,32 +1047,42 @@ def save(self, filename, format='netcdf'): f.close() @classmethod - def load(cls, filename, samplegeodict=None, resample=False, - method='linear', doPadding=False, padValue=np.nan): + def load( + cls, + filename, + samplegeodict=None, + resample=False, + method="linear", + doPadding=False, + padValue=np.nan, + ): """Create a GMTGrid object from a (possibly subsetted, resampled, or padded) GMT grid file. - :param gmtfilename: - Name of input file. - :param samplegeodict: - GeoDict used to specify subset bounds and resolution (if resample - is selected) - :param resample: - Boolean used to indicate whether grid should be resampled from the - file based on samplegeodict. - :param method: - If resample=True, resampling method to use ('nearest', 'linear', - 'cubic', 'quintic') - :param doPadding: - Boolean used to indicate whether, if samplegeodict is outside bounds - of grid, to pad values around the edges. - :param padValue: - Value to fill in around the edges if doPadding=True. - :returns: - GMTgrid instance (possibly subsetted, padded, or resampled) - :raises DataSetException: - * When sample bounds are outside (or too close to outside) the - bounds of the grid and doPadding=False. - * When the input file type is not recognized. + + Args: + filename (str): + Name of input file. + samplegeodict(:obj:`GeoDict`, optional): + GeoDict used to specify subset bounds and resolution (if resample + is selected) + resample (:obj:`bool`, optional): + Boolean used to indicate whether grid should be resampled from the + file based on samplegeodict. + method (:obj:`str`, optional): + If resample=True, resampling method to use ('nearest', 'linear', + 'cubic', 'quintic') + doPadding (:obj:`bool`, optional): + Boolean used to indicate whether, if samplegeodict is outside bounds + of grid, to pad values around the edges. + padValue (:obj:`float`, optional): + Value to fill in around the edges if doPadding=True. + Returns: + GMTgrid instance (possibly subsetted, padded, or resampled) + Raises: + DataSetException: + - When sample bounds are outside (or too close to outside) the + bounds of the grid and doPadding=False. + - When the input file type is not recognized. """ # get the geodict describing the source file, plus a boolean telling # us if the last column is a duplicate of the first column @@ -1021,9 +1103,9 @@ def load(cls, filename, samplegeodict=None, resample=False, # interpolation. if samplegeodict is not None: # parent static method - sampledict = cls.bufferBounds(samplegeodict, filegeodict, - resample=resample, - doPadding=doPadding) + sampledict = cls.bufferBounds( + samplegeodict, filegeodict, resample=resample, doPadding=doPadding + ) else: sampledict = filegeodict @@ -1037,8 +1119,10 @@ def load(cls, filename, samplegeodict=None, resample=False, except DataSetException: if doPadding: if not filegeodict.contains(sampledict): - data = np.ones((sampledict.ny, sampledict.nx), - dtype=np.float32)*padValue + data = ( + np.ones((sampledict.ny, sampledict.nx), dtype=np.float32) + * padValue + ) return cls(data=data, geodict=sampledict) sampledict = filegeodict.getIntersection(sampledict) @@ -1046,14 +1130,13 @@ def load(cls, filename, samplegeodict=None, resample=False, # sampledict.ymin, sampledict.ymax) data_range = cls.getDataRange( - filegeodict, sampledict, - first_column_duplicated=first_column_duplicated) + filegeodict, sampledict, first_column_duplicated=first_column_duplicated + ) data, geodict = cls.readFile(filename, data_range) - if (doPadding or resample) and 'float' not in str(data.dtype): + if (doPadding or resample) and "float" not in str(data.dtype): data = data.astype(np.float32) # parent static method - pad_dict = cls.getPadding( - filegeodict, samplegeodict, doPadding=doPadding) + pad_dict = cls.getPadding(filegeodict, samplegeodict, doPadding=doPadding) data, geodict = cls.padGrid(data, geodict, pad_dict) if doPadding: data[np.isinf(data)] = padValue @@ -1079,9 +1162,9 @@ def __getitem__(self, *args): col = args[0][1] ny = self.ny nx = self.nx - if row < 0 or row > ny-1: + if row < 0 or row > ny - 1: raise Exception("Row index out of bounds") - if col < 0 or col > nx-1: + if col < 0 or col > nx - 1: raise Exception("Row index out of bounds") idx = nx * row + col # offset = 0 @@ -1107,30 +1190,30 @@ def __getitem__(self, *args): colstep = 1 # error checking - if rowstart < 0 or rowstart > ny-1: + if rowstart < 0 or rowstart > ny - 1: raise Exception("Row index out of bounds") if rowend < 0 or rowend > ny: raise Exception("Row index out of bounds") - if colstart < 0 or colstart > nx-1: + if colstart < 0 or colstart > nx - 1: raise Exception("Col index out of bounds") if colend < 0 or colend > nx: raise Exception("Col index out of bounds") - colcount = (colend-colstart) - rowcount = (rowend-rowstart) - outrows = np.ceil(rowcount/rowstep) - outcols = np.ceil(colcount/colstep) + colcount = colend - colstart + rowcount = rowend - rowstart + outrows = np.ceil(rowcount / rowstep) + outcols = np.ceil(colcount / colstep) data = np.zeros([outrows, outcols], dtype=self.dtype) outrow = 0 for row in range(int(rowstart), int(rowend), int(rowstep)): # just go to the beginning of the row, we're going to read in # the whole line - idx = nx*row + idx = nx * row # offset = self.dwidth*idx # beginning of row - line = self.array[idx:idx+nx] + line = self.array[idx : idx + nx] data[outrow, :] = line[colstart:colend:colstep] - outrow = outrow+1 + outrow = outrow + 1 else: raise Exception("Unsupported __getitem__ input %s" % (str(args))) - return(data) + return data diff --git a/mapio/grid2d.py b/mapio/grid2d.py index efc60eb..989fcdb 100755 --- a/mapio/grid2d.py +++ b/mapio/grid2d.py @@ -41,15 +41,17 @@ def __init__(self, data=None, geodict=None): """ Construct a Grid object. - :param data: - A 2D numpy array (can be None). - :param geodict: - A GeoDict Object (or None) containing the following fields: - :returns: + Args: + data (:obj:`array`, optional): + A 2D numpy array (can be None). + geodict (:obj:`GeoDict`, optional): + A GeoDict Object (or None) containing the following fields: + Returns: A Grid2D object. - :raises DataSetException: - When data is not 2D or the number of rows and columns do not - match the geodict. + Raises: + DataSetException: + When data is not 2D or the number of rows and columns do not + match the geodict. """ if data is not None and geodict is not None: # complain if data is not 2D (squeeze 1d dimensions out) @@ -75,12 +77,13 @@ def checkFirstColumnDuplicated(geodict): """Check to see if the first column in a file described by geodict is duplicated by the last column. - :param geodict: - GeoDict object which may have duplicate column. - :returns: - Tuple containing: - - GeoDict object representing a grid with the last column removed. - - Boolean indicating whether the last column was a duplicate. + Args: + geodict (GeoDict): + GeoDict object which may have duplicate column. + Returns: + Tuple containing: + - GeoDict object representing a grid with the last column removed. + - Boolean indicating whether the last column was a duplicate. """ first_column_duplicated = False cols_per_degree = 1 / geodict.dx @@ -114,24 +117,25 @@ def getDataRange(fgeodict, sampledict, first_column_duplicated=False, padding=No """For a given set of input bounds and information about a file, determine the rows and columns for bounds. - :param fgeodict: - GeoDict object for a given file. - :param sampledict: - Sampling GeoDict object. - :param first_column_duplicated: - Boolean indicating whether the last column in a file is a - duplicate of the first column. - :returns: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. + Args: + fgeodict (GeoDict): + GeoDict object for a given file. + sampledict (GeoDict): + Sampling GeoDict object. + first_column_duplicated (bool): + indicating whether the last column in a file is a + duplicate of the first column. + Returns: + Dictionary containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. """ data_range = {} @@ -206,15 +210,16 @@ def verifyBounds(filegeodict, samplegeodict, resample=False): """Ensure that the two grids represented at least 1) intersect and 2) are aligned if resampling is True. - :param filegeodict: - GeoDict object describing file. - :param samplegeodict: - GeoDict object describing grid to use for sampling. - :param resample: - Boolean indicating whether we want to resample. - :raises: - DataSetException when geodicts do not intersect or if the grids - are not aligned. + Args: + filegeodict (GeoDict): + object describing file. + samplegeodict (GeoDict): + object describing grid to use for sampling. + resample (bool): + indicating whether we want to resample. + Raises: + DataSetException when geodicts do not intersect or if the grids + are not aligned. """ if samplegeodict is not None and not filegeodict.intersects(samplegeodict): msg = "Input samplegeodict must at least intersect with the " "file bounds" @@ -238,15 +243,16 @@ def bufferBounds( Buffer pixels shoud be at filegeodict resolution. - :param samplegeodict: - GeoDict object describing grid to use for sampling. - :param filegeodict: - GeoDict object describing file. - :param resample: - Boolean indicating whether we want to resample. - :param buffer_pixels: - Number of pixels to buffer bounds in any possible direction. - :returns: + Args: + samplegeodict (GeoDict): + object describing grid to use for sampling. + filegeodict (GeoDict): + object describing file. + resample (bool): + Boolean indicating whether we want to resample. + buffer_pixels (int): + Number of pixels to buffer bounds in any possible direction. + Returns: GeoDict which has been buffered by the appropriate number of pixels. """ if not resample: @@ -318,19 +324,20 @@ def bufferBounds( def padGrid(data, geodict, pad_dict): """Pad input data array with pixels specified by pad_dict on each side. - :param data: - 2D numpy array of data. - :param geodict: - GeoDict object describing data. - :param pad_dict: - A dictionary containing fields: - - padleft The number of padding pixels on the left edge. - - padright The number of padding pixels on the right edge. - - padbottom The number of padding pixels on the bottom edge. - - padtop The number of padding pixels on the top edge. - :returns: - Tuple of (data,geodict) where data has been padded and geodict - represents new padded data. + Args: + data (array): + 2D numpy array of data. + geodict (GeoDict): + GeoDict object describing data. + pad_dict (dict): + containing fields: + - padleft The number of padding pixels on the left edge. + - padright The number of padding pixels on the right edge. + - padbottom The number of padding pixels on the bottom edge. + - padtop The number of padding pixels on the top edge. + Returns: + Tuple of (data,geodict) where data has been padded and geodict + represents new padded data. """ if ( pad_dict["padleft"] == 0 @@ -389,22 +396,24 @@ def getPadding(filegeodict, samplegeodict, doPadding=False): """Determine how many pixels of padding there need to be on each side of requested grid. - :param filegeodict: - GeoDict object specifying the spatial information from a source file. - :param samplegeodict: - GeoDict object specifying the spatial information for a desired - sampling regime. - :param doPadding: - Boolean indicating that user does or does not want to add any padding. - :raises DataSetException: - When resampling is False and filegeodict and samplegeodict are - not pixel aligned. - :returns: - A dictionary containing fields: - - padleft The number of padding pixels on the left edge. - - padright The number of padding pixels on the right edge. - - padbottom The number of padding pixels on the bottom edge. - - padtop The number of padding pixels on the top edge. + Args: + filegeodict (GeoDict): + specifying the spatial information from a source file. + samplegeodict (GeoDict): + specifying the spatial information for a desired + sampling regime. + doPadding (bool): + Boolean indicating that user does or does not want to add any padding. + Raises: + DataSetException: + When resampling is False and filegeodict and samplegeodict are + not pixel aligned. + Returns: + A dictionary containing fields: + - padleft The number of padding pixels on the left edge. + - padright The number of padding pixels on the right edge. + - padbottom The number of padding pixels on the bottom edge. + - padtop The number of padding pixels on the top edge. """ if not doPadding: pad_dict = {"padleft": 0, "padright": 0, "padbottom": 0, "padtop": 0} @@ -471,8 +480,9 @@ def getPadding(filegeodict, samplegeodict, doPadding=False): def __repr__(self): """ String representation of a Grid2D object. - :returns: - String containing description of Grid2D object. + + Returns: + String containing description of Grid2D object. """ xmin, xmax, ymin, ymax = ( self._geodict.xmin, @@ -516,14 +526,15 @@ def _createSampleData(self, M, N): """Used for internal testing - create an NxN grid with lower left corner at 0.5,0.5, dx/dy = 1.0. - :param M: - Number of rows in output grid - :param N: - Number of columns in output grid - :returns: - GMTGrid object where data values are an NxN array of values - from 0 to N-squared minus 1, and geodict - lower left corner is at 0.5/0.5 and cell dimensions are 1.0. + Args: + M (int): + Number of rows in output grid + N (int): + Number of columns in output grid + Returns: + GMTGrid object where data values are an NxN array of values + from 0 to N-squared minus 1, and geodict + lower left corner is at 0.5/0.5 and cell dimensions are 1.0. """ data = np.arange(0, M * N).reshape(M, N) # arange gives int64 by default, not supported by netcdf3 @@ -556,19 +567,20 @@ def readFile(filename, data_range): """Read in data from the given file, at the pixels specified in data_range. - :param filename: - Name of file to read. - :param data_range: - Dictionary containing fields: - - iulx1 Upper left X of first (perhaps only) segment. - - iuly1 Upper left Y of first (perhaps only) segment. - - ilrx1 Lower right X of first (perhaps only) segment. - - ilry1 Lower right Y of first (perhaps only) segment. - (if bounds cross 180 meridian...) - - iulx2 Upper left X of second segment. - - iuly2 Upper left Y of second segment. - - ilrx2 Lower right X of second segment. - - ilry2 Lower right Y of second segment. + Args: + filename (str): + Name of file to read. + data_range (dict): + Dictionary containing fields: + - iulx1 Upper left X of first (perhaps only) segment. + - iuly1 Upper left Y of first (perhaps only) segment. + - ilrx1 Lower right X of first (perhaps only) segment. + - ilry1 Lower right Y of first (perhaps only) segment. + (if bounds cross 180 meridian...) + - iulx2 Upper left X of second segment. + - iuly2 Upper left Y of second segment. + - ilrx2 Lower right X of second segment. + - ilry2 Lower right Y of second segment. """ raise NotImplementedError( "readFile is only implemented in Grid2D " "subclasses." @@ -580,10 +592,11 @@ def copyFromGrid(cls, grid): Copy constructor - can be used to create an instance of any Grid2D subclass from another. - :param grid: - Any Grid2D instance. - :returns: - A copy of the data in that input Grid2D instance. + Args: + grid (Grid2D): + Any Grid2D instance. + Returns: + A copy of the data in that input Grid2D instance. """ if not isinstance(grid, Grid2D): raise DataSetException( @@ -604,21 +617,22 @@ def _createSections(self, bounds, geodict, firstColumnDuplicated, isScanLine=Fal """Given a grid that goes from -180 to 180 degrees, figure out the two pixel regions that up both sides of the subset. - :param bounds: - Tuple of (xmin,xmax,ymin,ymax) - :param geodict: - Geodict dictionary - :param firstColumnDuplicated: - Boolean indicating whether this is a global data set where the - first and last columns are identical. - :param isScanLine: - Boolean indicating whether this array is in scan line order - (pixel[0,0] is the geographic upper left). - :returns: - Two tuples of 4 elements each - (iulx,iuly,ilrx,ilry). The - first tuple defines the pixel offsets for the left - side of the subsetted region, and the second tuple - defines the pixel offsets for the right side. + Args: + bounds (tuple): + Tuple of (xmin,xmax,ymin,ymax) + geodict (Geodict): + Geodict dictionary + firstColumnDuplicated (bool): + Boolean indicating whether this is a global data set where the + first and last columns are identical. + isScanLine (bool): + Boolean indicating whether this array is in scan line order + (pixel[0,0] is the geographic upper left). + Returns: + Two tuples of 4 elements each - (iulx,iuly,ilrx,ilry). The + first tuple defines the pixel offsets for the left + side of the subsetted region, and the second tuple + defines the pixel offsets for the right side. """ (bxmin, bxmax, bymin, bymax) = bounds ulx = geodict.xmin @@ -653,20 +667,22 @@ def _createSections(self, bounds, geodict, firstColumnDuplicated, isScanLine=Fal def getData(self): """Return a reference to the data inside the Grid. - :returns: - A reference to a 2D numpy array. + Returns: + A reference to a 2D numpy array. """ return self._data # should we return a copy of the data? def setData(self, data): """Set the data inside the Grid. - :param data: - A 2D numpy array. - :raises: - DataSetException if the number of rows and columns do not match - the the internal GeoDict, or if the input - is not a numpy array. + + Args: + data (array): + A 2D numpy array. + Raises: + DataSetException if the number of rows and columns do not match + the the internal GeoDict, or if the input + is not a numpy array. """ if not isinstance(data, np.ndarray): raise DataSetException("setData() input is not a numpy array.") @@ -684,7 +700,8 @@ def setData(self, data): def getGeoDict(self): """ Return a reference to the geodict inside the Grid. - :returns: + + Returns:: A reference to a dictionary (see constructor). """ return self._geodict.copy() @@ -693,7 +710,7 @@ def getBounds(self): """ Return the lon/lat range of the data. - :returns: + Returns: Tuple of (lonmin,lonmax,latmin,latmax) """ return ( @@ -721,8 +738,9 @@ def applyNaN(self, force=False): losing precision. Otherwise, this method will raise an OverflowError unless the force option is set to True. - :param force: - Boolean indicating whether to override OverflowError (see Usage). + Args: + force (bool): + Boolean indicating whether to override OverflowError (see Usage). """ nodata = self._geodict.nodata if nodata is None or np.isnan(nodata) or np.isnan(self._data).any(): @@ -760,21 +778,23 @@ def applyNaN(self, force=False): def subdivide(self, finerdict, cellFill="max"): """Subdivide the cells of the host grid into finer-resolution cells. - :param finerdict: - GeoDict object defining a grid with a finer resolution than the - host grid. - :param cellFill: - String defining how to fill cells that span more than one host - grid cell. - Choices are: - 'max': Choose maximum value of host grid cells. - 'min': Choose minimum value of host grid cells. - 'mean': Choose mean value of host grid cells. - :returns: - Grid2D instance with host grid values subdivided onto finer grid. - :raises DataSetException: - When finerdict is not a) finer resolution or b) does not - intersect.x or cellFill is not valid. + Args: + finerdict (GeoDict): + GeoDict object defining a grid with a finer resolution than the + host grid. + cellFill (:obj:`str`, optional): + String defining how to fill cells that span more than one host + grid cell. + Choices are: + 'max': Choose maximum value of host grid cells. + 'min': Choose minimum value of host grid cells. + 'mean': Choose mean value of host grid cells. + Returns: + Grid2D instance with host grid values subdivided onto finer grid. + Raises: + DataSetException: + When finerdict is not a) finer resolution or b) does not + intersect.x or cellFill is not valid. """ fillvals = ["min", "max", "mean"] if cellFill not in fillvals: @@ -921,16 +941,22 @@ def subdivide(self, finerdict, cellFill="max"): def cut(self, xmin, xmax, ymin, ymax, align=False): """Cut out a section of Grid and return it. - :param xmin: Longitude coordinate of upper left pixel, must be - aligned with Grid. - :param xmax: Longitude coordinate of lower right pixel, must be - aligned with Grid. - :param ymin: Latitude coordinate of upper left pixel, must be - aligned with Grid. - :param ymax: Latitude coordinate of lower right pixel, must be - aligned with Grid. - :param align: Boolean indicating whether input boundaries - should be modified to align with host grid. + Args: + xmin (float): + Longitude coordinate of upper left pixel, must be + aligned with Grid. + xmax (float): + Longitude coordinate of lower right pixel, must be + aligned with Grid. + ymin (float): + Latitude coordinate of upper left pixel, must be + aligned with Grid. + ymax (float:): + Latitude coordinate of lower right pixel, must be + aligned with Grid. + align (bool): + Boolean indicating whether input boundaries + should be modified to align with host grid. """ td1 = GeoDict.createDictFromBox( xmin, xmax, ymin, ymax, self._geodict.dx, self._geodict.dy, inside=True @@ -958,18 +984,20 @@ def getValue(self, lat, lon, method="nearest", default=None): """Return numpy array at given latitude and longitude (using nearest neighbor). - :param lat: - Latitude (in decimal degrees) of desired data value. - :param lon: - Longitude (in decimal degrees) of desired data value. - :param method: - Interpolation method, one of ('nearest','linear','cubic','quintic') - :param default: - Default value to return when lat/lon is outside of grid bounds. - :return: + Args: + lat (float): + Latitude (in decimal degrees) of desired data value. + lon (float): + Longitude (in decimal degrees) of desired data value. + method (:obj:`str`, optional): + Interpolation method, one of ('nearest','linear','cubic','quintic') + default: + Default value to return when lat/lon is outside of grid bounds. + Returns: Value at input latitude,longitude position. - :raises DataSetException: - When lat/lon is outside of bounds and default is None. + Raises: + DataSetException: + When lat/lon is outside of bounds and default is None. """ if method == "nearest": @@ -1010,12 +1038,13 @@ 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. + Args: + row (int): + Row dimension index into internal data array. + col (int): + Column dimension index into internal data array. + Returns: + Tuple of latitude and longitude. """ return self._geodict.getLatLon(row, col) @@ -1023,15 +1052,16 @@ 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: - Input longitude. - :param returnFloat: - Boolean indicating whether floating point row/col coordinates - should be returned. - :returns: - Tuple of row and column. + Args: + lat (float): + Input latitude. + lon (float): + Input longitude. + returnFloat (bool): + Boolean indicating whether floating point row/col coordinates + should be returned. + Returns: + Tuple of row and column. """ return self._geodict.getRowCol(lat, lon, returnFloat) @@ -1056,7 +1086,6 @@ def _getInterpCoords(self, geodict): # if not, adjust one to the other hostxmin = self._geodict.xmin hostxmax = self._geodict.xmax - hostymin = self._geodict.ymin hostymax = self._geodict.ymax # there are three cases for longitude coordinates we need to @@ -1123,22 +1152,24 @@ def interpolate2(self, geodict, method="linear"): this method is 5 to 7 times faster than interpolateToGrid. - :param geodict: - geodict dictionary from another grid whose extents are - inside the extent of this grid. - :param method: - Optional interpolation method - ['linear', 'cubic','nearest'] - :raises DataSetException: - If the Grid object upon which this function is being called is - not completely contained by the grid to which this Grid is being - resampled. - :raises DataSetException: - If the method is not one of ['nearest','linear','cubic'] - If the resulting interpolated grid shape does not match input - geodict. - :returns: - A new instance of the Grid2D class or subclass with interpolated - data. + Args: + geodict (GeoDict): + geodict dictionary from another grid whose extents are + inside the extent of this grid. + method (:obj:`str`, optional): + Optional interpolation method - ['linear', 'cubic','nearest'] + Raises: + DataSetException: + If the Grid object upon which this function is being called is + not completely contained by the grid to which this Grid is being + resampled. + DataSetException: + If the method is not one of ['nearest','linear','cubic'] + If the resulting interpolated grid shape does not match input + geodict. + Returns: + A new instance of the Grid2D class or subclass with interpolated + data. """ if not self._geodict.contains(geodict): msg = "Input geodict not fully contained by host geodict." @@ -1196,22 +1227,24 @@ def interpolateToGrid(self, geodict, method="linear"): Given a geodict specifying another grid extent and resolution, resample current grid to match. - :param geodict: - geodict dictionary from another grid whose extents are inside - the extent of this grid. - :param method: - Optional interpolation method - ['linear', 'cubic','nearest'] - :raises DataSetException: - If the Grid object upon which this function is being called is - not completely contained by the grid to which this Grid is being - resampled. - :raises DataSetException: - If the method is not one of ['nearest','linear','cubic'] - If the resulting interpolated grid shape does not match input - geodict. - :returns: - A new instance of the Grid2D class or subclass with interpolated - data. + Args: + geodict (GeoDict): + geodict dictionary from another grid whose extents are inside + the extent of this grid. + method (:obj:`str`, optional):: + Optional interpolation method - ['linear', 'cubic','nearest'] + Raises: + DataSetException: + If the Grid object upon which this function is being called is + not completely contained by the grid to which this Grid is being + resampled. + DataSetException: + If the method is not one of ['nearest','linear','cubic'] + If the resulting interpolated grid shape does not match input + geodict. + Returns: + A new instance of the Grid2D class or subclass with interpolated + data. """ if method not in ["linear", "cubic", "nearest"]: raise DataSetException( @@ -1317,39 +1350,39 @@ def rasterizeFromGeometry( of a shape (point, line, polygon) inside a cell turns that cell "on". - :param shapes: - One of: - - One shapely geometry object (Point, Polygon, etc.) or a - sequence of such objects - - One GeoJSON like object or sequence of such objects. - (http://geojson.org/) - - A tuple of (geometry,value) or sequence of (geometry,value). - :param geodict: - GeoDict object which defines the grid onto which the shape values - should be "burned". - :param burnValue: - Optional value which will be used to set the value of the pixels - if there is no - value in the geometry field. - :param fillValue: - Optional value which will be used to fill the cells not touched - by any geometry. - :param mustContainCenter: - Optional boolean which indicates whether the geometry must touch - the center of the cell or merely be inside the cell in order to - set the value. - :raises DataSetException: - When geometry input is not a subclass of - shapely.geometry.base.BaseGeometry. - :returns: - Grid2D object. - This method is a thin wrapper around rasterio->features->rasterize(), - documented here: - https://github.com/mapbox/rasterio/blob/master/docs/features.rst - - which is itself a Python wrapper around the functionality found - in gdal_rasterize, documented here: - http://www.gdal.org/gdal_rasterize.html + Args: + shapes (geometry obj, GeoJson obj or tuple): + One of: + - One shapely geometry object (Point, Polygon, etc.) or a + sequence of such objects + - One GeoJSON like object or sequence of such objects. + (http://geojson.org/) + - A tuple of (geometry,value) or sequence of (geometry,value). + geodict (GeoDict): + GeoDict object which defines the grid onto which the shape values + should be "burned". + burnValue (:obj:`float`, optional): + Optional value which will be used to set the value of the pixels + if there is no value in the geometry field. + fillValue (:obj:`float`, optional): + Optional value which will be used to fill the cells not touched + by any geometry. + mustContainCenter (:obj:`bool`, optional):: + Optional boolean which indicates whether the geometry must touch + the center of the cell or merely be inside the cell in order to + set the value. + Raises: + DataSetException: + When geometry input is not a subclass of + shapely.geometry.base.BaseGeometry. + Returns: + Grid2D object. + This method is a thin wrapper around rasterio->features->rasterize(), + documented here: + https://github.com/mapbox/rasterio/blob/master/docs/features.rst + which is itself a Python wrapper around the functionality found + in gdal_rasterize, documented here: + http://www.gdal.org/gdal_rasterize.html """ # check the type of shapes # features.rasterize() documentation says this: @@ -1359,7 +1392,8 @@ def rasterizeFromGeometry( # create list of allowable types if sys.version_info.major == 2: - types = (int, float, long) + # long variable was previously here but was removed due to being undefined + types = (int, float) else: types = (int, float) @@ -1457,16 +1491,19 @@ def rasterizeFromGeometry( def project(self, projection, method="bilinear"): """Project Grid2D data into desired projection. - :param projection: - Valid proj4 projection string. - :param method: - One of the sampling methods described here: - https://mapbox.github.io/rasterio/topics/resampling.html#resampling-methods - :raises DataSetException: - If input projection is not a valid Proj4 string. - If method is not a valid resampling method found in above URL. - :returns: - Re-projected Grid2D object. + Args: + projection: + Valid proj4 projection string. + method: + One of the sampling methods described here: + mapbox.github.io + https://mapbox.github.io/rasterio/topics/resampling.html#resampling-methods + Raises: + DataSetException: + If input projection is not a valid Proj4 string. + If method is not a valid resampling method found in above URL. + Returns: + Re-projected Grid2D object. """ # hack to handle projections that wrap around the 180 meridian. crosses = self._geodict.xmax < self._geodict.xmin diff --git a/mapio/gridbase.py b/mapio/gridbase.py index b862c31..ef9784c 100755 --- a/mapio/gridbase.py +++ b/mapio/gridbase.py @@ -14,27 +14,38 @@ class Grid(DataSet): 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 def getFileGeoDict(filename): """ - Abstract method to return the bounding box, resolution, and shape of a file in whatever Grid format. - :param filename: - The path to the filename of whatever grid format this is being implemented in. - :returns: - A geodict specifying the bounding box, resolution, and shape of the data in a file. + Abstract method to return the bounding box, resolution, and shape of a file in + whatever Grid format. + + Args: + filename: + The path to the filename of whatever grid format this is being + implemented in. + Returns: + A geodict specifying the bounding box, resolution, and shape of the data in + a file. """ raise NotImplementedError @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: - The name of the file whose resolution/extent should be used. - :param geodict: - The geodict which is used as the base for finding the bounds for this file guaranteed to be inside of this geodict. - :raises NotImplementedError: - Always in base class + Abstract method to return a geodict for this file that is guaranteed to be + inside the input geodict defined, without resampling. + + Args: + filename: + The name of the file whose resolution/extent should be used. + geodict: + The geodict which is used as the base for finding the bounds for this + file guaranteed to be inside of this geodict. + Raises: + NotImplementedError: + Always in base class """ raise NotImplementedError @@ -42,9 +53,17 @@ def getBoundsWithin(filename, geodict): 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) + paddict.xmin, + paddict.xmax, + paddict.ymin, + paddict.ymax, + ) gxmin, gxmax, gymin, gymax = ( - geodict.xmin, geodict.xmax, geodict.ymin, geodict.ymax) + geodict.xmin, + geodict.xmax, + geodict.ymin, + geodict.ymax, + ) dx, dy = (geodict.dx, geodict.dy) ny, nx = (geodict.ny, geodict.nx) @@ -71,23 +90,22 @@ def _getPadding(cls, geodict, paddict, padvalue): # now figure out what the new bounds are outdict = {} - outdict['nx'] = int(nx) - outdict['ny'] = int(ny + bottompad.shape[0] + toppad.shape[0]) + 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['dx'] = dx - outdict['dy'] = 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']) + reqfields = set(["xmin", "xmax", "ymin", "ymax", "dx", "dy", "ny", "nx"]) if not reqfields.issubset(set(geodict.keys())): return False return True @@ -95,9 +113,12 @@ def checkGeoDict(cls, geodict): @abc.abstractmethod 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. + Abstract method to calculate average values for cells of larger size than the + current grid. + + Args: + geodict: + Geodict that defines the new coarser grid. """ raise NotImplementedError @@ -105,11 +126,14 @@ def blockmean(self, geodict): def loadFromCloud(cls, cloud, geodict): """ Create a grid from a Cloud instance (scattered XY data). - :param cloud: - A Cloud instance containing scattered XY data. - :param geodict: - A geodict object where ny/nx are optional (will be calculated from bounds/cell dimensions) - :returns: + + Args: + cloud: + A Cloud instance containing scattered XY data. + geodict: + A geodict object where ny/nx are optional (will be calculated from + bounds/cell dimensions) + Returns: An instance of a Grid object. """ raise NotImplementedError @@ -117,9 +141,7 @@ def loadFromCloud(cls, cloud, geodict): @staticmethod def getLatLonMesh(geodict): if geodict.xmax < geodict.xmin: - lons = np.linspace(geodict.xmin, - geodict.xmax + 360, - num=geodict.nx) + 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) @@ -131,38 +153,40 @@ def getGeoDict(self): """ Return a reference to the geodict inside the Grid - :returns: - A reference to a dictionary (see constructor). + 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): - """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. + """Return geographic coordinates (lat/lon decimal degrees) for given data row + and column. + + Args: + row: + Row dimension index into internal data array. + 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): - """Return data row and column from given geographic coordinates (lat/lon decimal degrees). - - :param lat: - Input latitude. - :param lon: - Input longitude. - :param returnFloat: - Boolean indicating whether floating point row/col coordinates should be returned. - :returns: - Tuple of row and column. + """Return data row and column from given geographic coordinates + (lat/lon decimal degrees). + + Args: + lat: + Input latitude. + lon: + Input longitude. + returnFloat: + Boolean indicating whether floating point row/col coordinates + should be returned. + Returns: + Tuple of row and column. """ - raise NotImplementedError( - 'getRowCol method not implemented in base class') + raise NotImplementedError("getRowCol method not implemented in base class") diff --git a/mapio/gridcontainer.py b/mapio/gridcontainer.py index 35a3c6d..9cd4e20 100644 --- a/mapio/gridcontainer.py +++ b/mapio/gridcontainer.py @@ -1,22 +1,12 @@ #!/usr/bin/env python - -# stdlib imports -from datetime import datetime -import collections -import time -import io -import copy - # third party imports -import h5py -import numpy as np from impactutils.io.container import HDFContainer # local imports from mapio.grid2d import Grid2D from mapio.geodict import GeoDict -GROUPS = {'grid': 'grids'} +GROUPS = {"grid": "grids"} class GridHDFContainer(HDFContainer): @@ -24,32 +14,34 @@ def setGrid(self, name, grid, metadata=None, compression=True): """Store a Grid2D object as a dataset. Args: - name (str): Name of the Grid2D object to be stored. - grid (Grid2D): Grid2D object to be stored. - metadata (dict): Simple dictionary (values of strings and numbers). - compression (bool): Boolean indicating whether dataset should be compressed - using the gzip algorithm. - + name (str): + Name of the Grid2D object to be stored. + grid (Grid2D): + Grid2D object to be stored. + metadata (dict): + Simple dictionary (values of strings and numbers). + compression (bool): + Boolean indicating whether dataset should be compressed + sing the gzip algorithm. Returns: HDF Dataset containing grid and metadata. """ if compression: - compression = 'gzip' + compression = "gzip" else: compression = None - grid_name = '%s' % name - if GROUPS['grid'] not in self._hdfobj: - grid_group = self._hdfobj.create_group(GROUPS['grid']) + grid_name = "%s" % name + if GROUPS["grid"] not in self._hdfobj: + grid_group = self._hdfobj.create_group(GROUPS["grid"]) else: - grid_group = self._hdfobj[GROUPS['grid']] + grid_group = self._hdfobj[GROUPS["grid"]] array_metadata = grid.getGeoDict().asDict() data = grid.getData() if metadata is not None: array_metadata.update(metadata) - dset = grid_group.create_dataset( - grid_name, data=data, compression=compression) + dset = grid_group.create_dataset(grid_name, data=data, compression=compression) for key, value in array_metadata.items(): dset.attrs[key] = value return dset @@ -61,15 +53,13 @@ def getGrid(self, name): Args: name (str): The name of the Grid2D object stored in the container. - Returns: (tuple) Grid2D object, and a dictionary of metadata. """ - grid_name = '%s' % name - grid_group = self._hdfobj[GROUPS['grid']] + grid_name = "%s" % name + grid_group = self._hdfobj[GROUPS["grid"]] if grid_name not in grid_group: - raise LookupError('Grid %s not in %s' - % (name, self.getFileName())) + raise LookupError("Grid %s not in %s" % (name, self.getFileName())) dset = grid_group[grid_name] data = dset[()] @@ -85,9 +75,9 @@ def getGrids(self): Returns: (list) List of names of Grid2D objects stored in container. """ - if GROUPS['grid'] not in self._hdfobj: + if GROUPS["grid"] not in self._hdfobj: return [] - grids = list(self._hdfobj[GROUPS['grid']].keys()) + grids = list(self._hdfobj[GROUPS["grid"]].keys()) return grids def dropGrid(self, name): @@ -99,18 +89,17 @@ def dropGrid(self, name): The name of the Grid2D object to be deleted. """ - mgrid = '%s' % name - grid_group = self._hdfobj[GROUPS['grid']] + mgrid = "%s" % name + grid_group = self._hdfobj[GROUPS["grid"]] if mgrid not in grid_group: - raise LookupError('Grid %s not in %s' - % (name, self._hdfobj.filename)) + raise LookupError("Grid %s not in %s" % (name, self._hdfobj.filename)) del grid_group[mgrid] def _split_dset_attrs(dset): geodict = {} metadata = {} - grid_keys = ['xmin', 'xmax', 'ymin', 'ymax', 'nx', 'ny', 'dx', 'dy'] + grid_keys = ["xmin", "xmax", "ymin", "ymax", "nx", "ny", "dx", "dy"] for key, value in dset.attrs.items(): if key in grid_keys: geodict[key] = value diff --git a/mapio/mapcity.py b/mapio/mapcity.py index e390731..5c8fe45 100644 --- a/mapio/mapcity.py +++ b/mapio/mapcity.py @@ -1,65 +1,73 @@ -#stdlib imports -import os.path -import warnings - -#local imports +# local imports from .city import Cities -from .dataset import DataSetException,DataSetWarning +from .dataset import DataSetException -#third party imports +# third party imports import matplotlib.font_manager -import matplotlib.pyplot as plt -import numpy as np + class MapCities(Cities): """ A subclass of Cities that can remove cities whose labels on a map intersect with larger cities. - """ - def __init__(self,dataframe): + """ + + def __init__(self, dataframe): """Construct a MapCities object from a pandas DataFrame. - :param dataframe: - pandas DataFrame, where each row represents a city. - Columns include: - - name Name of the city (required). - - lat Latitude of city (required). - - lon Longitude of city (required). - - pop Population of city (optional). - - iscap Boolean indicating capital status (optional). - - placement String indicating where city label - should be placed relative to city coordinates, - one of: E,W,N,S,NE,SE,SW,NW (optional). - -xoff Longitude offset for label relative to city coordinates (optional). - -yoff Latitude offset for label relative to city coordinates (optional). - :raises DataSetException: - When any of required columns are missing. - :returns: - MapCities instance. + + Args: + dataframe (dataframe): + pandas DataFrame, where each row represents a city. + Columns include: + - name Name of the city (required). + - lat Latitude of city (required). + - lon Longitude of city (required). + - pop Population of city (optional). + - iscap Boolean indicating capital status (optional). + - placement String indicating where city label + should be placed relative to city coordinates, + one of: E,W,N,S,NE,SE,SW,NW (optional). + -xoff Longitude offset for label relative to city coordinates + (optional). + -yoff Latitude offset for label relative to city coordinates + (optional). + Raises: + DataSetException: + When any of required columns are missing. + Returns: + MapCities instance. """ - self.SUGGESTED_FONTS = ['Bitstream Vera Sans','Times New Roman', - 'Courier New','Palatino LinoType', - 'Arial','Tahoma'] - self.DEFAULT_FONT = 'Times New Roman' + self.SUGGESTED_FONTS = [ + "Bitstream Vera Sans", + "Times New Roman", + "Courier New", + "Palatino LinoType", + "Arial", + "Tahoma", + ] + self.DEFAULT_FONT = "Times New Roman" self.DEFAULT_FONT_SIZE = 10.0 if len(set(dataframe.columns).intersection(set(self.REQFIELDS))) < 3: - raise DataSetException('Missing some of required keys: %s' % self.REQFIELDS) + raise DataSetException("Missing some of required keys: %s" % self.REQFIELDS) self._dataframe = dataframe.copy() self._fontlist = [f.name for f in matplotlib.font_manager.fontManager.ttflist] self._fontlist.sort() - + def limitByMapCollision(self): """Create a smaller Cities dataset by removing smaller cities whose bounding boxes collide with larger cities. - :returns: - New Cities instance where smaller colliding cities have been removed. + + Returns: + New Cities instance where smaller colliding cities have been removed. """ raise NotImplementedError - + def getFontList(self): """Return list of supported font names on this system. - :returns: - List of supported font names on this system. + + Returns: + List of supported font names on this system. """ return self._fontlist diff --git a/mapio/multihaz.py b/mapio/multihaz.py index 80ff1fb..2e82e4b 100755 --- a/mapio/multihaz.py +++ b/mapio/multihaz.py @@ -1,18 +1,17 @@ #!/usr/bin/env python -#python 3 compatibility +# python 3 compatibility from __future__ import print_function -#stdlib imports +# stdlib imports import sys import collections import datetime import time import os.path -#third party imports +# third party imports import h5py -from scipy.io import netcdf import numpy as np from .multiple import MultiGrid from .shake import ShakeGrid @@ -21,16 +20,17 @@ from .geodict import GeoDict - class MultiHazardGrid(MultiGrid): - def __init__(self,layers,geodict,origin,header,metadata=None): + def __init__(self, layers, geodict, origin, header, metadata=None): """Construct a MultiHazardGrid object. - - :param layers: - OrderedDict containing earthquake-induced hazard data layers (keys are 'pga', etc., values are 2D arrays of data). - :param geodict: + + Args: + layers (OrderedDict): + OrderedDict containing earthquake-induced hazard data layers (keys are 'pga', + etc., values are 2D arrays of data). + geodict (GeoDict): GeoDict object specifying the spatial extent,resolution and shape of the data. - :param origin: + origin (dict): Dictionary with elements: - id String of event ID (i.e., 'us2015abcd') - source String containing originating network ('us') @@ -39,168 +39,193 @@ def __init__(self,layers,geodict,origin,header,metadata=None): - lon Float event longitude - depth Float event depth - magnitude Datetime object representing event origin time. - :param header: + header (header): Dictionary with elements: - type Type of multi-layer earthquake induced hazard ('shakemap','gfe') - version Integer product version (1) - process_time Python datetime indicating when data was created (UTC) - code_version String version of code that created this file (i.e.,'4.0') - originator String representing network that created the hazard grid. - - product_id String representing the ID of the product (may be different from origin ID) + - product_id String representing the ID of the product (may be different + from origin ID) - map_status String, one of RELEASED, ?? - event_type String, one of ['ACTUAL','SCENARIO'] - :param metadata: + metadata (dict): Dictionary of dictionaries containing other metadata users wish to preserve. - :returns: A MultiHazardGrid object. + Returns: A MultiHazardGrid object. """ self._layers = collections.OrderedDict() self._geodict = geodict - for (layerkey,layerdata) in layers.items(): - self._layers[layerkey] = Grid2D(data=layerdata,geodict=geodict) - + for (layerkey, layerdata) in layers.items(): + self._layers[layerkey] = Grid2D(data=layerdata, geodict=geodict) + self.setHeader(header) self.setOrigin(origin) self.setMetadata(metadata) @classmethod - def getFileGeoDict(cls,filename): - """Get the spatial extent, resolution, and shape of grid inside MultiHazardGrid format HDF file. - - :param filename: - File name of HDF file. - :returns: - GeoDict specifying spatial extent, resolution, and shape of grid inside HDF file. + def getFileGeoDict(cls, filename): + """Get the spatial extent, resolution, and shape of grid inside MultiHazardGrid + format HDF file. + + Args: + filename (str): + File name of HDF file. + Returns: + GeoDict specifying spatial extent, resolution, and shape of grid inside + HDF file. """ f = h5py.File(filename, "r") geodict = {} - xvar = f['x'][:] - yvar = f['y'][:] - geodict['xmin'] = xvar[0] - geodict['xmax'] = xvar[-1] - geodict['ymin'] = yvar[0] - geodict['ymax'] = yvar[-1] - geodict['ny'] = len(yvar) - geodict['nx'] = len(xvar) - geodict['dx'] = xvar[1]-xvar[0] - geodict['dy'] = yvar[1]-yvar[0] + xvar = f["x"][:] + yvar = f["y"][:] + geodict["xmin"] = xvar[0] + geodict["xmax"] = xvar[-1] + geodict["ymin"] = yvar[0] + geodict["ymax"] = yvar[-1] + geodict["ny"] = len(yvar) + geodict["nx"] = len(xvar) + geodict["dx"] = xvar[1] - xvar[0] + geodict["dy"] = yvar[1] - yvar[0] f.close() return GeoDict(geodict) - - def _saveDict(self,group,mydict): + + def _saveDict(self, group, mydict): """ Recursively save dictionaries into groups in an HDF file. - :param group: - HDF group object to contain a given dictionary of data in HDF file. - :param mydict: - Dictionary of values to save in group. Dictionary can contain objects of the following types: - - str,unicode,int,float,long,list,tuple,np.ndarray,dict,datetime.datetime,collections.OrderedDict + + Args: + group (HDF): + HDF group object to contain a given dictionary of data in HDF file. + mydict (dict): + Dictionary of values to save in group. Dictionary can contain objects + of the following types: + - str,unicode,int,float,long,list,tuple,np.ndarray,dict, + datetime.datetime,collections.OrderedDict """ - ALLOWED = [str,int,float, - type(None), - list,tuple,np.ndarray, - dict,datetime.datetime, - collections.OrderedDict] + ALLOWED = [ + str, + int, + float, + type(None), + list, + tuple, + np.ndarray, + dict, + datetime.datetime, + collections.OrderedDict, + ] if sys.version_info.major == 2: - ALLOWED.append(unicode) - ALLOWED.append(long) - - for (key,value) in mydict.items(): + ALLOWED.append(str) + ALLOWED.append(int) + + for (key, value) in mydict.items(): tvalue = type(value) if tvalue not in ALLOWED: raise DataSetException('Unsupported metadata value type "%s"' % tvalue) - if not isinstance(value,dict): - if isinstance(value,datetime.datetime): + if not isinstance(value, dict): + if isinstance(value, datetime.datetime): value = time.mktime(value.timetuple()) group.attrs[key] = value else: subgroup = group.create_group(key) - self._saveDict(subgroup,value) + self._saveDict(subgroup, value) @classmethod - def _loadDict(cls,group): + def _loadDict(cls, group): """Recursively load dictionaries from groups in an HDF file. - - :param group: - HDF5 group object. - :returns: - Dictionary of metadata (possibly containing other dictionaries). + + Args: + group (HDF5): + HDF5 group object. + Returns: + Dictionary of metadata (possibly containing other dictionaries). """ tdict = {} - for (key,value) in group.attrs.items(): #attrs are NOT subgroups - if key.find('time') > -1: + for (key, value) in group.attrs.items(): # attrs are NOT subgroups + if key.find("time") > -1: value = value = datetime.datetime.utcfromtimestamp(value) tdict[key] = value - for (key,value) in group.items(): #these are going to be the subgroups + for (key, value) in group.items(): # these are going to be the subgroups tdict[key] = cls._loadDict(value) return tdict - - - def save(self,filename): + + def save(self, filename): """ - Save MultiHazardGrid object to HDF file. - Georeferencing information will be saved as datasets "x" and "y". Layers will be saved as - datasets named by layer keys. Dictionaries contained in "origin", and "header" will be saved in - groups of those same names. Dictionaries contained in the "metadata" dictionary will be contained - in a series of recursive groups under a group called "metadata". - - :param filename: - Output desired filename (HDF format). + Save MultiHazardGrid object to HDF file. + Georeferencing information will be saved as datasets "x" and "y". Layers will + be saved as datasets named by layer keys. Dictionaries contained in "origin", + and "header" will be saved in groups of those same names. Dictionaries + contained in the "metadata" dictionary will be contained in a series of + recursive groups under a group called "metadata". + + Args: + filename (str): + Output desired filename (HDF format). """ f = h5py.File(filename, "w") - #Add in some attributes that will help make this GMT friendly... - f.attrs['Conventions'] = 'COARDS, CF-1.5' - f.attrs['title'] = 'filename' - f.attrs['history'] = 'Created with python MultiHazardGrid.save(%s)' % filename - f.attrs['GMT_version'] = 'NA' - - #create two top-level groups that should always be present - header = f.create_group('header') - self._saveDict(header,self._header) - - origin = f.create_group('origin') - self._saveDict(origin,self._origin) - - #write out any other metadata, creating groups recursively as needed - metadata = f.create_group('metadata') - self._saveDict(metadata,self._metadata) - - xvar = np.linspace(self._geodict.xmin,self._geodict.xmax,self._geodict.nx) - yvar = np.linspace(self._geodict.ymin,self._geodict.ymax,self._geodict.ny) - x = f.create_dataset('x',data=xvar,compression='gzip',shape=xvar.shape,dtype=str(xvar.dtype)) - x.attrs['CLASS'] = 'DIMENSION_SCALE' - x.attrs['NAME'] = 'x' - x.attrs['_Netcdf4Dimid'] = 0 #no idea what this is - x.attrs['long_name'] = 'x' - x.attrs['actual_range'] = np.array((xvar[0],xvar[-1])) - - y = f.create_dataset('y',data=yvar,compression='gzip',shape=yvar.shape,dtype=str(yvar.dtype)) - y.attrs['CLASS'] = 'DIMENSION_SCALE' - y.attrs['NAME'] = 'y' - y.attrs['_Netcdf4Dimid'] = 1 #no idea what this is - y.attrs['long_name'] = 'y' - y.attrs['actual_range'] = np.array((yvar[0],yvar[-1])) - - for (layerkey,layer) in self._layers.items(): - dset = f.create_dataset(layerkey,data=layer.getData(),compression='gzip') - dset.attrs['long_name'] = layerkey - dset.attrs['actual_range'] = np.array((np.nanmin(layer._data),np.nanmax(layer._data))) - + # Add in some attributes that will help make this GMT friendly... + f.attrs["Conventions"] = "COARDS, CF-1.5" + f.attrs["title"] = "filename" + f.attrs["history"] = "Created with python MultiHazardGrid.save(%s)" % filename + f.attrs["GMT_version"] = "NA" + + # create two top-level groups that should always be present + header = f.create_group("header") + self._saveDict(header, self._header) + + origin = f.create_group("origin") + self._saveDict(origin, self._origin) + + # write out any other metadata, creating groups recursively as needed + metadata = f.create_group("metadata") + self._saveDict(metadata, self._metadata) + + xvar = np.linspace(self._geodict.xmin, self._geodict.xmax, self._geodict.nx) + yvar = np.linspace(self._geodict.ymin, self._geodict.ymax, self._geodict.ny) + x = f.create_dataset( + "x", data=xvar, compression="gzip", shape=xvar.shape, dtype=str(xvar.dtype) + ) + x.attrs["CLASS"] = "DIMENSION_SCALE" + x.attrs["NAME"] = "x" + x.attrs["_Netcdf4Dimid"] = 0 # no idea what this is + x.attrs["long_name"] = "x" + x.attrs["actual_range"] = np.array((xvar[0], xvar[-1])) + + y = f.create_dataset( + "y", data=yvar, compression="gzip", shape=yvar.shape, dtype=str(yvar.dtype) + ) + y.attrs["CLASS"] = "DIMENSION_SCALE" + y.attrs["NAME"] = "y" + y.attrs["_Netcdf4Dimid"] = 1 # no idea what this is + y.attrs["long_name"] = "y" + y.attrs["actual_range"] = np.array((yvar[0], yvar[-1])) + + for (layerkey, layer) in self._layers.items(): + dset = f.create_dataset(layerkey, data=layer.getData(), compression="gzip") + dset.attrs["long_name"] = layerkey + dset.attrs["actual_range"] = np.array( + (np.nanmin(layer._data), np.nanmax(layer._data)) + ) + f.close() @classmethod - def load(cls,filename): + def load(cls, filename): """ Load data from an HDF file into a MultiHazardGrid object. - - :param filename: - HDF file containing data and metadata for ShakeMap or Secondary Hazards data. - :returns: - MultiHazardGrid object. + + Args: + filename (str): + HDF file containing data and metadata for ShakeMap or Secondary Hazards + data. + Returns: + MultiHazardGrid object. """ f = h5py.File(filename, "r") - REQUIRED_GROUPS = ['origin','header'] - REQUIRED_DATASETS = ['x','y'] + REQUIRED_GROUPS = ["origin", "header"] + REQUIRED_DATASETS = ["x", "y"] for group in REQUIRED_GROUPS: if group not in f.keys(): raise DataSetException('Missing required metadata group "%s"' % group) @@ -209,151 +234,164 @@ def load(cls,filename): raise DataSetException('Missing required data set "%s"' % dset) header = {} - for (key,value) in f['header'].attrs.items(): - if key.find('time') > -1: + for (key, value) in f["header"].attrs.items(): + if key.find("time") > -1: value = datetime.datetime.utcfromtimestamp(value) header[key] = value origin = {} - for (key,value) in f['origin'].attrs.items(): - if key.find('time') > -1: + for (key, value) in f["origin"].attrs.items(): + if key.find("time") > -1: value = datetime.datetime.utcfromtimestamp(value) origin[key] = value - if 'metadata' in f.keys(): - metadata = cls._loadDict(f['metadata']) + if "metadata" in f.keys(): + metadata = cls._loadDict(f["metadata"]) geodict = {} - xvar = f['x'][:] - yvar = f['y'][:] - geodict['xmin'] = xvar[0] - geodict['xmax'] = xvar[-1] - geodict['ymin'] = yvar[0] - geodict['ymax'] = yvar[-1] - geodict['ny'] = len(yvar) - geodict['nx'] = len(xvar) - geodict['dx'] = xvar[1]-xvar[0] - geodict['dy'] = yvar[1]-yvar[0] + xvar = f["x"][:] + yvar = f["y"][:] + geodict["xmin"] = xvar[0] + geodict["xmax"] = xvar[-1] + geodict["ymin"] = yvar[0] + geodict["ymax"] = yvar[-1] + geodict["ny"] = len(yvar) + geodict["nx"] = len(xvar) + geodict["dx"] = xvar[1] - xvar[0] + geodict["dy"] = yvar[1] - yvar[0] gd = GeoDict(geodict) layers = collections.OrderedDict() - dictDict = {} for key in f.keys(): keytype = str(type(f[key])) - if keytype.find('Dataset') > -1: + if keytype.find("Dataset") > -1: if key in REQUIRED_DATASETS: continue layers[key] = f[key][:] f.close() - return cls(layers,gd,origin,header,metadata=metadata) + return cls(layers, gd, origin, header, metadata=metadata) - def setHeader(self,header): + def setHeader(self, header): """ Set the header dictionary. - - :param header: - Dictionary with elements: - - type Type of multi-layer earthquake induced hazard ('shakemap','gfe') - - version Integer product version (1) - - process_time Python datetime indicating when data was created. - - code_version String version of code that created this file (i.e.,'4.0') - - originator String representing network that created the hazard grid. - - product_id String representing the ID of the product (may be different from origin ID) - - map_status String, one of RELEASED, ?? - - event_type String, one of ['ACTUAL','SCENARIO'] + + Args: + header (dict): + Dictionary with elements: + - type Type of multi-layer earthquake induced hazard + ('shakemap','gfe') + - version Integer product version (1) + - process_time Python datetime indicating when data was created. + - code_version String version of code that created this file + (i.e.,'4.0') + - originator String representing network that created the hazard + grid. + - product_id String representing the ID of the product + (may be different from origin ID) + - map_status String, one of RELEASED, ?? + - event_type String, one of ['ACTUAL','SCENARIO'] """ - self._header = header.copy() #validate later + self._header = header.copy() # validate later - def setOrigin(self,origin): + def setOrigin(self, origin): """ Set the origin dictionary. - - Dictionary with elements: - - id String of event ID (i.e., 'us2015abcd') - - source String containing originating network ('us') - - time Float event magnitude - - lat Float event latitude - - lon Float event longitude - - depth Float event depth - - magnitude Datetime object representing event origin time. + + Args: + origin (dict): + Dictionary with elements: + - id String of event ID (i.e., 'us2015abcd') + - source String containing originating network ('us') + - time Float event magnitude + - lat Float event latitude + - lon Float event longitude + - depth Float event depth + - magnitude Datetime object representing event origin time. """ - self._origin = origin.copy() #validate later + self._origin = origin.copy() # validate later - def setMetadata(self,metadata): + def setMetadata(self, metadata): """ Set the metadata dictionary. - - :param metadata: - Dictionary of dictionaries of metadata. Each dictionary can contain any of the following types: - str,unicode,int,float,long,list,tuple,np.ndarray,dict,datetime.datetime,collections.OrderedDict. + + Args: + metadata (dict): + Dictionary of dictionaries of metadata. Each dictionary can contain + any of the following types: str,unicode,int,float,long,list,tuple, + np.ndarray,dict,datetime.datetime,collections.OrderedDict. """ self._metadata = metadata.copy() - + def getHeader(self): """ Return the header dictionary. - - :returns: - Header dictionary (see setHeader()). + + Returns: + Header dictionary (see setHeader()). """ return self._header.copy() def getOrigin(self): """ Return the origin dictionary. - - :returns: - Origin dictionary (see setOrigin()). + + Returns: + Origin dictionary (see setOrigin()). """ return self._origin.copy() def getMetadata(self): """ Return the dictionary of arbitrary metadata dictionaries. - - :returns: + + Returns: A dictionary of dictionaries containing arbitrary metadata. """ return self._metadata.copy() -if __name__ == '__main__': + +if __name__ == "__main__": shakefile = sys.argv[1] t1 = datetime.datetime.now() sgrid = ShakeGrid.load(shakefile) t2 = datetime.datetime.now() origin = {} - origin['id'] = sgrid._eventDict['event_id'] - origin['source'] = sgrid._eventDict['event_network'] - origin['time'] = sgrid._eventDict['event_timestamp'] - origin['lat'] = sgrid._eventDict['lat'] - origin['lon'] = sgrid._eventDict['lon'] - origin['depth'] = sgrid._eventDict['depth'] - origin['magnitude'] = sgrid._eventDict['magnitude'] + origin["id"] = sgrid._eventDict["event_id"] + origin["source"] = sgrid._eventDict["event_network"] + origin["time"] = sgrid._eventDict["event_timestamp"] + origin["lat"] = sgrid._eventDict["lat"] + origin["lon"] = sgrid._eventDict["lon"] + origin["depth"] = sgrid._eventDict["depth"] + origin["magnitude"] = sgrid._eventDict["magnitude"] header = {} - header['type'] = 'shakemap' - header['version'] = sgrid._shakeDict['shakemap_version'] - header['process_time'] = sgrid._shakeDict['process_timestamp'] - header['code_version'] = sgrid._shakeDict['code_version'] - header['originator'] = sgrid._shakeDict['shakemap_originator'] - header['product_id'] = sgrid._shakeDict['shakemap_id'] - header['map_status'] = sgrid._shakeDict['map_status'] - header['event_type'] = sgrid._shakeDict['shakemap_event_type'] + header["type"] = "shakemap" + header["version"] = sgrid._shakeDict["shakemap_version"] + header["process_time"] = sgrid._shakeDict["process_timestamp"] + header["code_version"] = sgrid._shakeDict["code_version"] + header["originator"] = sgrid._shakeDict["shakemap_originator"] + header["product_id"] = sgrid._shakeDict["shakemap_id"] + header["map_status"] = sgrid._shakeDict["map_status"] + header["event_type"] = sgrid._shakeDict["shakemap_event_type"] layers = collections.OrderedDict() - for (layername,layerdata) in sgrid.getData().items(): + for (layername, layerdata) in sgrid.getData().items(): layers[layername] = layerdata.getData() - tdict = {'name':'fred','family':{'wife':'wilma','daughter':'pebbles'}} - mgrid = MultiHazardGrid(layers,sgrid.getGeoDict(),origin,header,metadata={'flintstones':tdict}) - mgrid.save('test.hdf') + tdict = {"name": "fred", "family": {"wife": "wilma", "daughter": "pebbles"}} + mgrid = MultiHazardGrid( + layers, sgrid.getGeoDict(), origin, header, metadata={"flintstones": tdict} + ) + mgrid.save("test.hdf") t3 = datetime.datetime.now() - mgrid2 = MultiHazardGrid.load('test.hdf') + mgrid2 = MultiHazardGrid.load("test.hdf") t4 = datetime.datetime.now() - xmlmb = os.path.getsize(shakefile)/float(1e6) - hdfmb = os.path.getsize('test.hdf')/float(1e6) - xmltime = (t2-t1).seconds + (t2-t1).microseconds/float(1e6) - hdftime = (t4-t3).seconds + (t4-t3).microseconds/float(1e6) - print('Input XML file size: %.2f MB (loading time %.3f seconds)' % (xmlmb,xmltime)) - print('Output HDF file size: %.2f MB (loading time %.3f seconds)' % (hdfmb,hdftime)) - os.remove('test.hdf') + xmlmb = os.path.getsize(shakefile) / float(1e6) + hdfmb = os.path.getsize("test.hdf") / float(1e6) + xmltime = (t2 - t1).seconds + (t2 - t1).microseconds / float(1e6) + hdftime = (t4 - t3).seconds + (t4 - t3).microseconds / float(1e6) + print("Input XML file size: %.2f MB (loading time %.3f seconds)" % (xmlmb, xmltime)) + print( + "Output HDF file size: %.2f MB (loading time %.3f seconds)" % (hdfmb, hdftime) + ) + os.remove("test.hdf") diff --git a/mapio/multiple.py b/mapio/multiple.py index c756cb8..a8088c5 100755 --- a/mapio/multiple.py +++ b/mapio/multiple.py @@ -1,5 +1,4 @@ - -#python 3 compatibility +# python 3 compatibility from __future__ import print_function from .gridbase import Grid @@ -8,32 +7,38 @@ import abc from collections import OrderedDict + class MultiGrid(Grid): - def __init__(self,layers,descriptions=None): - """ - Construct a semi-abstract MultiGrid object, which can contain many 2D layers of gridded data, all at the - same resolution and with the same extent. - - :param layers: - OrderedDict of Grid2D objects. - :param descriptions: - list of layer descriptions, or None - :raises DataSetException: - When: - - length of descriptions (when not None) does not match length of layers - - input layers is not an OrderedDict - """ - if not isinstance(layers,OrderedDict): - raise DataSetException('Input layers must be of type OrderedDict.') + def __init__(self, layers, descriptions=None): + """ + Construct a semi-abstract MultiGrid object, which can contain many 2D layers of + gridded data, all at the same resolution and with the same extent. + + Args: + layers (OrderedDict): + OrderedDict of Grid2D objects. + descriptions (list): + list of layer descriptions, or None + Raises: + DataSetException: + When: + - length of descriptions (when not None) does not match length of + layers + - input layers is not an OrderedDict + """ + if not isinstance(layers, OrderedDict): + raise DataSetException("Input layers must be of type OrderedDict.") if descriptions is None: - descriptions = ['' for l in layers.keys()] + descriptions = ["" for ll in layers.keys()] if len(descriptions) != len(layers): - raise DataSetException('List of descriptions does not match length of layers.') + raise DataSetException( + "List of descriptions does not match length of layers." + ) lnames = list(layers.keys()) self._layers = OrderedDict() self._descriptions = OrderedDict() - for i in range(0,len(lnames)): + for i in range(0, len(lnames)): layername = lnames[i] layer = layers[layername] desc = descriptions[i] @@ -42,57 +47,64 @@ def __init__(self,layers,descriptions=None): self._descriptions[layername] = desc self._geodict = geodict - @abc.abstractmethod - def save(self,filename): + def save(self, filename): """ Save layers of data to a file - - :param filename: - File to save data to. + + Args: + filename (str): + File to save data to. """ raise NotImplementedError("save method not implemented in MultiGrid") - #subclassed implementation should be a @classmethod + # subclassed implementation should be a @classmethod @abc.abstractmethod - def load(self,filename): + def load(self, filename): """ Load layers of data from a file. - - :param filename: - File to load data from. + + Args: + filename (str): + File to load data from. """ raise NotImplementedError("save method not implemented in MultiGrid") - - def setLayer(self,name,data,desc=None): + + def setLayer(self, name, data, desc=None): """ Add a 2D layer of data to a MultiGrid object. - :param name: - String which will be used to retrieve the data. - :param data: - 2D numpy array of data - :param desc: - Optional text description of layer - :raises DataSetException: - If the data layer dimensions don't match the geodict. + Args: + name (str): + String which will be used to retrieve the data. + data (array): + 2D numpy array of data + desc (:obj:`str`, optional): + Optional text description of layer + Raises: + DataSetException: + If the data layer dimensions don't match the geodict. """ - nr,nc = data.shape + nr, nc = data.shape if nr != self._geodict.ny or nc != self._geodict.nx: - raise DataSetException("Data layer dimensions don't match those already in the grid") - self._layers[name] = Grid2D(data,self._geodict.copy()) + raise DataSetException( + "Data layer dimensions don't match those already in the grid" + ) + self._layers[name] = Grid2D(data, self._geodict.copy()) self._descriptions[name] = desc - def getLayer(self,name): + def getLayer(self, name): """ Retrieve the 2D associated with a layer name. - :param name: - Name of data layer. - :returns: - Grid2D object. - :raises DataSetException: - When name is not found in list of layer names. + Args: + name (str): + Name of data layer. + Returns: + Grid2D object. + Raises: + DataSetException: + When name is not found in list of layer names. """ if name not in self._layers.keys(): raise DataSetException('Layer "%s" not in list of layers.' % name) @@ -102,17 +114,18 @@ def getData(self): """ Return the OrderedDict of data layers contained in MultiGrid. - :returns: - OrderedDict of Grid2D objects. + Returns: + OrderedDict of Grid2D objects. """ return self._layers - def setData(self,layers,descriptions=None): + def setData(self, layers, descriptions=None): """ Return the OrderedDict of data layers contained in MultiGrid. - :param layers: - OrderedDict of Grid2D objects. + Args: + layers (OrderedDict): + OrderedDict of Grid2D objects. """ self._layers = layers layernames = layers.keys() @@ -120,135 +133,160 @@ def setData(self,layers,descriptions=None): def getGeoDict(self): """ - Return the geodict object which defines the extent and resolution of all the grids. - - :returns: - geodict dictionary (see constructor) + Return the geodict object which defines the extent and resolution of all the + grids. + + Returns: + geodict dictionary (see constructor) """ return self._geodict def getBounds(self): """ Return the lat/lon range of the data. - - :returns: - Tuple of (lonmin,lonmax,latmin,latmax) - """ - return (self._geodict.xmin,self._geodict.xmax,self._geodict.ymin,self._geodict.ymax) - - def trim(self,geodict,resample=False,method='linear'): - """ - Trim all data layers to a smaller set of bounds, resampling if requested. If not resampling, - data will be trimmed to smallest grid boundary possible. - - :param geodict: - GeoDict used to specify subset bounds and resolution (if resample is selected) - :param resample: - Boolean indicating whether the data should be resampled to *exactly* match input bounds. - :param method: - If resampling, method used, one of ('linear','nearest','cubic','quintic') - """ - for (layername,layer) in self._layers.items(): - layer.trim(geodict,resample=resample,method=method) + + Returns: + Tuple of (lonmin,lonmax,latmin,latmax) + """ + return ( + self._geodict.xmin, + self._geodict.xmax, + self._geodict.ymin, + self._geodict.ymax, + ) + + def trim(self, geodict, resample=False, method="linear"): + """ + Trim all data layers to a smaller set of bounds, resampling if requested. If + not resampling, data will be trimmed to smallest grid boundary possible. + + Args: + geodict (geodict): + GeoDict used to specify subset bounds and resolution (if resample is + selected) + resample: (bool) + Boolean indicating whether the data should be resampled to *exactly* + match input bounds. + method (:obj:`str`, optional): + If resampling, method used, one of ('linear','nearest','cubic','quintic' + ) + """ + for (layername, layer) in self._layers.items(): + layer.trim(geodict, resample=resample, method=method) self._geodict = layer.getGeoDict().copy() def getLayerNames(self): """ Return the list of layer names contained in the MultiGrid. - - :returns: - List of layer names. + + Returns: + List of layer names. """ return self._layers.keys() - - def getValue(self,lat,lon,layername,method='nearest',default=None): + + def getValue(self, lat, lon, layername, method="nearest", default=None): """Return numpy array at given latitude and longitude (using nearest neighbor). - - :param lat: - Latitude (in decimal degrees) of desired data value. - :param lon: - Longitude (in decimal degrees) of desired data value. - :param layername: - Name of layer from which to retrieve data. - :param method: - Interpolation method, one of ('nearest','linear','cubic','quintic') - :param default: - Default value to return when lat/lon is outside of grid bounds. - :return: + + Args: + lat (float): + Latitude (in decimal degrees) of desired data value. + lon (float): + Longitude (in decimal degrees) of desired data value. + layername (str): + Name of layer from which to retrieve data. + method (:obj:`str`, optional): + Interpolation method, one of ('nearest','linear','cubic','quintic') + default (:obj:`float`, optional): + Default value to return when lat/lon is outside of grid bounds. + Return: Value at input latitude,longitude position. """ - return self._layers[layername].getValue(lat,lon,method=method,default=default) + return self._layers[layername].getValue( + lat, lon, method=method, default=default + ) + + def getLatLon(self, row, col): + """Return geographic coordinates (lat/lon decimal degrees) for given data row + and column. - 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. + Args: + row (int): + Row dimension index into internal data array. + col (int): + Column dimension index into internal data array. + Returns: + Tuple of latitude and longitude. """ layernames = self._layers.keys() - return self._layers[layernames[0]].getLatLon(row,col) - - 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: - Input longitude. - :param returnFloat: - Boolean indicating whether floating point row/col coordinates should be returned. - :returns: - Tuple of row and column. + return self._layers[layernames[0]].getLatLon(row, col) + + def getRowCol(self, lat, lon, returnFloat=False): + """Return data row and column from given geographic coordinates (lat/lon + decimal degrees). + + Args: + lat (float): + Input latitude. + lon (float): + Input longitude. + returnFloat (bool): + Boolean indicating whether floating point row/col coordinates should be + returned. + Returns: + Tuple of row and column. """ layernames = self._layers.keys() - return self._layers[layernames[0]].getRowCol(lat,lon,returnFloat=returnFloat) + return self._layers[layernames[0]].getRowCol(lat, lon, returnFloat=returnFloat) - def subdivide(self,finerdict,cellFill='max'): + def subdivide(self, finerdict, cellFill="max"): """Subdivide the cells of the host grid into finer-resolution cells. - :param finerdict: - GeoDict object defining a grid with a finer resolution than the host grid. - :param cellFill: - String defining how to fill cells that span more than one host grid cell. - Choices are: - 'max': Choose maximum value of host grid cells. - 'min': Choose minimum value of host grid cells. - 'mean': Choose mean value of host grid cells. - :returns: - MultiGrid instance with host grid values subdivided onto finer grid. - :raises DataSetException: - When finerdict is not a) finer resolution or b) does not intersect.x or cellFill is not valid. + Args: + finerdict (geodict): + GeoDict object defining a grid with a finer resolution than the host + grid. + cellFill (:obj:`str`, optional): + String defining how to fill cells that span more than one host grid + cell. + Choices are: + 'max': Choose maximum value of host grid cells. + 'min': Choose minimum value of host grid cells. + 'mean': Choose mean value of host grid cells. + Returns: + MultiGrid instance with host grid values subdivided onto finer grid. + Raises: + DataSetException: + When finerdict is not a) finer resolution or b) does not intersect.x or + cellFill is not valid. """ layers = OrderedDict() - for (layername,layer) in self._layers.items(): - layers[layername] = layer.subdivide(finerdict,cellFill=cellFill) + for (layername, layer) in self._layers.items(): + layers[layername] = layer.subdivide(finerdict, cellFill=cellFill) return MultiGrid(layers) - - def interpolateToGrid(self,geodict,method='linear'): - """ - Given a geodict specifying another grid extent and resolution, resample all grids to match. - - :param geodict: - geodict dictionary from another grid whose extents are inside the extent of this grid. - :param method: - Optional interpolation method - ['linear', 'cubic','quintic','nearest'] - :raises DataSetException: - If the Grid object upon which this function is being called is not completely contained by the - grid to which this Grid is being resampled. - :raises DataSetException: - If the resulting interpolated grid shape does not match input geodict. + + def interpolateToGrid(self, geodict, method="linear"): + """ + Given a geodict specifying another grid extent and resolution, resample all + grids to match. + + Arggs: + geodict (geodict): + geodict dictionary from another grid whose extents are inside the + extent of this grid. + method (:obj:`str`, optional): + Optional interpolation method - ['linear', 'cubic','quintic','nearest'] + Raises: + DataSetException: + If the Grid object upon which this function is being called is not + completely contained by the grid to which this Grid is being resampled. + DataSetException: + If the resulting interpolated grid shape does not match input geodict. This function modifies the internal griddata and geodict object variables. """ layers = OrderedDict() - for (layername,layer) in self._layers.items(): - #layer.interpolateToGrid(geodict,method=method) - layers[layername] = layer.interpolateToGrid(geodict,method=method) - #self._geodict = layer.getGeoDict().copy() + for (layername, layer) in self._layers.items(): + # layer.interpolateToGrid(geodict,method=method) + layers[layername] = layer.interpolateToGrid(geodict, method=method) + # self._geodict = layer.getGeoDict().copy() return MultiGrid(layers) - -