diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py index 85731aeee..87b69c19f 100755 --- a/tools/RAiDER/cli/validators.py +++ b/tools/RAiDER/cli/validators.py @@ -72,7 +72,7 @@ def get_heights(args, out, station_file, bounding_box=None): elif os.path.exists(args.dem): out['dem'] = args['dem'] if bounding_box is not None: - dem_bounds = rio_extents(rio_profile(args.dem)) + dem_bounds = rio_extents(rio_profile(args.dem)[0]) lats = dem_bounds[:2] lons = dem_bounds[2:] if isOutside( diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py index 2329ac2ae..0df558856 100755 --- a/tools/RAiDER/delay.py +++ b/tools/RAiDER/delay.py @@ -118,7 +118,7 @@ def tropo_delay(dt, wetFilename, hydroFilename, args): ########################################################### # Load the downloaded model file for CRS information - wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"] + wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")[0]["crs"] if wm_proj is None: print("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84") wm_proj = CRS.from_epsg(4326) @@ -282,7 +282,7 @@ def tropo_delay_cube(dt, wf, args, model_file=None): logger.debug(f'Output height range is {min(heights)} to {max(heights)}') # Load CRS from weather model file - wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"] + wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")[0]["crs"] if wm_proj is None: print("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84") wm_proj = CRS.from_epsg(4326) diff --git a/tools/RAiDER/llreader.py b/tools/RAiDER/llreader.py index d72566669..52c323ba2 100644 --- a/tools/RAiDER/llreader.py +++ b/tools/RAiDER/llreader.py @@ -88,6 +88,7 @@ def __init__(self, lat_file, lon_file=None, hgt_file=None, convention='isce'): else: self._latfile = lat_file self._lonfile = lon_file + self._gt = None self._proj, self._bounding_box, _ = bounds_from_latlon_rasters(lat_file, lon_file) # keep track of the height file @@ -131,11 +132,10 @@ class GeocodedFile(AOI): '''Parse a Geocoded file for coordinates''' def __init__(self, filename, is_dem=False): AOI.__init__(self) - self._filename = filename - self.p = rio_profile(filename) - self._bounding_box = rio_extents(self.p) - self._is_dem = is_dem - _, self._proj, self._gt = rio_stats(filename) + self._is_dem = is_dem + self.p, self._filename = rio_profile(filename) + self._bounding_box = rio_extents(self.p) + _, self._proj, self._gt = rio_stats(self._filename) def type(self): diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py index 9d183fae5..fe285a173 100644 --- a/tools/RAiDER/models/hrrr.py +++ b/tools/RAiDER/models/hrrr.py @@ -134,7 +134,7 @@ def _makeDataCubes(self, filename, out=None): out = filename # Get profile information from gdal - prof = rio_profile(str(filename)) + prof, filename = rio_profile(str(filename)) # Now get bounds S, N, W, E = self._ll_bounds diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py index dc9204000..f5736193e 100755 --- a/tools/RAiDER/utilFcns.py +++ b/tools/RAiDER/utilFcns.py @@ -23,6 +23,9 @@ from RAiDER.logger import logger +pbar = None + + def projectDelays(delay, inc): '''Project zenith delays to LOS''' return delay / cosd(inc) @@ -120,7 +123,13 @@ def rio_profile(fname): f"{fname} does not contain geotransform information" ) - return profile + elif profile['crs'].to_epsg() != 4326: + fname_wgs84 = rio_reproject_geocoded(fname) + with rasterio.open(fname_wgs84) as src: + profile = src.profile + fname = fname_wgs84 + + return profile, fname def rio_extents(profile): @@ -244,7 +253,7 @@ def writeArrayToRaster(array, filename, noDataValue=0., fmt='ENVI', proj=None, g ''' array_shp = np.shape(array) if array.ndim != 2: - raise RuntimeError('writeArrayToRaster: cannot write an array of shape {} to a raster image'.format(array_shp)) + raise RuntimeError(f'writeArrayToRaster: cannot write an array of shape {array_shp} to a raster image') # Data type if "complex" in str(array.dtype): @@ -439,6 +448,12 @@ def writeDelays(aoi, wetDelay, hydroDelay, lats, lons, wetDelay[np.isnan(wetDelay)] = ndv hydroDelay[np.isnan(hydroDelay)] = ndv + try: + proj, gt = aoi._proj, aoi._gt + except: + proj, gt = None, None + + # Do different things, depending on the type of input if aoi.type() == 'station_file': try: @@ -451,6 +466,7 @@ def writeDelays(aoi, wetDelay, hydroDelay, lats, lons, df['totalDelay'] = wetDelay + hydroDelay df.to_csv(wetFilename, index=False) + elif outformat == 'hdf5': writeResultsToHDF5( lats, @@ -966,9 +982,6 @@ def read_EarthData_loginInfo(filepath=None): return urs_usr, urs_pwd -pbar = None - - def show_progress(block_num, block_size, total_size): global pbar if pbar is None: @@ -1075,3 +1088,37 @@ def calcgeoh(lnsp, t, q, z, a, b, R_d, num_levels): z_h += TRd * dlogP return geopotential, pressurelvs, geoheight + + +def rio_reproject_geocoded(fname, epsg=4326): + """ Reproject the input geocoded file to 4326 + + Write a new geocoded file with with suffix .WGS84 + Required so weather model bounds align + https://rasterio.readthedocs.io/en/latest/topics/reproject.html + """ + dst_crs = f'EPSG:{epsg}' + dst_fname = f'{fname}.WGS84' + + with rasterio.open(fname) as src: + transform, width, height = rasterio.warp.calculate_default_transform( + src.crs, dst_crs, src.width, src.height, *src.bounds) + kwargs = src.meta.copy() + kwargs.update({ + 'crs': dst_crs, + 'transform': transform, + 'width': width, + 'height': height + }) + + with rasterio.open(dst_fname, 'w', **kwargs) 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, + dst_transform=transform, + dst_crs=dst_crs, + resampling=rasterio.warp.Resampling.nearest) + return dst_fname