|
| 1 | +import logging |
| 2 | +import netrc |
| 3 | +import os |
| 4 | +from copy import deepcopy |
| 5 | + |
| 6 | +from sardem import utils |
| 7 | +from sardem.constants import DEFAULT_RES |
| 8 | + |
| 9 | +_NISAR_BASE_URL = "https://nisar.asf.earthdatacloud.nasa.gov/NISAR/DEM/v1.2" |
| 10 | +NISAR_VRTS = { |
| 11 | + "EPSG4326": f"{_NISAR_BASE_URL}/EPSG4326/EPSG4326.vrt", |
| 12 | + "EPSG3031": f"{_NISAR_BASE_URL}/EPSG3031/EPSG3031.vrt", # South pole |
| 13 | + "EPSG3413": f"{_NISAR_BASE_URL}/EPSG3413/EPSG3413.vrt", # North pole |
| 14 | +} |
| 15 | +# Latitude thresholds for switching to polar stereographic VRTs |
| 16 | +_SOUTH_POLE_LAT = -60.0 |
| 17 | +_NORTH_POLE_LAT = 60.0 |
| 18 | + |
| 19 | +EARTHDATA_HOST = "urs.earthdata.nasa.gov" |
| 20 | + |
| 21 | +logger = logging.getLogger("sardem") |
| 22 | +utils.set_logger_handler(logger) |
| 23 | + |
| 24 | + |
| 25 | +def _check_earthdata_credentials() -> None: |
| 26 | + """Verify that ~/.netrc has credentials for urs.earthdata.nasa.gov. |
| 27 | +
|
| 28 | + Raises |
| 29 | + ------ |
| 30 | + RuntimeError |
| 31 | + If credentials are missing or incomplete. |
| 32 | + """ |
| 33 | + netrc_path = os.path.expanduser("~/.netrc") |
| 34 | + try: |
| 35 | + nrc = netrc.netrc(netrc_path) |
| 36 | + auth = nrc.authenticators(EARTHDATA_HOST) |
| 37 | + assert auth is not None and auth[0] and auth[2] |
| 38 | + except (OSError, netrc.NetrcParseError, AssertionError): |
| 39 | + raise RuntimeError( |
| 40 | + "NASA Earthdata credentials not found in ~/.netrc for" |
| 41 | + f" {EARTHDATA_HOST}.\n" |
| 42 | + "The NISAR DEM requires a free NASA Earthdata account.\n" |
| 43 | + "Sign up at: https://urs.earthdata.nasa.gov/users/new\n" |
| 44 | + "Then add to ~/.netrc:\n" |
| 45 | + f" machine {EARTHDATA_HOST}\n" |
| 46 | + " login <username>\n" |
| 47 | + " password <password>" |
| 48 | + ) |
| 49 | + |
| 50 | + |
| 51 | +def _configure_gdal_auth() -> None: |
| 52 | + """Set GDAL config options for cookie-based NASA Earthdata authentication.""" |
| 53 | + from osgeo import gdal |
| 54 | + |
| 55 | + cookie_file = os.path.join(utils.get_cache_dir(), "earthdata_cookies.txt") |
| 56 | + gdal.SetConfigOption("GDAL_HTTP_COOKIEFILE", cookie_file) |
| 57 | + gdal.SetConfigOption("GDAL_HTTP_COOKIEJAR", cookie_file) |
| 58 | + gdal.SetConfigOption("GDAL_HTTP_AUTH", "BASIC") |
| 59 | + gdal.SetConfigOption("GDAL_HTTP_NETRC", "YES") |
| 60 | + |
| 61 | + |
| 62 | +def _select_vrt(bbox: tuple) -> tuple[str, str | None]: |
| 63 | + """Choose the appropriate NISAR VRT and output SRS based on bbox latitude. |
| 64 | +
|
| 65 | + Parameters |
| 66 | + ---------- |
| 67 | + bbox : tuple |
| 68 | + (left, bottom, right, top) in degrees. |
| 69 | +
|
| 70 | + Returns |
| 71 | + ------- |
| 72 | + vrt_url : str |
| 73 | + The ``/vsicurl/`` URL to the selected VRT. |
| 74 | + dst_srs : str or None |
| 75 | + Target SRS for ``gdal.Warp``. ``None`` when the VRT is already |
| 76 | + EPSG:4326; ``"EPSG:4326"`` when reprojecting from a polar VRT. |
| 77 | + """ |
| 78 | + _left, bottom, _right, top = bbox |
| 79 | + if top <= _SOUTH_POLE_LAT: |
| 80 | + key = "EPSG3031" |
| 81 | + elif bottom >= _NORTH_POLE_LAT: |
| 82 | + key = "EPSG3413" |
| 83 | + else: |
| 84 | + key = "EPSG4326" |
| 85 | + url = "/vsicurl/" + NISAR_VRTS[key] |
| 86 | + dst_srs = "EPSG:4326" if key != "EPSG4326" else None |
| 87 | + logger.info("Selected NISAR VRT: %s", key) |
| 88 | + return url, dst_srs |
| 89 | + |
| 90 | + |
| 91 | +def download_and_stitch( |
| 92 | + output_name: str, |
| 93 | + bbox: tuple, |
| 94 | + xrate: int = 1, |
| 95 | + yrate: int = 1, |
| 96 | + vrt_filename: str | None = None, |
| 97 | + output_format: str = "GTiff", |
| 98 | + output_type: str = "float32", |
| 99 | +) -> None: |
| 100 | + """Download the NISAR DEM via its global VRT. |
| 101 | +
|
| 102 | + The NISAR DEM is a Copernicus-derived DEM prepared by JPL with |
| 103 | + EGM2008-to-WGS84 conversion pre-applied and ocean gaps filled. |
| 104 | + No vertical datum conversion is needed. |
| 105 | +
|
| 106 | + Parameters |
| 107 | + ---------- |
| 108 | + output_name : str |
| 109 | + Path for the output DEM file. |
| 110 | + bbox : tuple |
| 111 | + (left, bottom, right, top) in degrees. |
| 112 | + xrate : int |
| 113 | + Column upsampling rate. |
| 114 | + yrate : int |
| 115 | + Row upsampling rate. |
| 116 | + vrt_filename : str, optional |
| 117 | + Override the VRT source (for testing). Defaults to the NISAR VRT URL. |
| 118 | + output_format : str |
| 119 | + GDAL output format (default ``"GTiff"``). |
| 120 | + output_type : str |
| 121 | + GDAL output data type (default ``"float32"``). |
| 122 | +
|
| 123 | + References |
| 124 | + ---------- |
| 125 | + https://nisar.asf.earthdatacloud.nasa.gov/NISAR/DEM/v1.2/EPSG4326/ |
| 126 | + """ |
| 127 | + from osgeo import gdal |
| 128 | + |
| 129 | + gdal.UseExceptions() |
| 130 | + |
| 131 | + dst_srs = None |
| 132 | + if vrt_filename is None: |
| 133 | + _check_earthdata_credentials() |
| 134 | + _configure_gdal_auth() |
| 135 | + vrt_filename, dst_srs = _select_vrt(bbox) |
| 136 | + |
| 137 | + xres = DEFAULT_RES / xrate |
| 138 | + yres = DEFAULT_RES / yrate |
| 139 | + resamp = "bilinear" if (xrate > 1 or yrate > 1) else "nearest" |
| 140 | + |
| 141 | + option_dict = dict( |
| 142 | + format=output_format, |
| 143 | + outputBounds=utils.align_bounds_to_pixel_grid(bbox), |
| 144 | + dstSRS=dst_srs, |
| 145 | + xRes=xres, |
| 146 | + yRes=yres, |
| 147 | + outputType=gdal.GetDataTypeByName(output_type.title()), |
| 148 | + resampleAlg=resamp, |
| 149 | + multithread=True, |
| 150 | + warpMemoryLimit=5000, |
| 151 | + warpOptions=["NUM_THREADS=4"], |
| 152 | + ) |
| 153 | + |
| 154 | + logger.info("Creating %s", output_name) |
| 155 | + logger.info("Fetching remote tiles...") |
| 156 | + try: |
| 157 | + cmd = _gdal_cmd_from_options(vrt_filename, output_name, option_dict) |
| 158 | + logger.info("Running GDAL command:") |
| 159 | + logger.info(cmd) |
| 160 | + except Exception: |
| 161 | + logger.info("Running gdal.Warp with options:") |
| 162 | + logger.info(option_dict) |
| 163 | + |
| 164 | + option_dict["callback"] = gdal.TermProgress |
| 165 | + gdal.Warp(output_name, vrt_filename, options=gdal.WarpOptions(**option_dict)) |
| 166 | + |
| 167 | + |
| 168 | +def _gdal_cmd_from_options(src: str, dst: str, option_dict: dict) -> str: |
| 169 | + """Build an equivalent ``gdalwarp`` CLI string for debugging.""" |
| 170 | + from osgeo import gdal |
| 171 | + |
| 172 | + opts = deepcopy(option_dict) |
| 173 | + opts["options"] = "__RETURN_OPTION_LIST__" |
| 174 | + opt_list = gdal.WarpOptions(**opts) |
| 175 | + out_opt_list = deepcopy(opt_list) |
| 176 | + for idx, o in enumerate(opt_list): |
| 177 | + if o.endswith("srs"): |
| 178 | + out_opt_list[idx + 1] = '"{}"'.format(out_opt_list[idx + 1]) |
| 179 | + return "gdalwarp {} {} {}".format(src, dst, " ".join(out_opt_list)) |
0 commit comments