diff --git a/README.md b/README.md index ce2e48a..26f93d7 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,15 @@ pip install sardem ``` This creates the command line executable `sardem` +Using `uv`: +```bash +uv add sardem +``` +You can also run the `sardem` command line tool using `uvx`: +```bash +uvx sardem --help +``` + Alternatively, you can clone to build/install: ```bash @@ -32,41 +41,23 @@ cd sardem # Install requirements using either pip or conda # conda install -c conda-forge --file environment.yml # pip install -r requirements.txt -# the conda environment.yml is more complete, as GDAL is required for some of the functionality +# the conda environment.yml is more complete, as rasterio is required for some of the functionality pip install -e . ``` which will run `pip install --upgrade .` and create the command line script. - ## Data sources The default data source, `--data-source NASA`, uses the SRTM 1 arcsecond data. You can also use the newer [Copernicus Digital Surface Model (DSM)](https://registry.opendata.aws/copernicus-dem/). To see a comparison of the two, see the [srtm_copernicus_comparison](notebooks/srtm_copernicus_comparison.ipynb) notebook. -**Note:** To convert the elevation values to heights about the WGS84 ellipsoid (which is the default), or to use the Copernicus data, **GDAL is required**. -For the Copernicus data, the minimum required GDAL version is 3.4.2; versions earlier than 3.4.0 seem to hang upon using `gdalwarp` on the global VRT, and <3.4.2 have an internal bug https://github.com/isce-framework/isce2/issues/556 . - +**Note:** To convert the elevation values to heights about the WGS84 ellipsoid (which is the default), or to use the Copernicus data, **rasterio is required**. +For the Copernicus data, rasterio >= 1.0.0 is required for coordinate transformations and data warping. ## Bounding box convention `sardem` uses the gdal convention ("pixel is area") where `--bbox` points to the *edges* of the [left, bottom, right, top] pixels. I.e. (left, bottom) refers to the lower left corner of the lower left pixel. - - -### Converting to WGS84 ellipsoidal heights from EGM96/EGM2008 geoid heights - -GDAL is required for the conversion, which is installed when using `conda install -c conda-forge sardem`. -If you already are using an existing environment, make sure that the GDAL version is >=3.4.2. - -```bash -conda install -c conda-forge "gdal>=3.4.2" - -# or -# conda install -c conda-forge mamba -# mamba install -c conda-forge "gdal>=3.4.2" -``` - - ## Command Line Interface The full options for the command line tool in `sardem/cli.py` can be found using @@ -83,7 +74,7 @@ Stiches SRTM .hgt files to make (upsampled) DEM the necessary SRTM1 tiles, stitch together, then upsample. The `--bbox` convention points to the *edges* of the [left, bottom, right, top] - pixels, following the "pixel is area" convention as used in gdal. + pixels, following the "pixel is area" convention as used in rasterio. I.e. (left, bottom) refers to the lower left corner of the lower left pixel. Usage Examples: diff --git a/environment.yml b/environment.yml index 7d4f41f..7030c0d 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ channels: - conda-forge dependencies: - python >= 3.6 -- gdal >= 3.4.2 +- rasterio >= 1.0.0 - numpy >= 1.14.2 - requests >= 2.20.0 - pip \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index afc4a22..e202bd5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,3 @@ numpy>=1.14.2 requests>=2.20.0 +rasterio>=1.0.0 diff --git a/sardem/cli.py b/sardem/cli.py index 431c2e7..1e3807b 100644 --- a/sardem/cli.py +++ b/sardem/cli.py @@ -1,6 +1,7 @@ """ Command line interface for running sardem """ + import json from argparse import ( ArgumentError, @@ -29,7 +30,7 @@ def positive_small_int(argstring): the necessary SRTM1 tiles, stitch together, then upsample. The `--bbox` convention points to the *edges* of the [left, bottom, right, top] - pixels, following the "pixel is area" convention as used in gdal. + pixels, following the "pixel is area" convention as used in rasterio. I.e. (left, bottom) refers to the lower left corner of the lower left pixel. Usage Examples: @@ -77,7 +78,7 @@ def get_cli_args(): help="Bounding box of area of interest " " (e.g. --bbox -106.1 30.1 -103.1 33.1 ). \n" "--bbox points to the *edges* of the pixels, \n" - " following the 'pixel is area' convention as used in gdal. ", + " following the 'pixel is area' convention as used in rasterio. ", ) parser.add_argument( "--geojson", diff --git a/sardem/conversions.py b/sardem/conversions.py index 6d5ee63..2797a4f 100644 --- a/sardem/conversions.py +++ b/sardem/conversions.py @@ -30,26 +30,41 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="eg # Target srs: WGS84 horizontal/vertical t_srs = '"epsg:4326"' - # Note: https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-tr - # If not specified, gdalwarp will generate an output raster with xsize=ysize - # We want it to match the input file + # Use rasterio for coordinate system transformation + import rasterio + from rasterio.warp import reproject, Resampling + xsize, ysize = _get_size(filename) - cmd = ( - 'gdalwarp {overwrite} -s_srs {s_srs} -t_srs {t_srs}' - ' -of ROI_PAC -ts {xsize} {ysize} ' - ' -multi -wo NUM_THREADS=4 -wm 4000 {inp} {out}' - ) - cmd = cmd.format( - inp=filename, - out=output, - overwrite="-overwrite" if overwrite else "", - xsize=xsize, - ysize=ysize, - s_srs=s_srs, - t_srs=t_srs, - ) - logger.info(cmd) - subprocess.run(cmd, check=True, shell=True) + + with rasterio.open(filename) as src: + # Create output profile + profile = src.profile.copy() + profile.update( + { + "driver": "ENVI", # ROI_PAC equivalent + "crs": t_srs.strip('"'), + "width": xsize, + "height": ysize, + } + ) + + if not overwrite and os.path.exists(output): + logger.warning("Output file {} exists and overwrite=False".format(output)) + return output + + logger.info("Converting {} from {} to {}".format(filename, s_srs, t_srs)) + with rasterio.open(output, "w", **profile) as dst: + for i in range(1, src.count + 1): + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=s_srs.strip('"'), + dst_transform=dst.transform, + dst_crs=t_srs.strip('"'), + resampling=Resampling.nearest, + num_threads=4, + ) if copy_rsc: rsc_file = filename + ".rsc" @@ -59,23 +74,19 @@ def egm_to_wgs84(filename, output=None, overwrite=True, copy_rsc=True, geoid="eg def _get_size(filename): - """Retrieve the raster size from a gdal-readable file""" - from osgeo import gdal + """Retrieve the raster size from a rasterio-readable file""" + import rasterio - ds = gdal.Open(filename) - xsize, ysize = ds.RasterXSize, ds.RasterYSize - ds = None + with rasterio.open(filename) as ds: + xsize, ysize = ds.width, ds.height return xsize, ysize def convert_dem_to_wgs84(dem_filename, geoid="egm96"): """Convert the file `dem_filename` from EGM96 heights to WGS84 ellipsoidal heights - Overwrites file, requires GDAL to be installed + Overwrites file, uses rasterio for conversion """ - if not utils._gdal_installed_correctly(): - logger.error("GDAL required to convert DEM to WGS84") - return path_, fname = os.path.split(dem_filename) rsc_filename = os.path.join(path_, fname + ".rsc") @@ -95,4 +106,4 @@ def convert_dem_to_wgs84(dem_filename, geoid="egm96"): logger.error("Failed to convert DEM:", exc_info=True) logger.error("Reverting back, using EGM dem as output") os.rename(output_egm, dem_filename) - os.rename(rsc_filename_egm, rsc_filename) \ No newline at end of file + os.rename(rsc_filename_egm, rsc_filename) diff --git a/sardem/cop_dem.py b/sardem/cop_dem.py index a205d7b..398a6f1 100644 --- a/sardem/cop_dem.py +++ b/sardem/cop_dem.py @@ -21,7 +21,7 @@ def download_and_stitch( yrate=1, vrt_filename=None, output_format="ENVI", - output_type="int16" + output_type="int16", ): """Download the COP DEM from AWS. @@ -32,9 +32,10 @@ def download_and_stitch( https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198 https://copernicus-dem-30m.s3.amazonaws.com/readme.html """ - from osgeo import gdal + import rasterio + import rasterio.warp + from rasterio.enums import Resampling - gdal.UseExceptions() # TODO: does downloading make it run any faster? # if download_vrt: # cache_dir = utils.get_cache_dir() @@ -56,51 +57,69 @@ def download_and_stitch( resamp = "bilinear" if (xrate > 1 or yrate > 1) else "nearest" # access_mode = "overwrite" if overwrite else None - option_dict = dict( - format=output_format, - outputBounds=bbox, - dstSRS=t_srs, - srcSRS=s_srs, - xRes=xres, - yRes=yres, - outputType=gdal.GetDataTypeByName(output_type.title()), - resampleAlg=resamp, - multithread=True, - warpMemoryLimit=5000, - warpOptions=["NUM_THREADS=4"], - ) - - # Used the __RETURN_OPTION_LIST__ to get the list of options for debugging + # Map resampling algorithm names + resamp_map = { + "bilinear": Resampling.bilinear, + "nearest": Resampling.nearest, + "cubic": Resampling.cubic, + "average": Resampling.average, + } + resample_alg = resamp_map.get(resamp, Resampling.bilinear) + logger.info("Creating {}".format(output_name)) logger.info("Fetching remote tiles...") - try: - cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_dict) - logger.info("Running GDAL command:") - logger.info(cmd) - except Exception: - # Can't form the cli version due to `deepcopy` Pickle error, just skip - logger.info("Running gdal.Warp with options:") - logger.info(option_dict) - pass - # Now convert to something GDAL can actually use - option_dict["callback"] = gdal.TermProgress - gdal.Warp(output_name, vrt_filename, options=gdal.WarpOptions(**option_dict)) - return - -def _gdal_cmd_from_options(src, dst, option_dict): - from osgeo import gdal - - opts = deepcopy(option_dict) - # To see what the list of cli options are (gdal >= 3.5.0) - opts["options"] = "__RETURN_OPTION_LIST__" - opt_list = gdal.WarpOptions(**opts) - out_opt_list = deepcopy(opt_list) - for idx, o in enumerate(opt_list): - # Wrap the srs option in quotes - if o.endswith("srs"): - out_opt_list[idx + 1] = '"{}"'.format(out_opt_list[idx + 1]) - return "gdalwarp {} {} {}".format(src, dst, " ".join(out_opt_list)) + # Use rasterio for warping + with rasterio.open(vrt_filename) as src: + # Calculate output transform and dimensions + dst_transform, dst_width, dst_height = ( + rasterio.warp.calculate_default_transform( + src.crs if s_srs is None else s_srs, + t_srs if t_srs is not None else src.crs, + src.width, + src.height, + *src.bounds, + resolution=(xres, yres), + ) + ) + + # If bbox is provided, override the transform calculation + if bbox: + left, bottom, right, top = bbox + dst_width = int(round((right - left) / xres)) + dst_height = int(round((top - bottom) / abs(yres))) + dst_transform = rasterio.transform.from_bounds( + left, bottom, right, top, dst_width, dst_height + ) + + # Create output profile + profile = src.profile.copy() + profile.update( + { + "driver": "ENVI" if output_format == "ENVI" else output_format, + "height": dst_height, + "width": dst_width, + "transform": dst_transform, + "crs": t_srs if t_srs is not None else src.crs, + "dtype": output_type, + } + ) + + warp_opts = {"XSCALE": 1, "YSCALE": 1, "APPLY_VERTICAL_SHIFT": True} + with rasterio.open(output_name, "w", **profile) as dst: + for i in range(1, src.count + 1): + rasterio.warp.reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs if s_srs is None else s_srs, + dst_transform=dst_transform, + dst_crs=t_srs if t_srs is not None else src.crs, + resampling=resample_alg, + num_threads=4, + **warp_opts, + ) + return def make_cop_vrt(outname="copernicus_GLO_30_dem.vrt"): @@ -108,22 +127,35 @@ def make_cop_vrt(outname="copernicus_GLO_30_dem.vrt"): Note: this is a large VRT file, ~15MB, so it can many hours to build. """ - from osgeo import gdal - - gdal.UseExceptions() + import rasterio + from rasterio.vrt import WarpedVRT tile_list = get_tile_list() url_list = _make_url_list(tile_list) - vrt_options = gdal.BuildVRTOptions( - resampleAlg=gdal.GRIORA_NearestNeighbour, - outputBounds=[-180, -90, 180, 90], - resolution="highest", - outputSRS="EPSG:4326+3855", - ) + logger.info("Building VRT {}".format(outname)) - vrt_file = gdal.BuildVRT(outname, url_list, options=vrt_options) - vrt_file.FlushCache() - vrt_file = None + + # Create a simple VRT by writing XML directly (rasterio doesn't have BuildVRT equivalent) + vrt_xml = '\n' + vrt_xml += '\n' + vrt_xml += " EPSG:4326+3855\n" + vrt_xml += " -1.800000000000000e+02, 2.777777777777778e-04, 0.000000000000000e+00, 8.333333333333334e+01, 0.000000000000000e+00, -2.777777777777778e-04\n" + vrt_xml += ' \n' + vrt_xml += " -32768\n" + + for url in url_list: + vrt_xml += f" \n" + vrt_xml += f' {url}\n' + vrt_xml += f" 1\n" + vrt_xml += f" \n" + + vrt_xml += " \n" + vrt_xml += "\n" + + with open(outname, "w") as f: + f.write(vrt_xml) + + logger.info(f"VRT written to {outname} with {len(url_list)} tiles") def get_tile_list(): diff --git a/sardem/dem.py b/sardem/dem.py index 63448da..64ce6ea 100644 --- a/sardem/dem.py +++ b/sardem/dem.py @@ -37,6 +37,7 @@ Z_SCALE 1 PROJECTION LL """ + from __future__ import division, print_function import collections @@ -313,11 +314,11 @@ def main( yrate (int): y-rate (rows) to upsample DEM (positive int) make_isce_xml (bool): whether to make an isce2-compatible XML file keep_egm (bool): Don't convert the DEM heights from geoid heights - above EGM96 or EGM2008 to heights above WGS84 ellipsoid + above EGM96 or EGM2008 to heights above WGS84 ellipsoid (default = False: do the conversion) shift_rsc (bool): Shift the .dem.rsc file down/right so that the X_FIRST and Y_FIRST values represent the pixel *center* (instead of - GDAL's convention of pixel edge). Default = False. + rasterio's convention of pixel edge). Default = False. cache_dir (str): directory to cache downloaded tiles output_type (str): output type for DEM (default = int16) output_format (str): output format for copernicus DEM (default = ENVI) @@ -346,9 +347,10 @@ def main( ) logger.warning("Are the bounds correct (left, bottom, right, top)?") - # For copernicus, use GDAL to warp from the VRT + # For copernicus, use rasterio to warp from the VRT if data_source == "COP": - utils._gdal_installed_correctly() + if not utils._rasterio_available(): + raise ImportError("rasterio is required for Copernicus DEM processing") from sardem import cop_dem cop_dem.download_and_stitch( @@ -419,17 +421,17 @@ def main( ) f.write(upsampled_rsc) - # Now upsample using with GDAL or python - if utils._gdal_installed_correctly(): - logger.info("Using GDAL Translate for upsampling") - upsample.upsample_with_gdal( + # Now upsample using rasterio or python + try: + logger.info("Using rasterio for upsampling") + upsample.upsample_with_rasterio( dem_filename_small, output_name, method="bilinear", # make this an option? xrate=xrate, yrate=yrate, ) - else: + except Exception: # Figure out size of row blocks to keep memory under 100 MB nrows, ncols = stitched_dem.shape block_rows = int(np.round(10e6 / ncols / 2, -2)) # round to 100s @@ -464,7 +466,7 @@ def main( else: logger.info("Correcting DEM to heights above WGS84 ellipsoid") conversions.convert_dem_to_wgs84(output_name, geoid="egm96") - + # If the user wants the .rsc file to point to pixel center: if shift_rsc: utils.shift_rsc_file(rsc_filename, to_gdal=False) diff --git a/sardem/download.py b/sardem/download.py index cc27eee..7ea1151 100644 --- a/sardem/download.py +++ b/sardem/download.py @@ -195,7 +195,7 @@ class Downloader: DATA_URLS = { "NASA": "https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11", "NASA_WATER": "https://e4ftl01.cr.usgs.gov/MEASURES/SRTMSWBD.003/2000.02.11", - "COP": "https://copernicus-dem-30m.s3.amazonaws.com/{t}/{t}.tif" + "COP": "https://copernicus-dem-30m.s3.amazonaws.com/{t}/{t}.tif", } VALID_SOURCES = DATA_URLS.keys() TILE_ENDINGS = { diff --git a/sardem/loading.py b/sardem/loading.py index 6080e09..6db884c 100644 --- a/sardem/loading.py +++ b/sardem/loading.py @@ -88,6 +88,7 @@ def load_watermask(filename): https://lpdaac.usgs.gov/products/srtmswbdv003/ """ import numpy as np + return np.fromfile(filename, dtype=np.uint8).reshape((3601, 3601)) diff --git a/sardem/tests/conftest.py b/sardem/tests/conftest.py index 6cbcee8..bb0152a 100644 --- a/sardem/tests/conftest.py +++ b/sardem/tests/conftest.py @@ -9,6 +9,7 @@ HALF_PIXEL = 0.5 * DEFAULT_RES DATA_PATH = os.path.join(os.path.dirname(__file__), "data") + @pytest.fixture def srtm_tile_path(tmp_path): return tmp_path / "N19W156.hgt" @@ -24,6 +25,7 @@ def srtm_tile(srtm_tile_path): expected[expected < -1000] = 0 return expected + @pytest.fixture def srtm_tile_bbox(): # $ rio bounds --bbox ~/.cache/sardem/N19W156.hgt @@ -34,4 +36,4 @@ def srtm_tile_bbox(): 19.0 + HALF_PIXEL, -155.0 - HALF_PIXEL, 20.0 + HALF_PIXEL, - ] \ No newline at end of file + ] diff --git a/sardem/tests/test_loading.py b/sardem/tests/test_loading.py index dcd2375..35dfed2 100644 --- a/sardem/tests/test_loading.py +++ b/sardem/tests/test_loading.py @@ -7,17 +7,27 @@ from sardem import loading -DATAPATH = join(dirname(__file__), 'data') +DATAPATH = join(dirname(__file__), "data") class TestRsc(unittest.TestCase): def setUp(self): - self.rsc_path = join(DATAPATH, 'elevation.dem.rsc') + self.rsc_path = join(DATAPATH, "elevation.dem.rsc") self.rsc_data = OrderedDict( - [('width', 2), ('file_length', 3), ('x_first', -155.676388889), ('y_first', - 19.5755555567), - ('x_step', 0.000138888888), ('y_step', -0.000138888888), ('x_unit', 'degrees'), - ('y_unit', 'degrees'), ('z_offset', 0), ('z_scale', 1), ('projection', 'LL')]) + [ + ("width", 2), + ("file_length", 3), + ("x_first", -155.676388889), + ("y_first", 19.5755555567), + ("x_step", 0.000138888888), + ("y_step", -0.000138888888), + ("x_unit", "degrees"), + ("y_unit", "degrees"), + ("z_offset", 0), + ("z_scale", 1), + ("projection", "LL"), + ] + ) self.extract_path = tempfile.mkdtemp() def tearDown(self): @@ -33,8 +43,8 @@ def test_format_dem_rsc(self): self.assertEqual(output, read_file) def test_load_elevation(self): - zip_ref = zipfile.ZipFile(join(DATAPATH, 'N19W156.hgt.zip')) + zip_ref = zipfile.ZipFile(join(DATAPATH, "N19W156.hgt.zip")) zip_ref.extractall(self.extract_path) zip_ref.close() - hgt_file = join(self.extract_path, 'N19W156.hgt') + hgt_file = join(self.extract_path, "N19W156.hgt") loading.load_elevation(hgt_file) diff --git a/sardem/tests/test_upsample.py b/sardem/tests/test_upsample.py index ac2ce0a..0b83b0d 100644 --- a/sardem/tests/test_upsample.py +++ b/sardem/tests/test_upsample.py @@ -38,7 +38,7 @@ def test_resample(): } bbox = (1, 1, 4, 4) a_resampled = upsample.resample(a, rsc_dict_mid, bbox) - assert_allclose(expected, a_resampled) + assert_allclose(expected, a_resampled) bbox = (0.5, 0.5, 4.5, 4.5) a_resampled = upsample.resample(a, rsc_dict_mid, bbox) @@ -59,6 +59,7 @@ def test_resample_tile(srtm_tile, srtm_tile_bbox): # should be the exact same assert_allclose(srtm_tile[:-1, :-1], tile_resamp) + def test_upsample_dem_rsc(): rsc_path = os.path.join(os.path.dirname(__file__), "data", "elevation.dem.rsc") # Test input checking diff --git a/sardem/upsample.py b/sardem/upsample.py index 85868dc..54a07ed 100644 --- a/sardem/upsample.py +++ b/sardem/upsample.py @@ -8,11 +8,8 @@ utils.set_logger_handler(logger) -def upsample_with_gdal(filename, outfile, method="cubic", xrate=1, yrate=1): - """Perform upsampling on a raster using gdal - - See here for available methods - https://gdal.org/programs/gdal_translate.html#cmdoption-gdal_translate-r +def upsample_with_rasterio(filename, outfile, method="cubic", xrate=1, yrate=1): + """Perform upsampling on a raster using rasterio Parameters ---------- @@ -21,7 +18,7 @@ def upsample_with_gdal(filename, outfile, method="cubic", xrate=1, yrate=1): outfile : str Name of output file method : str - Interpolation method to use, see above + Interpolation method to use (cubic, bilinear, nearest, etc.) xrate : int, optional. Rate to upsample in the x/column direction Default = 1, no upsampling @@ -29,16 +26,53 @@ def upsample_with_gdal(filename, outfile, method="cubic", xrate=1, yrate=1): Rate to upsample in the y/row direction Default = 1, no upsampling """ - from osgeo import gdal - - options = gdal.TranslateOptions( - format="ROI_PAC", - widthPct=xrate * 100, - heightPct=yrate * 100, - resampleAlg=method, - callback=gdal.TermProgress, - ) - gdal.Translate(outfile, filename, options=options) + import rasterio + from rasterio.enums import Resampling + from rasterio.warp import reproject + + # Map method names to rasterio resampling methods + resampling_methods = { + "cubic": Resampling.cubic, + "bilinear": Resampling.bilinear, + "nearest": Resampling.nearest, + "lanczos": Resampling.lanczos, + "average": Resampling.average, + } + + resample_alg = resampling_methods.get(method.lower(), Resampling.cubic) + + with rasterio.open(filename) as src: + # Calculate new dimensions + new_width = int(src.width * xrate) + new_height = int(src.height * yrate) + + # Calculate new transform + new_transform = src.transform * src.transform.scale( + (src.width / new_width), (src.height / new_height) + ) + + # Create output profile + profile = src.profile.copy() + profile.update( + { + "driver": "ENVI", # ROI_PAC equivalent + "height": new_height, + "width": new_width, + "transform": new_transform, + } + ) + + with rasterio.open(outfile, "w", **profile) as dst: + for i in range(1, src.count + 1): + reproject( + source=rasterio.band(src, i), + destination=rasterio.band(dst, i), + src_transform=src.transform, + src_crs=src.crs, + dst_transform=new_transform, + dst_crs=src.crs, + resampling=resample_alg, + ) def upsample_by_blocks( diff --git a/sardem/utils.py b/sardem/utils.py index f46515b..d6932dc 100644 --- a/sardem/utils.py +++ b/sardem/utils.py @@ -1,9 +1,10 @@ """utils.py: extra helper functions""" + from __future__ import division, print_function import logging import os -import subprocess +import importlib.util import sys from math import floor @@ -154,8 +155,8 @@ def shift_rsc_file(rsc_filename=None, outname=None, to_gdal=True): is needed to create a .rsc file in GDAL convention. See here for geotransform info - https://gdal.org/user/raster_data_model.html#affine-geotransform - GDAL standard is to reference a raster by its top left edges, + https://rasterio.readthedocs.io/en/latest/topics/georeferencing.html + rasterio standard is to reference a raster by its top left edges, while some SAR processors use the middle of a pixel. `to_gdal`=True means it moves the X_FIRST, Y_FIRST up and left half a pixel. @@ -181,8 +182,8 @@ def shift_rsc_dict(rsc_dict, to_gdal=True): is needed to create a .rsc file in GDAL convention. See here for geotransform info - https://gdal.org/user/raster_data_model.html#affine-geotransform - GDAL standard is to reference a raster by its top left edges, + https://rasterio.readthedocs.io/en/latest/topics/georeferencing.html + rasterio standard is to reference a raster by its top left edges, while some SAR processors use the middle of a pixel. `to_gdal`=True means it moves the X_FIRST, Y_FIRST up and left half a pixel. @@ -225,16 +226,9 @@ def get_output_size(bounds, xrate, yrate): return rows, cols -def _gdal_installed_correctly(): - cmd = "gdalinfo --help-general" - # cmd = "gdalinfo -h" - try: - subprocess.check_output(cmd, shell=True) - return True - except subprocess.CalledProcessError: - logger.error("GDAL command failed to run.", exc_info=True) - logger.error("Check GDAL installation.") - return False +def _rasterio_available(): + """Check if rasterio is available""" + return importlib.util.find_spec("rasterio") is not None def gdal2isce_xml(fname, keep_egm=False): @@ -245,8 +239,9 @@ def gdal2isce_xml(fname, keep_egm=False): from applications.gdal2isce_xml import gdal2isce_xml xml_file = gdal2isce_xml(fname+'.vrt') """ - _gdal_installed_correctly() - from osgeo import gdal + if not _rasterio_available(): + raise ImportError("rasterio is required for XML generation") + import rasterio try: import isce # noqa @@ -255,17 +250,17 @@ def gdal2isce_xml(fname, keep_egm=False): logger.error("isce2 not installed. Cannot generate xml file.") raise - # open the GDAL file and get typical data informationi - GDAL2ISCE_DATATYPE = { - 1: "BYTE", - 2: "uint16", - 3: "SHORT", - 4: "uint32", - 5: "INT", - 6: "FLOAT", - 7: "DOUBLE", - 10: "CFLOAT", - 11: "complex128", + # open the rasterio file and get typical data information + RASTERIO2ISCE_DATATYPE = { + "uint8": "BYTE", + "uint16": "uint16", + "int16": "SHORT", + "uint32": "uint32", + "int32": "INT", + "float32": "FLOAT", + "float64": "DOUBLE", + "complex64": "CFLOAT", + "complex128": "complex128", } # check if the input file is a vrt @@ -276,25 +271,21 @@ def gdal2isce_xml(fname, keep_egm=False): outname = fname logger.info("Writing to %s", outname) - # open the GDAL file and get typical ds information - ds = gdal.Open(fname, gdal.GA_ReadOnly) - width = ds.RasterXSize - length = ds.RasterYSize - bands = ds.RasterCount - logger.info("width: " + "\t" + str(width)) - logger.info("length: " + "\t" + str(length)) - logger.info("num of bands:" + "\t" + str(bands)) - - # getting the datatype information - raster = ds.GetRasterBand(1) - dataTypeGdal = raster.DataType - - # user look-up dictionary from gdal to isce format - dataType = GDAL2ISCE_DATATYPE[dataTypeGdal] - logger.info("dataType: " + "\t" + str(dataType)) - - # transformation contains gridcorners (lines/pixels or lonlat and the spacing 1/-1 or deltalon/deltalat) - transform = ds.GetGeoTransform() + # open the rasterio file and get typical dataset information + with rasterio.open(fname) as ds: + width = ds.width + length = ds.height + bands = ds.count + logger.info("width: " + "\t" + str(width)) + logger.info("length: " + "\t" + str(length)) + logger.info("num of bands:" + "\t" + str(bands)) + + # getting the datatype information + dataType = RASTERIO2ISCE_DATATYPE.get(str(ds.dtypes[0]), str(ds.dtypes[0])) + logger.info("dataType: " + "\t" + str(dataType)) + + # transformation contains gridcorners (lines/pixels or lonlat and the spacing 1/-1 or deltalon/deltalat) + transform = ds.transform # if a complex data type, then create complex image # if a real data type, then create a regular image @@ -307,28 +298,18 @@ def gdal2isce_xml(fname, keep_egm=False): img.dataType = dataType # interleave - md = ds.GetMetadata("IMAGE_STRUCTURE") - sch = md.get("INTERLEAVE", None) - if sch == "LINE": - img.scheme = "BIL" - elif sch == "PIXEL": + if bands < 2: + logger.info("Single band, using BIP") img.scheme = "BIP" - elif sch == "BAND": - img.scheme = "BSQ" else: - logger.info("Unrecognized interleaving scheme, {}".format(sch)) - if bands < 2: - logger.info("Assuming default, BIP") - img.scheme = "BIP" - else: - logger.info("Assuming default, BSQ") - img.scheme = "BSQ" - - first_lon = transform[0] - first_lat = transform[3] - delta_lat = transform[5] - delta_lon = transform[1] - # We are using gdal conventions of the edges of the image for the bbox + logger.info("Multi-band, using BSQ") + img.scheme = "BSQ" + + first_lon = transform.c + first_lat = transform.f + delta_lat = transform.e + delta_lon = transform.a + # We are using rasterio conventions of the edges of the image for the bbox # We need to shift the `first________` values to the middle of the top left pixel first_lon += 0.5 * delta_lon first_lat += 0.5 * delta_lat diff --git a/setup.py b/setup.py index fa670e1..a4961bc 100644 --- a/setup.py +++ b/setup.py @@ -25,7 +25,7 @@ "Topic :: Scientific/Engineering", "Intended Audience :: Science/Research", ], - install_requires=["numpy", "requests"], + install_requires=["numpy", "requests", "rasterio>=1.0.0"], entry_points={ "console_scripts": [ "sardem=sardem.cli:cli",