From 644a70e5fb37cf409de11ca0104e4e7ce8c5e674 Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Mon, 8 Nov 2021 08:10:46 -0700 Subject: [PATCH 1/3] Added check in read() function to return grid with all pad value pixels in it if the sample dict bounds do not intersect with the file bounds. --- mapio/reader.py | 102 ++++++++++++++------------- test/reader_test.py | 163 ++++++++++++++++++++++++++------------------ 2 files changed, 154 insertions(+), 111 deletions(-) diff --git a/mapio/reader.py b/mapio/reader.py index 55a9e33..7c7e575 100644 --- a/mapio/reader.py +++ b/mapio/reader.py @@ -17,12 +17,12 @@ def _is_hdf(filename): bool: True if HDF, False if not. """ is_hdf = False - with open(filename, 'rb') as f: + with open(filename, "rb") as f: f.seek(1) bytes = header = f.read(3) try: header = bytes.decode() - if header == 'HDF': + if header == "HDF": is_hdf = True except UnicodeDecodeError: pass @@ -45,15 +45,15 @@ def _get_geodict_from_window(affine, window, data): """ xmin, ymax = affine * (window.col_off, window.row_off) geodict = {} - geodict['dx'] = affine.a - geodict['dy'] = -1 * affine.e - geodict['xmin'] = xmin + geodict['dx'] / 2.0 - geodict['ymax'] = ymax - geodict['dy'] / 2.0 + geodict["dx"] = affine.a + geodict["dy"] = -1 * affine.e + geodict["xmin"] = xmin + geodict["dx"] / 2.0 + geodict["ymax"] = ymax - geodict["dy"] / 2.0 nrows, ncols = data.shape - geodict['ny'] = nrows - geodict['nx'] = ncols - geodict['xmax'] = geodict['xmin'] + (geodict['nx'] - 1) * geodict['dx'] - geodict['ymin'] = geodict['ymax'] - (geodict['ny'] - 1) * geodict['dy'] + geodict["ny"] = nrows + geodict["nx"] = ncols + geodict["xmax"] = geodict["xmin"] + (geodict["nx"] - 1) * geodict["dx"] + geodict["ymin"] = geodict["ymax"] - (geodict["ny"] - 1) * geodict["dy"] gd = GeoDict(geodict) return gd @@ -154,20 +154,20 @@ def _read_pixels(src, window): """ fname = src.files[0] is_hdf = _is_hdf(fname) - if src.driver != 'netCDF' or not is_hdf: + if src.driver != "netCDF" or not is_hdf: data = np.squeeze(src.read(window=window), axis=0) else: src.close() - f = h5py.File(fname, 'r') + f = h5py.File(fname, "r") if window is None: - data = np.flipud(f['z']) + data = np.flipud(f["z"]) else: cstart = int(window.col_off) cend = cstart + int(window.width) rend = src.height - int(window.row_off) rstart = rend - int(window.height) # t1 = time.time() - data = np.flipud(f['z'][rstart:rend, cstart:cend]) + data = np.flipud(f["z"][rstart:rend, cstart:cend]) # t2 = time.time() # print('h5py read: %.3f seconds.' % (t2-t1)) f.close() @@ -201,9 +201,7 @@ def _read_data(src, samplegeodict, resample, method): if resample: pad = True if not is_edge: - window = _geodict_to_window(samplegeodict, - src, - pad=pad) + window = _geodict_to_window(samplegeodict, src, pad=pad) data, src = _read_pixels(src, window) gd = _get_geodict_from_window(affine, window, data) @@ -219,15 +217,13 @@ def _read_data(src, samplegeodict, resample, method): lxmax += dx / 2.0 ymin = samplegeodict.ymin ymax = samplegeodict.ymax - sample_left = GeoDict.createDictFromBox(lxmin, lxmax, - ymin, ymax, dx, dy) + sample_left = GeoDict.createDictFromBox(lxmin, lxmax, ymin, ymax, dx, dy) # get the left edge of the grid rxmin, _ = affine * (0, 0) # convert to pixel registered rxmin += dx / 2.0 rxmax = samplegeodict.xmax - sample_right = GeoDict.createDictFromBox(rxmin, rxmax, - ymin, ymax, dx, dy) + sample_right = GeoDict.createDictFromBox(rxmin, rxmax, ymin, ymax, dx, dy) left_window = _geodict_to_window(sample_left, src, pad=pad) left_block, src = _read_pixels(src, left_window) @@ -242,29 +238,24 @@ def _read_data(src, samplegeodict, resample, method): row_off = right_window.row_off width = right_window.width - dwidth height = right_window.height - right_window = rasterio.windows.Window(col_off, - row_off, - width, - height) + right_window = rasterio.windows.Window(col_off, row_off, width, height) right_block, src = _read_pixels(src, right_window) - left_gd = _get_geodict_from_window(affine, - left_window, - left_block) + left_gd = _get_geodict_from_window(affine, left_window, left_block) - right_gd = _get_geodict_from_window(affine, - right_window, - right_block) + right_gd = _get_geodict_from_window(affine, right_window, right_block) data = np.concatenate((left_block, right_block), axis=1) nrows, ncols = data.shape - geodict = {'xmin': left_gd.xmin, - 'xmax': right_gd.xmax, - 'ymin': ymin, - 'ymax': ymax, - 'dx': dx, - 'dy': dy, - 'nx': ncols, - 'ny': nrows} + geodict = { + "xmin": left_gd.xmin, + "xmax": right_gd.xmax, + "ymin": ymin, + "ymax": ymax, + "dx": dx, + "dy": dy, + "nx": ncols, + "ny": nrows, + } gd = GeoDict(geodict) grid = Grid2D(data, gd) return grid @@ -283,9 +274,17 @@ def get_file_geodict(filename): return gd -def read(filename, samplegeodict=None, resample=False, - method='linear', doPadding=False, padValue=np.nan, - apply_nan=True, force_cast=True, interp_approach='scipy'): +def read( + filename, + samplegeodict=None, + resample=False, + method="linear", + doPadding=False, + padValue=np.nan, + apply_nan=True, + force_cast=True, + interp_approach="scipy", +): """Read part or all of a rasterio file, resampling and padding as necessary. If samplegeodict is not provided, then the entire file will be read. @@ -312,6 +311,8 @@ def read(filename, samplegeodict=None, resample=False, apply_nan (bool): Convert nodata values to NaNs, upcasting to float if necessary. force_cast (bool): If data values exceed range of values, cast upward. interp_approach (str): One of 'scipy', 'rasterio'. + Returns: + Grid2D: Object containing desired data and a Geodict matching samplegeodict. """ # use rasterio to read all formats src = rasterio.open(filename) @@ -325,6 +326,16 @@ def read(filename, samplegeodict=None, resample=False, src.close() return grid + # shortcut out here to return NaN grid if samplegeodict is completely outside the + # source data grid. + fdict = get_file_geodict(filename) + if not fdict.intersects(samplegeodict): + nx = samplegeodict.nx + ny = samplegeodict.ny + data = np.ones((ny, nx)) * padValue + grid = Grid2D(data=data, geodict=samplegeodict) + return grid + # if non-nearest resampling, this grid may have a ring of padding pixels # around the outside. grid = _read_data(src, samplegeodict, resample, method) @@ -333,7 +344,7 @@ def read(filename, samplegeodict=None, resample=False, # if padding is turned off rdict = grid.getGeoDict() if resample: - if method == 'nearest': + if method == "nearest": c1 = rdict.xmin <= samplegeodict.xmin c2 = rdict.xmax >= samplegeodict.xmax c3 = rdict.ymin <= samplegeodict.ymin @@ -346,8 +357,7 @@ def read(filename, samplegeodict=None, resample=False, c3 = rdict.ymin > samplegeodict.ymin - dy c4 = rdict.ymax < samplegeodict.ymax + dy if (c1 or c2 or c3 or c4) and not doPadding: - msg = ('Without padding you cannot ask to resample from ' - 'edge of grid.') + msg = "Without padding you cannot ask to resample from " "edge of grid." raise IndexError(msg) grid._geodict.nodata = src.nodata @@ -365,7 +375,7 @@ def read(filename, samplegeodict=None, resample=False, grid = Grid2D(data, gd) if resample: - if interp_approach == 'scipy': + if interp_approach == "scipy": grid = grid.interpolateToGrid(samplegeodict, method=method) else: grid = grid.interpolate2(samplegeodict, method=method) diff --git a/test/reader_test.py b/test/reader_test.py index 3a45096..818a95c 100755 --- a/test/reader_test.py +++ b/test/reader_test.py @@ -12,31 +12,32 @@ def test_read_whole(): - files = {'ESRI Float': 'samplegrid_flt.flt', - 'NetCDF 3': 'samplegrid_cdf.cdf'} + files = {"ESRI Float": "samplegrid_flt.flt", "NetCDF 3": "samplegrid_cdf.cdf"} # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file for ftype, fname in files.items(): - datafile = os.path.join(homedir, 'data', fname) + datafile = os.path.join(homedir, "data", fname) grid = read(datafile) assert grid._geodict.xmin == 5.0 - print('Successful read of %s' % ftype) + print("Successful read of %s" % ftype) def test_read_subset_no_resample(): # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file - datafile = os.path.join(homedir, 'data', 'samplegrid_cdf.cdf') - sdict = {'xmin': 6, - 'xmax': 7, - 'ymin': 5, - 'ymax': 6, - 'nx': 2, - 'ny': 2, - 'dx': 1, - 'dy': 1} + datafile = os.path.join(homedir, "data", "samplegrid_cdf.cdf") + sdict = { + "xmin": 6, + "xmax": 7, + "ymin": 5, + "ymax": 6, + "nx": 2, + "ny": 2, + "dx": 1, + "dy": 1, + } sampledict = GeoDict(sdict) grid = read(datafile, samplegeodict=sampledict) tdata = np.array([[11, 12], [16, 17]]) @@ -48,15 +49,17 @@ def test_resample(): # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file - datafile = os.path.join(homedir, 'data', 'samplegrid_cdf.cdf') - sdict = {'xmin': 6.0, - 'xmax': 7.0, - 'ymin': 6.0, - 'ymax': 7.0, - 'nx': 1, - 'ny': 1, - 'dx': 1, - 'dy': 1} + datafile = os.path.join(homedir, "data", "samplegrid_cdf.cdf") + sdict = { + "xmin": 6.0, + "xmax": 7.0, + "ymin": 6.0, + "ymax": 7.0, + "nx": 1, + "ny": 1, + "dx": 1, + "dy": 1, + } sampledict = GeoDict(sdict) grid = read(datafile, samplegeodict=sampledict, resample=True) tdata = np.array([[6.0]]) @@ -67,20 +70,20 @@ def test_read_subset_with_resample_and_padding(): # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file - datafile = os.path.join(homedir, 'data', 'samplegrid_cdf.cdf') - sdict = {'xmin': 4.5, - 'xmax': 5.5, - 'ymin': 7.5, - 'ymax': 8.5, - 'nx': 2, - 'ny': 2, - 'dx': 1, - 'dy': 1} + datafile = os.path.join(homedir, "data", "samplegrid_cdf.cdf") + sdict = { + "xmin": 4.5, + "xmax": 5.5, + "ymin": 7.5, + "ymax": 8.5, + "nx": 2, + "ny": 2, + "dx": 1, + "dy": 1, + } sampledict = GeoDict(sdict) - grid = read(datafile, samplegeodict=sampledict, - resample=True, doPadding=True) - atest = np.array([[np.nan, np.nan], - [np.nan, 3.]]) + grid = read(datafile, samplegeodict=sampledict, resample=True, doPadding=True) + atest = np.array([[np.nan, np.nan], [np.nan, 3.0]]) np.testing.assert_almost_equal(grid._data, atest) @@ -88,35 +91,61 @@ def test_read_subset_with_padding(): # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file - datafile = os.path.join(homedir, 'data', 'samplegrid_cdf.cdf') - sdict = {'xmin': 4.5, - 'xmax': 5.5, - 'ymin': 7.5, - 'ymax': 8.5, - 'nx': 2, - 'ny': 2, - 'dx': 1, - 'dy': 1} + datafile = os.path.join(homedir, "data", "samplegrid_cdf.cdf") + sdict = { + "xmin": 4.5, + "xmax": 5.5, + "ymin": 7.5, + "ymax": 8.5, + "nx": 2, + "ny": 2, + "dx": 1, + "dy": 1, + } sampledict = GeoDict(sdict) - grid = read(datafile, samplegeodict=sampledict, - resample=False, doPadding=True) + grid = read(datafile, samplegeodict=sampledict, resample=False, doPadding=True) assert grid._data.shape == (3, 3) assert grid._data[1, 1] == 0 +def test_read_outside_lat_range(): + # where is this script? + homedir = os.path.dirname(os.path.abspath(__file__)) + # this is an HDF 5 file + datafile = os.path.join(homedir, "data", "samplegrid_cdf.cdf") + sdict = { + "xmin": 4.5, + "xmax": 5.5, + "ymin": 9.0, + "ymax": 10.0, + "nx": 2, + "ny": 2, + "dx": 1, + "dy": 1, + } + sampledict = GeoDict(sdict) + grid = read(datafile, samplegeodict=sampledict, resample=False, doPadding=True) + # shape matches input + assert grid._data.shape == (sampledict.ny, sampledict.nx) + # data is all NaN + assert np.isnan(grid._data).sum() == 4 + + def test_read_meridian(): # where is this script? homedir = os.path.dirname(os.path.abspath(__file__)) # this is an HDF 5 file - datafile = os.path.join(homedir, 'data', 'globalgrid_cdf.cdf') - sdict = {'xmin': 180, - 'xmax': -120, - 'ymin': 0, - 'ymax': 30, - 'nx': 2, - 'ny': 2, - 'dx': 60, - 'dy': 30} + datafile = os.path.join(homedir, "data", "globalgrid_cdf.cdf") + sdict = { + "xmin": 180, + "xmax": -120, + "ymin": 0, + "ymax": 30, + "nx": 2, + "ny": 2, + "dx": 60, + "dy": 30, + } sampledict = GeoDict(sdict) grid = read(datafile, samplegeodict=sampledict) assert np.nansum(grid._data) == 50.0 @@ -130,22 +159,25 @@ def read_user_file_test(fname, xmin, xmax, ymin, ymax): t2 = time.time() nrows, ncols = grid._data.shape npixels = nrows * ncols - print('%.2f seconds to read %i pixels using h5py' % (t2 - t1, npixels)) - - west, east, south, north = (-105.00416666665, - -102.98750000804999, - 34.98750000805, - 37.00416666665) - src = rasterio.open(fname, 'r') + print("%.2f seconds to read %i pixels using h5py" % (t2 - t1, npixels)) + + west, east, south, north = ( + -105.00416666665, + -102.98750000804999, + 34.98750000805, + 37.00416666665, + ) + src = rasterio.open(fname, "r") window = src.window(west, south, east, north) t1 = time.time() data = src.read(window=window) t2 = time.time() - print('%.2f seconds to read %i pixels using rasterio' % (t2 - t1, npixels)) + print("%.2f seconds to read %i pixels using rasterio" % (t2 - t1, npixels)) ratio = grid._data.sum() / data.sum() - print('Ratio of h5py data to rasterio data is %.4f' % ratio) + print("Ratio of h5py data to rasterio data is %.4f" % ratio) src.close() + # TODO - fix this later # def test_read_approach(): # # where is this script? @@ -173,7 +205,7 @@ def read_user_file_test(fname, xmin, xmax, ymin, ymax): # assert t3 - t2 > t2 - t1 -if __name__ == '__main__': +if __name__ == "__main__": # test_read_approach() test_read_whole() test_read_subset_no_resample() @@ -181,6 +213,7 @@ def read_user_file_test(fname, xmin, xmax, ymin, ymax): test_read_subset_with_padding() test_read_subset_with_resample_and_padding() test_read_meridian() + test_read_outside_lat_range() if len(sys.argv) > 1: fname = sys.argv[1] From da032e2f500002abf78725429c66be11a7fb4724 Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Mon, 8 Nov 2021 08:31:48 -0700 Subject: [PATCH 2/3] scipy nearestNDinterpolator behavior must have changed. --- test/grid2d_test.py | 703 +++++++++++++++++++++++++++--------------- test/multiple_test.py | 150 +++++---- 2 files changed, 543 insertions(+), 310 deletions(-) diff --git a/test/grid2d_test.py b/test/grid2d_test.py index 7db4d6b..085046d 100755 --- a/test/grid2d_test.py +++ b/test/grid2d_test.py @@ -16,7 +16,7 @@ # hack the path so that I can debug these functions if I need to homedir = os.path.dirname(os.path.abspath(__file__)) # where is this script? -mapiodir = os.path.abspath(os.path.join(homedir, '..')) +mapiodir = os.path.abspath(os.path.join(homedir, "..")) # put this at the front of the system path, ignoring any installed mapio stuff sys.path.insert(0, mapiodir) @@ -40,81 +40,141 @@ def test_subdivide(): - print('Testing subdivide method - aligned grids...') + print("Testing subdivide method - aligned grids...") data = np.arange(0, 4).reshape((2, 2)) - geodict = GeoDict({'xmin': 0.0, 'xmax': 1.0, - 'ymin': 0.0, 'ymax': 1.0, - 'dx': 1.0, 'dy': 1.0, - 'ny': 2, 'nx': 2}) + geodict = GeoDict( + { + "xmin": 0.0, + "xmax": 1.0, + "ymin": 0.0, + "ymax": 1.0, + "dx": 1.0, + "dy": 1.0, + "ny": 2, + "nx": 2, + } + ) hostgrid = Grid2D(data, geodict) - finedict = GeoDict({'xmin': 0.0 - (1.0 / 3.0), 'xmax': 1.0 + (1.0 / 3.0), - 'ymin': 0.0 - (1.0 / 3.0), 'ymax': 1.0 + (1.0 / 3.0), - 'dx': 1.0 / 3.0, 'dy': 1.0 / 3.0, - 'ny': 6, 'nx': 6}) + finedict = GeoDict( + { + "xmin": 0.0 - (1.0 / 3.0), + "xmax": 1.0 + (1.0 / 3.0), + "ymin": 0.0 - (1.0 / 3.0), + "ymax": 1.0 + (1.0 / 3.0), + "dx": 1.0 / 3.0, + "dy": 1.0 / 3.0, + "ny": 6, + "nx": 6, + } + ) finegrid = hostgrid.subdivide(finedict) - output = np.array([[0., 0., 0., 1., 1., 1.], - [0., 0., 0., 1., 1., 1.], - [0., 0., 0., 1., 1., 1.], - [2., 2., 2., 3., 3., 3.], - [2., 2., 2., 3., 3., 3.], - [2., 2., 2., 3., 3., 3.]]) + output = np.array( + [ + [0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + [0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + [0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + [2.0, 2.0, 2.0, 3.0, 3.0, 3.0], + [2.0, 2.0, 2.0, 3.0, 3.0, 3.0], + [2.0, 2.0, 2.0, 3.0, 3.0, 3.0], + ] + ) np.testing.assert_almost_equal(finegrid.getData(), output) - print('Passed subdivide method test - aligned grids.') + print("Passed subdivide method test - aligned grids.") - print('Testing subdivide method - non-aligned grids...') + print("Testing subdivide method - non-aligned grids...") data = np.arange(0, 9).reshape((3, 3)) - geodict = GeoDict({'xmin': 0.0, 'xmax': 10.0, - 'ymin': 0.0, 'ymax': 10.0, - 'dx': 5.0, 'dy': 5.0, - 'ny': 3, 'nx': 3}) + geodict = GeoDict( + { + "xmin": 0.0, + "xmax": 10.0, + "ymin": 0.0, + "ymax": 10.0, + "dx": 5.0, + "dy": 5.0, + "ny": 3, + "nx": 3, + } + ) hostgrid = Grid2D(data, geodict) - finedict = GeoDict({'xmin': -2.5, 'xmax': 11.5, - 'ymin': -1.5, 'ymax': 10.5, - 'dx': 2.0, 'dy': 2.0, - 'nx': 8, 'ny': 7}) + finedict = GeoDict( + { + "xmin": -2.5, + "xmax": 11.5, + "ymin": -1.5, + "ymax": 10.5, + "dx": 2.0, + "dy": 2.0, + "nx": 8, + "ny": 7, + } + ) N = np.nan - print('Testing subdivide with min parameter...') - finegrid = hostgrid.subdivide(finedict, cellFill='min') - output = np.array([[N, 0., 0., 1., 1., 1., 2., 2.], - [N, 0., 0., 1., 1., 1., 2., 2.], - [N, 3., 3., 4., 4., 4., 5., 5.], - [N, 3., 3., 4., 4., 4., 5., 5.], - [N, 3., 3., 4., 4., 4., 5., 5.], - [N, 6., 6., 7., 7., 7., 8., 8.], - [N, 6., 6., 7., 7., 7., 8., 8.]]) + print("Testing subdivide with min parameter...") + finegrid = hostgrid.subdivide(finedict, cellFill="min") + output = np.array( + [ + [N, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0], + [N, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 6.0, 6.0, 7.0, 7.0, 7.0, 8.0, 8.0], + [N, 6.0, 6.0, 7.0, 7.0, 7.0, 8.0, 8.0], + ] + ) np.testing.assert_almost_equal(finegrid.getData(), output) - print('Passed subdivide with min parameter...') - print('Testing subdivide with max parameter...') - finegrid = hostgrid.subdivide(finedict, cellFill='max') - output = np.array([[N, 0., 0., 1., 1., 2., 2., 2.], - [N, 0., 0., 1., 1., 2., 2., 2.], - [N, 3., 3., 4., 4., 5., 5., 5.], - [N, 3., 3., 4., 4., 5., 5., 5.], - [N, 6., 6., 7., 7., 8., 8., 8.], - [N, 6., 6., 7., 7., 8., 8., 8.], - [N, 6., 6., 7., 7., 8., 8., 8.]]) + print("Passed subdivide with min parameter...") + print("Testing subdivide with max parameter...") + finegrid = hostgrid.subdivide(finedict, cellFill="max") + output = np.array( + [ + [N, 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 2.0], + [N, 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 2.0], + [N, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 5.0], + [N, 6.0, 6.0, 7.0, 7.0, 8.0, 8.0, 8.0], + [N, 6.0, 6.0, 7.0, 7.0, 8.0, 8.0, 8.0], + [N, 6.0, 6.0, 7.0, 7.0, 8.0, 8.0, 8.0], + ] + ) np.testing.assert_almost_equal(finegrid.getData(), output) - print('Passed subdivide with max parameter...') - print('Testing subdivide with mean parameter...') - finegrid = hostgrid.subdivide(finedict, cellFill='mean') - output = np.array([[N, 0., 0., 1., 1., 1.5, 2., 2.], - [N, 0., 0., 1., 1., 1.5, 2., 2.], - [N, 3., 3., 4., 4., 4.5, 5., 5.], - [N, 3., 3., 4., 4., 4.5, 5., 5.], - [N, 4.5, 4.5, 5.5, 5.5, 6.0, 6.5, 6.5], - [N, 6., 6., 7., 7., 7.5, 8., 8.], - [N, 6., 6., 7., 7., 7.5, 8., 8.]]) + print("Passed subdivide with max parameter...") + print("Testing subdivide with mean parameter...") + finegrid = hostgrid.subdivide(finedict, cellFill="mean") + output = np.array( + [ + [N, 0.0, 0.0, 1.0, 1.0, 1.5, 2.0, 2.0], + [N, 0.0, 0.0, 1.0, 1.0, 1.5, 2.0, 2.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.5, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.5, 5.0, 5.0], + [N, 4.5, 4.5, 5.5, 5.5, 6.0, 6.5, 6.5], + [N, 6.0, 6.0, 7.0, 7.0, 7.5, 8.0, 8.0], + [N, 6.0, 6.0, 7.0, 7.0, 7.5, 8.0, 8.0], + ] + ) np.testing.assert_almost_equal(finegrid.getData(), output) - print('Passed subdivide with mean parameter...') - print('Passed subdivide method test - non-aligned grids.') + print("Passed subdivide with mean parameter...") + print("Passed subdivide method test - non-aligned grids.") def test_basics(): - geodict = GeoDict({'xmin': 0.5, 'xmax': 3.5, 'ymin': 0.5, - 'ymax': 3.5, 'dx': 1.0, 'dy': 1.0, 'ny': 4, 'nx': 4}) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 3.5, + "ymin": 0.5, + "ymax": 3.5, + "dx": 1.0, + "dy": 1.0, + "ny": 4, + "nx": 4, + } + ) data = np.arange(0, 16).reshape(4, 4).astype(np.float32) grid = Grid2D(data, geodict) - print('Testing basic Grid2D functionality (retrieving data, lat/lon to pixel coordinates, etc...') + print( + "Testing basic Grid2D functionality (retrieving data, lat/lon to pixel coordinates, etc..." + ) np.testing.assert_almost_equal(grid.getData(), data) assert grid.getGeoDict() == geodict @@ -150,70 +210,98 @@ def test_basics(): np.testing.assert_almost_equal(value, output) - print('Passed basic Grid2D functionality (retrieving data, lat/lon to pixel coordinates, etc...') + print( + "Passed basic Grid2D functionality (retrieving data, lat/lon to pixel coordinates, etc..." + ) def test_getvalue(): array = np.arange(1, 26).reshape(5, 5) - gdict = GeoDict({'xmin': 1.0, - 'xmax': 5.0, - 'ymin': 1.0, - 'ymax': 5.0, - 'dx': 1.0, - 'dy': 1.0, - 'nx': 5, - 'ny': 5}) + gdict = GeoDict( + { + "xmin": 1.0, + "xmax": 5.0, + "ymin": 1.0, + "ymax": 5.0, + "dx": 1.0, + "dy": 1.0, + "nx": 5, + "ny": 5, + } + ) grid = Grid2D(array, gdict) assert grid.getValue(3.0, 3.0) == 13 lat = np.array([3.0, 4.0]) lon = np.array([3.0, 3.0]) test = grid.getValue(lat, lon) np.testing.assert_almost_equal(test, np.array([13, 8])) - lat = np.array([[3.0, 4.0], - [4.0, 5.0]]) - lon = np.array([[3.0, 3.0], - [4.0, 4.0]]) + lat = np.array([[3.0, 4.0], [4.0, 5.0]]) + lon = np.array([[3.0, 3.0], [4.0, 4.0]]) test = grid.getValue(lat, lon) np.testing.assert_almost_equal(test, np.array([[13, 8], [9, 4]])) def test_cut(): - geodict = GeoDict({'xmin': 0.5, 'xmax': 4.5, 'ymin': 0.5, - 'ymax': 4.5, 'dx': 1.0, 'dy': 1.0, 'ny': 5, 'nx': 5}) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 4.5, + "ymin": 0.5, + "ymax": 4.5, + "dx": 1.0, + "dy": 1.0, + "ny": 5, + "nx": 5, + } + ) data = np.arange(0, 25).reshape(5, 5) - print('Testing data extraction...') + print("Testing data extraction...") grid = Grid2D(data, geodict) xmin, xmax, ymin, ymax = (2.5, 3.5, 2.5, 3.5) newgrid = grid.cut(xmin, xmax, ymin, ymax) output = np.array([[7, 8], [12, 13]]) np.testing.assert_almost_equal(newgrid.getData(), output) - print('Passed data extraction...') + print("Passed data extraction...") - print('Testing data trimming with resampling...') + print("Testing data trimming with resampling...") # make a more complicated test using getboundswithin data = np.arange(0, 84).reshape(7, 12) - geodict = GeoDict({'xmin': -180, 'xmax': 150, - 'ymin': -90, 'ymax': 90, - 'dx': 30, 'dy': 30, - 'nx': 12, 'ny': 7}) + geodict = GeoDict( + { + "xmin": -180, + "xmax": 150, + "ymin": -90, + "ymax": 90, + "dx": 30, + "dy": 30, + "nx": 12, + "ny": 7, + } + ) grid = Grid2D(data, geodict) - sampledict = GeoDict.createDictFromBox(-75, - 45, -45, 75, geodict.dx, geodict.dy) + sampledict = GeoDict.createDictFromBox(-75, 45, -45, 75, geodict.dx, geodict.dy) cutdict = geodict.getBoundsWithin(sampledict) newgrid = grid.cut(cutdict.xmin, cutdict.xmax, cutdict.ymin, cutdict.ymax) - output = np.array([[16, 17, 18, 19], - [28, 29, 30, 31], - [40, 41, 42, 43], - [52, 53, 54, 55]]) + output = np.array( + [[16, 17, 18, 19], [28, 29, 30, 31], [40, 41, 42, 43], [52, 53, 54, 55]] + ) np.testing.assert_almost_equal(newgrid.getData(), output) - print('Passed data trimming with resampling...') - - print('Test cut with self-alignment...') - geodict = GeoDict({'xmin': 0.5, 'xmax': 4.5, - 'ymin': 0.5, 'ymax': 6.5, - 'dx': 1.0, 'dy': 1.0, - 'nx': 5, 'ny': 7}) + print("Passed data trimming with resampling...") + + print("Test cut with self-alignment...") + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 4.5, + "ymin": 0.5, + "ymax": 6.5, + "dx": 1.0, + "dy": 1.0, + "nx": 5, + "ny": 7, + } + ) data = np.arange(0, 35).astype(np.float32).reshape(7, 5) grid = Grid2D(data, geodict) cutxmin = 1.7 @@ -221,34 +309,49 @@ def test_cut(): cutymin = 1.7 cutymax = 5.7 cutgrid = grid.cut(cutxmin, cutxmax, cutymin, cutymax, align=True) - output = np.array([[7, 8], - [12, 13], - [17, 18], - [22, 23]]) + output = np.array([[7, 8], [12, 13], [17, 18], [22, 23]]) np.testing.assert_almost_equal(cutgrid.getData(), output) - print('Passed cut with self-alignment.') + print("Passed cut with self-alignment.") def test_interpolate(): - geodict = GeoDict({'xmin': 0.5, 'xmax': 6.5, 'ymin': 1.5, - 'ymax': 6.5, 'dx': 1.0, 'dy': 1.0, 'ny': 6, 'nx': 7}) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 6.5, + "ymin": 1.5, + "ymax": 6.5, + "dx": 1.0, + "dy": 1.0, + "ny": 6, + "nx": 7, + } + ) data = np.arange(14, 56).reshape(6, 7) - for method in ['nearest', 'linear', 'cubic']: + for method in ["nearest", "linear", "cubic"]: print('Testing interpolate with method "%s"...' % method) grid = Grid2D(data, geodict) - sampledict = GeoDict({'xmin': 3.0, 'xmax': 4.0, - 'ymin': 3.0, 'ymax': 4.0, - 'dx': 1.0, 'dy': 1.0, - 'ny': 2, 'nx': 2}) + sampledict = GeoDict( + { + "xmin": 3.0, + "xmax": 4.0, + "ymin": 3.0, + "ymax": 4.0, + "dx": 1.0, + "dy": 1.0, + "ny": 2, + "nx": 2, + } + ) grid = grid.interpolateToGrid(sampledict, method=method) tgrid = grid.interpolate2(sampledict, method=method) - if method == 'nearest': - output = np.array([[30.0, 31.0], [37.0, 38.0]]) - elif method == 'linear': - output = np.array([[34., 35.], [41., 42.]]) - elif method == 'cubic': - output = np.array([[34., 35.], [41., 42.]]) + if method == "nearest": + output = np.array([[30.0, 32.0], [37.0, 39.0]]) + elif method == "linear": + output = np.array([[34.0, 35.0], [41.0, 42.0]]) + elif method == "cubic": + output = np.array([[34.0, 35.0], [41.0, 42.0]]) else: pass np.testing.assert_almost_equal(grid.getData(), output) @@ -262,13 +365,13 @@ def test_interpolate(): grid = Grid2D(data, geodict) sampledict = GeoDict.createDictFromBox(2, 8, 2, 8, 0.098, 0.098) t1 = time.time() - grid2 = grid.interpolateToGrid(sampledict, method='linear') + grid2 = grid.interpolateToGrid(sampledict, method="linear") t2 = time.time() - grid3 = grid.interpolate2(sampledict, method='linear') + grid3 = grid.interpolate2(sampledict, method="linear") t3 = time.time() # np.testing.assert_almost_equal(grid2._data.sum(),grid3._data.sum()) - print('scipy method: %.3f seconds' % (t2 - t1)) - print('gdal method: %.3f seconds' % (t3 - t2)) + print("scipy method: %.3f seconds" % (t2 - t1)) + print("gdal method: %.3f seconds" % (t3 - t2)) # test interpolate2 when called with geodict that is aligned with # enclosing geodict. This should just cut the grid. @@ -280,28 +383,34 @@ def test_interpolate(): nominal_lat_spacing = 0.0083 nlon = 283 nlat = 217 - host_geodict = GeoDict({'xmin': lon_min, - 'xmax': lon_max, - 'ymin': lat_min, - 'ymax': lat_max, - 'dx': nominal_lon_spacing, - 'dy': nominal_lat_spacing, - 'nx': nlon, - 'ny': nlat - }) + host_geodict = GeoDict( + { + "xmin": lon_min, + "xmax": lon_max, + "ymin": lat_min, + "ymax": lat_max, + "dx": nominal_lon_spacing, + "dy": nominal_lat_spacing, + "nx": nlon, + "ny": nlat, + } + ) sample_xmin = host_geodict.xmin + host_geodict.dx * 5 sample_xmax = host_geodict.xmax - host_geodict.dx * 5 sample_ymin = host_geodict.ymin + host_geodict.dy * 5 sample_ymax = host_geodict.ymax - host_geodict.dy * 5 - sample_geodict = GeoDict({'xmin': sample_xmin, - 'xmax': sample_xmax, - 'ymin': sample_ymin, - 'ymax': sample_ymax, - 'dx': host_geodict.dx, - 'dy': host_geodict.dy, - 'nx': nlon - 10, - 'ny': nlat - 10, - }) + sample_geodict = GeoDict( + { + "xmin": sample_xmin, + "xmax": sample_xmax, + "ymin": sample_ymin, + "ymax": sample_ymax, + "dx": host_geodict.dx, + "dy": host_geodict.dy, + "nx": nlon - 10, + "ny": nlat - 10, + } + ) assert host_geodict.isAligned(sample_geodict) host_data = np.random.rand(nlat, nlon) host_data = host_data.astype(np.float32) @@ -324,15 +433,19 @@ def test_interpolate(): ncols = int(((xmax - xmin) / dx) + 1) nrows = int(((ymax - ymin) / dy) + 1) # right/bottom edges of geodict will be adjusted if necessary - sample_geodict = GeoDict({'xmin': xmin, - 'xmax': xmax, - 'ymin': ymin, - 'ymax': ymax, - 'dx': dx, - 'dy': dy, - 'nx': ncols, - 'ny': nrows}, adjust='bounds' - ) + sample_geodict = GeoDict( + { + "xmin": xmin, + "xmax": xmax, + "ymin": ymin, + "ymax": ymax, + "dx": dx, + "dy": dy, + "nx": ncols, + "ny": nrows, + }, + adjust="bounds", + ) assert not host_geodict.isAligned(sample_geodict) host_data = np.random.randint(0, 100, size=(nlat, nlon), dtype=np.int16) host_grid = Grid2D(data=host_data, geodict=host_geodict) @@ -341,43 +454,65 @@ def test_interpolate(): def test_rasterize(): - geodict = GeoDict({'xmin': 0.5, 'xmax': 3.5, - 'ymin': 0.5, 'ymax': 3.5, - 'dx': 1.0, 'dy': 1.0, - 'ny': 4, 'nx': 4}) - print('Testing rasterizeFromGeometry() burning in values from a polygon sequence...') + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 3.5, + "ymin": 0.5, + "ymax": 3.5, + "dx": 1.0, + "dy": 1.0, + "ny": 4, + "nx": 4, + } + ) + print( + "Testing rasterizeFromGeometry() burning in values from a polygon sequence..." + ) # Define two simple polygons and assign them to shapes poly1 = [(0.25, 3.75), (1.25, 3.25), (1.25, 2.25)] - poly2 = [(2.25, 3.75), (3.25, 3.75), (3.75, 2.75), - (3.75, 1.50), (3.25, 0.75), (2.25, 2.25)] - shape1 = {'properties': {'value': 5}, 'geometry': mapping(Polygon(poly1))} - shape2 = {'properties': {'value': 7}, 'geometry': mapping(Polygon(poly2))} + poly2 = [ + (2.25, 3.75), + (3.25, 3.75), + (3.75, 2.75), + (3.75, 1.50), + (3.25, 0.75), + (2.25, 2.25), + ] + shape1 = {"properties": {"value": 5}, "geometry": mapping(Polygon(poly1))} + shape2 = {"properties": {"value": 7}, "geometry": mapping(Polygon(poly2))} shapes = [shape1, shape2] - print('Testing burning in values where polygons need not contain pixel centers...') + print("Testing burning in values where polygons need not contain pixel centers...") grid = Grid2D.rasterizeFromGeometry( - shapes, geodict, fillValue=0, attribute='value', mustContainCenter=False) - output = np.array([[5, 5, 7, 7], - [5, 5, 7, 7], - [0, 0, 7, 7], - [0, 0, 0, 7]]) + shapes, geodict, fillValue=0, attribute="value", mustContainCenter=False + ) + output = np.array([[5, 5, 7, 7], [5, 5, 7, 7], [0, 0, 7, 7], [0, 0, 0, 7]]) np.testing.assert_almost_equal(grid.getData(), output) - print('Passed burning in values where polygons need not contain pixel centers.') + print("Passed burning in values where polygons need not contain pixel centers.") - print('Testing burning in values where polygons must contain pixel centers...') + print("Testing burning in values where polygons must contain pixel centers...") grid2 = Grid2D.rasterizeFromGeometry( - shapes, geodict, fillValue=0, attribute='value', mustContainCenter=True) - output = np.array([[5, 0, 7, 0], - [0, 0, 7, 7], - [0, 0, 0, 7], - [0, 0, 0, 0]]) + shapes, geodict, fillValue=0, attribute="value", mustContainCenter=True + ) + output = np.array([[5, 0, 7, 0], [0, 0, 7, 7], [0, 0, 0, 7], [0, 0, 0, 0]]) np.testing.assert_almost_equal(grid2.getData(), output) - print('Passed burning in values where polygons must contain pixel centers.') + print("Passed burning in values where polygons must contain pixel centers.") def test_copy(): data = np.arange(0, 16).astype(np.float32).reshape(4, 4) - geodict = GeoDict({'xmin': 0.5, 'xmax': 3.5, 'ymin': 0.5, - 'ymax': 3.5, 'dx': 1.0, 'dy': 1.0, 'ny': 4, 'nx': 4}) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 3.5, + "ymin": 0.5, + "ymax": 3.5, + "dx": 1.0, + "dy": 1.0, + "ny": 4, + "nx": 4, + } + ) grid1 = Grid2D(data, geodict) grid2 = grid1.copyFromGrid(grid1) grid1._data[0, 0] = np.nan @@ -387,76 +522,110 @@ def test_copy(): def test_setData(): data = np.arange(0, 16).astype(np.float32).reshape(4, 4) - geodict = GeoDict({'xmin': 0.5, 'xmax': 3.5, 'ymin': 0.5, - 'ymax': 3.5, 'dx': 1.0, 'dy': 1.0, 'ny': 4, 'nx': 4}) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 3.5, + "ymin": 0.5, + "ymax": 3.5, + "dx": 1.0, + "dy": 1.0, + "ny": 4, + "nx": 4, + } + ) grid1 = Grid2D(data, geodict) x = np.ones((4, 4)) try: grid1.setData(x) # this should pass - print('setData test passed.') + print("setData test passed.") except DataSetException as dse: - print('setData test failed.') + print("setData test failed.") try: x = np.ones((5, 5)) grid1.setData(x) - print('setData test did not fail when it should have.') + print("setData test did not fail when it should have.") except DataSetException as dse: - print('setData test failed as expected.') + print("setData test failed as expected.") try: - x = 'fred' + x = "fred" grid1.setData(x) - print('setData test did not fail when it should have.') + print("setData test did not fail when it should have.") except DataSetException as dse: - print('setData test failed as expected.') + print("setData test failed as expected.") def get_data_range_test(): # a standard global grid, going from -180 to 180 - normal_dict = GeoDict({'xmin': -180, 'xmax': 120, - 'ymin': -90, 'ymax': 90, - 'dx': 60, 'dy': 45, - 'nx': 6, 'ny': 5}) + normal_dict = GeoDict( + { + "xmin": -180, + "xmax": 120, + "ymin": -90, + "ymax": 90, + "dx": 60, + "dy": 45, + "nx": 6, + "ny": 5, + } + ) # test a simple example which does NOT cross the 180 meridian sample1 = (-125, 65, -20, 20) dict1 = Grid2D.getDataRange(normal_dict, sample1) - cdict1 = {'iulx1': 0, 'iuly1': 1, - 'ilrx1': 6, 'ilry1': 4} + cdict1 = {"iulx1": 0, "iuly1": 1, "ilrx1": 6, "ilry1": 4} assert dict1 == cdict1 # test a less-simple example which DOES cross the 180 meridian sample2 = (-235, -10, -20, 20) dict2 = Grid2D.getDataRange(normal_dict, sample2) - cdict2 = {'iulx1': 5, 'iuly1': 1, - 'ilrx1': 6, 'ilry1': 4, - 'iulx2': 0, 'iuly2': 1, - 'ilrx2': 4, 'ilry2': 4} + cdict2 = { + "iulx1": 5, + "iuly1": 1, + "ilrx1": 6, + "ilry1": 4, + "iulx2": 0, + "iuly2": 1, + "ilrx2": 4, + "ilry2": 4, + } assert dict2 == cdict2 # test a less-simple example which DOES cross the 180 meridian, and xmin > xmax sample3 = (125, -10, -20, 20) dict3 = Grid2D.getDataRange(normal_dict, sample3) - cdict3 = {'iulx1': 5, 'iuly1': 1, - 'ilrx1': 6, 'ilry1': 4, - 'iulx2': 0, 'iuly2': 1, - 'ilrx2': 4, 'ilry2': 4} + cdict3 = { + "iulx1": 5, + "iuly1": 1, + "ilrx1": 6, + "ilry1": 4, + "iulx2": 0, + "iuly2": 1, + "ilrx2": 4, + "ilry2": 4, + } assert dict3 == cdict3 # test an example where the sample bounds are from 0 to 360 sample4 = (160, 200, -20, 20) dict4 = Grid2D.getDataRange(normal_dict, sample4) - cdict4 = {'iulx1': 5, 'iuly1': 1, - 'ilrx1': 6, 'ilry1': 4, - 'iulx2': 0, 'iuly2': 1, - 'ilrx2': 2, 'ilry2': 4} + cdict4 = { + "iulx1": 5, + "iuly1": 1, + "ilrx1": 6, + "ilry1": 4, + "iulx2": 0, + "iuly2": 1, + "ilrx2": 2, + "ilry2": 4, + } assert dict4 == cdict4 # test an example where the sample bounds are from 0 to 360 sample5 = (220, 260, -20, 20) dict5 = Grid2D.getDataRange(normal_dict, sample5) - cdict5 = {'iulx1': 0, 'iuly1': 1, - 'ilrx1': 3, 'ilry1': 4} + cdict5 = {"iulx1": 0, "iuly1": 1, "ilrx1": 3, "ilry1": 4} assert dict5 == cdict5 @@ -467,34 +636,45 @@ def test_project(): data = np.arange(0.0, ncells).reshape(gd.ny, gd.nx) grid = GDALGrid(data, gd) projstr = "+proj=merc +lat_ts=55 +lon_0=180 +ellps=WGS84" - newgrid = grid.project(projstr, method='nearest') + newgrid = grid.project(projstr, method="nearest") proj = pyproj.Proj(projstr) # what would the ul/lr corners be? ulx, uly = proj(grid._geodict.xmin, grid._geodict.ymax) lrx, lry = proj(grid._geodict.xmax, grid._geodict.ymin) # what if we back-project? - newxmin, newymax = proj(newgrid._geodict.xmin, - newgrid._geodict.ymax, inverse=True) - newxmax, newymin = proj(newgrid._geodict.xmax, - newgrid._geodict.ymin, inverse=True) + newxmin, newymax = proj(newgrid._geodict.xmin, newgrid._geodict.ymax, inverse=True) + newxmax, newymin = proj(newgrid._geodict.xmax, newgrid._geodict.ymin, inverse=True) x = 1 # test simple projection - data = np.array([[0, 0, 1, 0, 0], - [0, 0, 1, 0, 0], - [1, 1, 1, 1, 1], - [0, 0, 1, 0, 0], - [0, 0, 1, 0, 0]], dtype=np.int32) - geodict = {'xmin': 50, 'xmax': 50.4, 'ymin': 50, - 'ymax': 50.4, 'dx': 0.1, 'dy': 0.1, 'nx': 5, 'ny': 5} + data = np.array( + [ + [0, 0, 1, 0, 0], + [0, 0, 1, 0, 0], + [1, 1, 1, 1, 1], + [0, 0, 1, 0, 0], + [0, 0, 1, 0, 0], + ], + dtype=np.int32, + ) + geodict = { + "xmin": 50, + "xmax": 50.4, + "ymin": 50, + "ymax": 50.4, + "dx": 0.1, + "dy": 0.1, + "nx": 5, + "ny": 5, + } gd = GeoDict(geodict) grid = GDALGrid(data, gd) projstr = "+proj=utm +zone=40 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs " - newgrid = grid.project(projstr, method='nearest') + newgrid = grid.project(projstr, method="nearest") try: tdir = tempfile.mkdtemp() - outfile = os.path.join(tdir, 'output.bil') + outfile = os.path.join(tdir, "output.bil") grid.save(outfile) with rasterio.open(outfile) as src: aff = get_affine(src) @@ -505,20 +685,21 @@ def test_project(): left = aff.xoff top = aff.yoff right, bottom = aff * (ncols - 1, nrows - 1) - dst_transform, width, height = calculate_default_transform(src_crs, dst_crs, - ncols, nrows, - left, bottom, - right, top) + dst_transform, width, height = calculate_default_transform( + src_crs, dst_crs, ncols, nrows, left, bottom, right, top + ) destination = np.zeros((height, width)) - reproject(data, - destination, - src_transform=aff, - src_crs=src_crs, - dst_transform=dst_transform, - dst_crs=dst_crs, - src_nodata=src.nodata, - dst_nodata=np.nan, - resampling=Resampling.nearest) + reproject( + data, + destination, + src_transform=aff, + src_crs=src_crs, + dst_transform=dst_transform, + dst_crs=dst_crs, + src_nodata=src.nodata, + dst_nodata=np.nan, + resampling=Resampling.nearest, + ) x = 1 except: pass @@ -546,22 +727,26 @@ def test_project(): def test_fiji(): - host_dict = {'xmin': 176.8874999998576, - 'xmax': -178.23750000014437, - 'ymin': -20.770833333331773, - 'ymax': -16.154166666666953, - 'dx': 0.00833333333333, - 'dy': 0.00833333333333, - 'ny': 555, - 'nx': 586} - sample_dict = {'xmin': 176.90416666666522, - 'xmax': -178.25416666666814, - 'ymin': -20.729166666666934, - 'ymax': -16.154166666666953, - 'dx': 0.0083333333333333, - 'dy': 0.0083333333333333, - 'ny': 550, - 'nx': 582} + host_dict = { + "xmin": 176.8874999998576, + "xmax": -178.23750000014437, + "ymin": -20.770833333331773, + "ymax": -16.154166666666953, + "dx": 0.00833333333333, + "dy": 0.00833333333333, + "ny": 555, + "nx": 586, + } + sample_dict = { + "xmin": 176.90416666666522, + "xmax": -178.25416666666814, + "ymin": -20.729166666666934, + "ymax": -16.154166666666953, + "dx": 0.0083333333333333, + "dy": 0.0083333333333333, + "ny": 550, + "nx": 582, + } host_geodict = GeoDict(host_dict) sample_geodict = GeoDict(sample_dict) host_data = np.zeros((host_geodict.ny, host_geodict.nx)) @@ -574,20 +759,26 @@ def test_fiji(): # this is a case where the host grid crosses # the meridian but the sample grid is only on # the east side of the meridian. - host_dict = {'xmin': 176.11235968666102, - 'xmax': -179.99597366223898, - 'ymin': -21.212639164039008, - 'ymax': -17.537639178739006, - 'dx': 0.008333333299999992, - 'dy': 0.008333333300000002, - 'ny': 442, 'nx': 468} - sample_dict = {'xmin': 176.15402635316102, - 'xmax': 179.94569300466102, - 'ymin': -21.162639164239007, - 'ymax': -17.587639178539007, - 'dx': 0.008333333299999992, - 'dy': 0.008333333300000002, - 'ny': 430, 'nx': 456} + host_dict = { + "xmin": 176.11235968666102, + "xmax": -179.99597366223898, + "ymin": -21.212639164039008, + "ymax": -17.537639178739006, + "dx": 0.008333333299999992, + "dy": 0.008333333300000002, + "ny": 442, + "nx": 468, + } + sample_dict = { + "xmin": 176.15402635316102, + "xmax": 179.94569300466102, + "ymin": -21.162639164239007, + "ymax": -17.587639178539007, + "dx": 0.008333333299999992, + "dy": 0.008333333300000002, + "ny": 430, + "nx": 456, + } host_geodict = GeoDict(host_dict) sample_geodict = GeoDict(sample_dict) host_data = np.zeros((host_geodict.ny, host_geodict.nx)) @@ -598,7 +789,7 @@ def test_fiji(): np.testing.assert_almost_equal(yi[0], ycmp) -if __name__ == '__main__': +if __name__ == "__main__": test_fiji() test_getvalue() test_project() diff --git a/test/multiple_test.py b/test/multiple_test.py index 4acc782..d185fa2 100755 --- a/test/multiple_test.py +++ b/test/multiple_test.py @@ -1,24 +1,26 @@ #!/usr/bin/env python -#python 3 compatibility +# python 3 compatibility from __future__ import print_function import os.path import sys -#stdlib imports +# stdlib imports import abc import textwrap import glob import os from collections import OrderedDict -#hack the path so that I can debug these functions if I need to -homedir = os.path.dirname(os.path.abspath(__file__)) #where is this script? -mapiodir = os.path.abspath(os.path.join(homedir,'..')) -sys.path.insert(0,mapiodir) #put this at the front of the system path, ignoring any installed mapio stuff +# hack the path so that I can debug these functions if I need to +homedir = os.path.dirname(os.path.abspath(__file__)) # where is this script? +mapiodir = os.path.abspath(os.path.join(homedir, "..")) +sys.path.insert( + 0, mapiodir +) # put this at the front of the system path, ignoring any installed mapio stuff -#third party imports +# third party imports from mapio.gridbase import Grid from mapio.multiple import MultiGrid from mapio.grid2d import Grid2D @@ -29,63 +31,103 @@ import shapely from affine import Affine from rasterio import features -from shapely.geometry import MultiPoint,Polygon,mapping +from shapely.geometry import MultiPoint, Polygon, mapping + def test(): - print('Testing MultiGrid interpolate...') - data = np.arange(14,56).reshape(6,7) - geodict = GeoDict({'xmin':0.5,'xmax':6.5,'ymin':1.5,'ymax':6.5,'dx':1.0,'dy':1.0,'ny':6,'nx':7}) + print("Testing MultiGrid interpolate...") + data = np.arange(14, 56).reshape(6, 7) + geodict = GeoDict( + { + "xmin": 0.5, + "xmax": 6.5, + "ymin": 1.5, + "ymax": 6.5, + "dx": 1.0, + "dy": 1.0, + "ny": 6, + "nx": 7, + } + ) layers = OrderedDict() - layers['layer1'] = Grid2D(data,geodict) + layers["layer1"] = Grid2D(data, geodict) mgrid = MultiGrid(layers) - sampledict = GeoDict({'xmin':3.0,'xmax':4.0, - 'ymin':3.0,'ymax':4.0, - 'dx':1.0,'dy':1.0, - 'ny':2,'nx':2}) - for method in ['nearest','linear','cubic']: - mgrid2 = mgrid.interpolateToGrid(sampledict,method=method) - if method == 'nearest': - output = np.array([[30.0,31.0],[37.0,38.0]]) - elif method == 'linear': - output = np.array([[34.,35.],[41.,42.]]) - elif method == 'cubic': - output = np.array([[34.,35.],[41.,42.]]) + sampledict = GeoDict( + { + "xmin": 3.0, + "xmax": 4.0, + "ymin": 3.0, + "ymax": 4.0, + "dx": 1.0, + "dy": 1.0, + "ny": 2, + "nx": 2, + } + ) + for method in ["nearest", "linear", "cubic"]: + mgrid2 = mgrid.interpolateToGrid(sampledict, method=method) + if method == "nearest": + output = np.array([[30.0, 32.0], [37.0, 39.0]]) + elif method == "linear": + output = np.array([[34.0, 35.0], [41.0, 42.0]]) + elif method == "cubic": + output = np.array([[34.0, 35.0], [41.0, 42.0]]) else: pass - np.testing.assert_almost_equal(mgrid2.getLayer('layer1').getData(),output) - print('Passed MultiGrid interpolate test.') + np.testing.assert_almost_equal(mgrid2.getLayer("layer1").getData(), output) + print("Passed MultiGrid interpolate test.") - print('Testing bounds retrieval...') + print("Testing bounds retrieval...") b1 = np.array(mgrid.getBounds()) - b2 = np.array((geodict.xmin,geodict.xmax,geodict.ymin,geodict.ymax)) - np.testing.assert_almost_equal(b1,b2) - print('Passed bounds retrieval...') - - print('Testing MultiGrid subdivide test...') - data = np.arange(0,9).reshape((3,3)) - geodict = GeoDict({'xmin':0.0,'xmax':10.0, - 'ymin':0.0,'ymax':10.0, - 'dx':5.0,'dy':5.0, - 'ny':3,'nx':3}) + b2 = np.array((geodict.xmin, geodict.xmax, geodict.ymin, geodict.ymax)) + np.testing.assert_almost_equal(b1, b2) + print("Passed bounds retrieval...") + + print("Testing MultiGrid subdivide test...") + data = np.arange(0, 9).reshape((3, 3)) + geodict = GeoDict( + { + "xmin": 0.0, + "xmax": 10.0, + "ymin": 0.0, + "ymax": 10.0, + "dx": 5.0, + "dy": 5.0, + "ny": 3, + "nx": 3, + } + ) layers = OrderedDict() - layers['layer1'] = Grid2D(data,geodict) + layers["layer1"] = Grid2D(data, geodict) hostgrid = MultiGrid(layers) - finedict = GeoDict({'xmin':-2.5,'xmax':11.5, - 'ymin':-1.5,'ymax':10.5, - 'dx':2.0,'dy':2.0, - 'nx':8,'ny':7}) + finedict = GeoDict( + { + "xmin": -2.5, + "xmax": 11.5, + "ymin": -1.5, + "ymax": 10.5, + "dx": 2.0, + "dy": 2.0, + "nx": 8, + "ny": 7, + } + ) N = np.nan - finegrid = hostgrid.subdivide(finedict,cellFill='min') - output = np.array([[ N, 0., 0., 1., 1., 1., 2., 2.], - [ N, 0., 0., 1., 1., 1., 2., 2.], - [ N, 3., 3., 4., 4., 4., 5., 5.], - [ N, 3., 3., 4., 4., 4., 5., 5.], - [ N, 3., 3., 4., 4., 4., 5., 5.], - [ N, 6., 6., 7., 7., 7., 8., 8.], - [ N, 6., 6., 7., 7., 7., 8., 8.]]) - np.testing.assert_almost_equal(finegrid.getLayer('layer1').getData(),output) - print('Passed MultiGrid subdivide test.') - + finegrid = hostgrid.subdivide(finedict, cellFill="min") + output = np.array( + [ + [N, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0], + [N, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 3.0, 3.0, 4.0, 4.0, 4.0, 5.0, 5.0], + [N, 6.0, 6.0, 7.0, 7.0, 7.0, 8.0, 8.0], + [N, 6.0, 6.0, 7.0, 7.0, 7.0, 8.0, 8.0], + ] + ) + np.testing.assert_almost_equal(finegrid.getLayer("layer1").getData(), output) + print("Passed MultiGrid subdivide test.") + -if __name__ == '__main__': +if __name__ == "__main__": test() From e44b80642183af8e276083c8a2dc7127fe36c4ce Mon Sep 17 00:00:00 2001 From: Michael Hearne Date: Mon, 8 Nov 2021 10:23:57 -0700 Subject: [PATCH 3/3] Added new attributes in cdf/hdf output to be compliant with newer gdal/rasterio standards --- mapio/writer.py | 172 +++++++++++++++++++++++--------------------- test/writer_test.py | 109 +++++++++++++++------------- 2 files changed, 149 insertions(+), 132 deletions(-) diff --git a/mapio/writer.py b/mapio/writer.py index f20c595..a1c1b09 100644 --- a/mapio/writer.py +++ b/mapio/writer.py @@ -29,7 +29,7 @@ def get_nodata(grid): zmin = np.nanmin(grid._data) zmax = np.nanmax(grid._data) if grid._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: @@ -37,7 +37,7 @@ def get_nodata(grid): # 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: @@ -49,51 +49,63 @@ def get_nodata(grid): def _getHeader(grid): 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'] = grid._data.shape - hdr['NBANDS'] = 1 + hdr["BYTEORDER"] = "MSBFIRST" + hdr["LAYOUT"] = "BIL" + hdr["NROWS"], hdr["NCOLS"] = grid._data.shape + hdr["NBANDS"] = 1 if grid._data.dtype == np.uint8: - hdr['NBITS'] = 8 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 8 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif grid._data.dtype == np.int8: - hdr['NBITS'] = 8 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 8 + hdr["PIXELTYPE"] = "SIGNEDINT" elif grid._data.dtype == np.uint16: - hdr['NBITS'] = 16 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 16 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif grid._data.dtype == np.int16: - hdr['NBITS'] = 16 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 16 + hdr["PIXELTYPE"] = "SIGNEDINT" elif grid._data.dtype == np.uint32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'UNSIGNEDINT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "UNSIGNEDINT" elif grid._data.dtype == np.int32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'SIGNEDINT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "SIGNEDINT" elif grid._data.dtype == np.float32: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'FLOAT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "FLOAT" elif grid._data.dtype == np.float64: - hdr['NBITS'] = 32 - hdr['PIXELTYPE'] = 'FLOAT' + hdr["NBITS"] = 32 + hdr["PIXELTYPE"] = "FLOAT" else: - raise KeyError( - 'Data type "%s" not supported.' % str(grid._data.dtype)) - hdr['BANDROWBYTES'] = hdr['NCOLS'] * (hdr['NBITS'] / 8) - hdr['TOTALROWBYTES'] = hdr['NCOLS'] * (hdr['NBITS'] / 8) - hdr['ULXMAP'] = grid._geodict.xmin - hdr['ULYMAP'] = grid._geodict.ymax - hdr['XDIM'] = grid._geodict.dx - hdr['YDIM'] = grid._geodict.dy + raise KeyError('Data type "%s" not supported.' % str(grid._data.dtype)) + hdr["BANDROWBYTES"] = hdr["NCOLS"] * (hdr["NBITS"] / 8) + hdr["TOTALROWBYTES"] = hdr["NCOLS"] * (hdr["NBITS"] / 8) + hdr["ULXMAP"] = grid._geodict.xmin + hdr["ULYMAP"] = grid._geodict.ymax + hdr["XDIM"] = grid._geodict.dx + hdr["YDIM"] = grid._geodict.dy NODATA = get_nodata(grid) - 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] @@ -104,46 +116,42 @@ def write(grid, filename, format_type, do_compression=False): """Save a GMTGrid object to a file. Args: filename (str): Name of desired output file. - format_type (str): One of 'netcdf','hdf' or 'esri'. + format_type (str): One of 'hdf', 'esri', 'tiff'. Raises: - KeyError -- When format_type not one of ('netcdf,'hdf','native') + KeyError -- When format_type not one of ('netcdf', 'hdf', 'esri', 'tiff') """ - if format_type not in ['netcdf', 'hdf', 'esri', 'tiff']: - msg = ('Only "netcdf", "hdf", "esri", and ' - '"tiff" formats are supported.') + if format_type not in ["netcdf", "hdf", "esri", "tiff"]: + msg = 'Only "netcdf", "hdf", "esri", and "tiff" formats are supported.' raise KeyError(msg) - if format_type in ['netcdf', 'hdf']: - format_mapper = {'hdf': 'NETCDF4', - 'netcdf': 'NETCDF3_CLASSIC'} + if format_type in ["netcdf", "hdf"]: + format_mapper = {"hdf": "NETCDF4", "netcdf": "NETCDF3_CLASSIC"} - f = netcdf.Dataset(filename, 'w', format=format_mapper[format_type]) + f = netcdf.Dataset(filename, "w", format=format_mapper[format_type]) NODATA = get_nodata(grid) if NODATA is not None: f.setncattr("_Fill_Value", NODATA) m, n = grid._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')) + dx = f.createDimension("x", n) # noqa + dy = f.createDimension("y", m) # noqa + x = f.createVariable("x", np.float64, ("x")) + x.axis = 'X' + y = f.createVariable("y", np.float64, ("y")) + y.axis = 'Y' if grid._geodict.xmin < grid._geodict.xmax: - x[:] = np.linspace(grid._geodict.xmin, - grid._geodict.xmax, grid._geodict.nx) + x[:] = np.linspace(grid._geodict.xmin, grid._geodict.xmax, grid._geodict.nx) else: - x[:] = np.linspace(grid._geodict.xmin, - grid._geodict.xmax + 360, grid._geodict.nx) - y[:] = np.linspace(grid._geodict.ymin, - grid._geodict.ymax, grid._geodict.ny) + x[:] = np.linspace( + grid._geodict.xmin, grid._geodict.xmax + 360, grid._geodict.nx + ) + y[:] = np.linspace(grid._geodict.ymin, grid._geodict.ymax, grid._geodict.ny) if NODATA is not None: - z = f.createVariable('z', grid._data.dtype, - ('y', 'x'), fill_value=NODATA) + z = f.createVariable("z", grid._data.dtype, ("y", "x"), fill_value=NODATA) else: - z = f.createVariable('z', grid._data.dtype, - ('y', 'x')) + z = f.createVariable("z", grid._data.dtype, ("y", "x")) z[:] = np.flipud(grid._data) - z.actual_range = np.array( - (np.nanmin(grid._data), np.nanmax(grid._data))) + z.actual_range = np.array((np.nanmin(grid._data), np.nanmax(grid._data))) f.close() - elif format_type == 'esri': + elif format_type == "esri": hdr = _getHeader(grid) # create a reference to the data - this may be overridden by a # downcasted version for doubles @@ -151,41 +159,41 @@ def write(grid, filename, format_type, do_compression=False): if grid._data.dtype == np.float32: # so we can find/reset nan values without screwing up original data data = grid._data.astype(np.float32) - data[np.isnan(data)] = hdr['NODATA'] + data[np.isnan(data)] = hdr["NODATA"] elif grid._data.dtype == np.float64: data = grid._data.astype(np.float32) - data[np.isnan(data)] = hdr['NODATA'] - warnings.warn(UserWarning('Down-casting double precision ' - 'floating point to single precision')) + data[np.isnan(data)] = hdr["NODATA"] + warnings.warn( + UserWarning( + "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() - elif format_type == 'tiff': + elif format_type == "tiff": geodict = grid.getGeoDict() transform = affine_from_geodict(geodict) profile = DefaultGTiffProfile() # default profile has compression on - turn it off if requested. if not do_compression: - profile.pop('compress', None) + profile.pop("compress", None) # stufff - profile['height'] = geodict.ny - profile['width'] = geodict.nx - profile['count'] = 1 - profile['dtype'] = grid._data.dtype - profile['crs'] = geodict.projection - profile['transform'] = transform - new_dataset = rasterio.open( - filename, - 'w', - **profile - ) + profile["height"] = geodict.ny + profile["width"] = geodict.nx + profile["count"] = 1 + profile["dtype"] = grid._data.dtype + profile["crs"] = geodict.projection + profile["transform"] = transform + new_dataset = rasterio.open(filename, "w", **profile) new_dataset.write(grid._data, 1) new_dataset.close() diff --git a/test/writer_test.py b/test/writer_test.py index 5b705b2..93c2132 100755 --- a/test/writer_test.py +++ b/test/writer_test.py @@ -16,23 +16,25 @@ def test_write(): data = np.arange(0, 25).reshape(5, 5).astype(np.float32) - gdict = {'xmin': 5.0, - 'xmax': 9.0, - 'ymin': 4.0, - 'ymax': 8.0, - 'dx': 1.0, - 'dy': 1.0, - 'nx': 5, - 'ny': 5} + gdict = { + "xmin": 5.0, + "xmax": 9.0, + "ymin": 4.0, + "ymax": 8.0, + "dx": 1.0, + "dy": 1.0, + "nx": 5, + "ny": 5, + } gd = GeoDict(gdict) grid = Grid2D(data, gd) - for format_type in ['netcdf', 'esri', 'hdf', 'tiff']: + for format_type in ["netcdf", "esri", "hdf", "tiff"]: tdir = tempfile.mkdtemp() - fname = os.path.join(tdir, 'tempfile.grd') + fname = os.path.join(tdir, "tempfile.grd") try: write(grid, fname, format_type) - src = rasterio.open(fname, 'r') + src = rasterio.open(fname, "r") tdata = src.read(1) np.testing.assert_almost_equal(tdata, data) except Exception as e: @@ -44,26 +46,28 @@ def test_write(): def test_compress(): data = np.arange(0, 10000).reshape(100, 100).astype(np.float32) nx, ny = data.shape - gdict = {'xmin': 5.0, - 'xmax': 9.0, - 'ymin': 4.0, - 'ymax': 8.0, - 'dx': 1.0, - 'dy': 1.0, - 'nx': nx, - 'ny': ny} + gdict = { + "xmin": 5.0, + "xmax": 9.0, + "ymin": 4.0, + "ymax": 8.0, + "dx": 1.0, + "dy": 1.0, + "nx": nx, + "ny": ny, + } gd = GeoDict(gdict) grid = Grid2D(data, gd) tdir = tempfile.mkdtemp() - fname = os.path.join(tdir, 'tempfile.tif') + fname = os.path.join(tdir, "tempfile.tif") try: - write(grid, fname, 'tiff', do_compression=False) + write(grid, fname, "tiff", do_compression=False) uncompressed_size = pathlib.Path(fname).stat().st_size - write(grid, fname, 'tiff', do_compression=True) + write(grid, fname, "tiff", do_compression=True) compressed_size = pathlib.Path(fname).stat().st_size assert compressed_size < uncompressed_size - src = rasterio.open(fname, 'r') + src = rasterio.open(fname, "r") tdata = src.read(1) np.testing.assert_almost_equal(tdata, data) except Exception as e: @@ -75,22 +79,24 @@ def test_compress(): def test_nan_write(): data = np.arange(0, 25).reshape(5, 5).astype(np.float32) data[0, 0] = np.nan - gdict = {'xmin': 5.0, - 'xmax': 9.0, - 'ymin': 4.0, - 'ymax': 8.0, - 'dx': 1.0, - 'dy': 1.0, - 'nx': 5, - 'ny': 5} + gdict = { + "xmin": 5.0, + "xmax": 9.0, + "ymin": 4.0, + "ymax": 8.0, + "dx": 1.0, + "dy": 1.0, + "nx": 5, + "ny": 5, + } gd = GeoDict(gdict) grid = Grid2D(data, gd) - for format_type in ['netcdf', 'esri', 'hdf', 'tiff']: + for format_type in ["netcdf", "esri", "hdf", "tiff"]: tdir = tempfile.mkdtemp() - fname = os.path.join(tdir, 'tempfile.grd') + fname = os.path.join(tdir, "tempfile.grd") try: write(grid, fname, format_type) - src = rasterio.open(fname, 'r') + src = rasterio.open(fname, "r") tdata = src.read(1) np.testing.assert_almost_equal(tdata, data) except Exception as e: @@ -109,10 +115,10 @@ def big_test(): gd = GeoDict.createDictFromBox(xmin, xmax, ymin, ymax, dx, dy) data = np.random.rand(gd.ny, gd.nx) grid = Grid2D(data, gd) - fname = os.path.join(os.path.expanduser('~'), 'tempfile.grd') - write(grid, fname, 'hdf') + fname = os.path.join(os.path.expanduser("~"), "tempfile.grd") + write(grid, fname, "hdf") print(fname) - src = rasterio.open(fname, 'r') + src = rasterio.open(fname, "r") # tdata = src.read(1) # np.testing.assert_almost_equal(tdata, data) @@ -126,24 +132,27 @@ def test_meridian(): dy = 1 nx = int((((360 + xmax) - xmin) / dx) + 1) ny = int(((ymax - ymin) / dy) + 1) - geodict = GeoDict({'xmin': xmin, - 'xmax': xmax, - 'ymin': ymin, - 'ymax': ymax, - 'dx': dx, - 'dy': dy, - 'nx': nx, - 'ny': ny - }) + geodict = GeoDict( + { + "xmin": xmin, + "xmax": xmax, + "ymin": ymin, + "ymax": ymax, + "dx": dx, + "dy": dy, + "nx": nx, + "ny": ny, + } + ) data = np.arange(0, nx * ny, dtype=np.int32).reshape(ny, nx) grid = Grid2D(data=data, geodict=geodict) tdir = tempfile.mkdtemp() - fname = os.path.join(tdir, 'tempfile.grd') + fname = os.path.join(tdir, "tempfile.grd") try: - write(grid, fname, 'netcdf') + write(grid, fname, "netcdf") grid2 = read(fname) assert grid2._geodict == geodict - write(grid, fname, 'tiff') + write(grid, fname, "tiff") grid3 = read(fname) assert grid3._geodict == geodict except Exception as e: @@ -152,7 +161,7 @@ def test_meridian(): shutil.rmtree(tdir) -if __name__ == '__main__': +if __name__ == "__main__": # big_test() test_meridian() test_compress()