Skip to content

return delay files in geo coords if geocoded aoi #408

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
bbuzzanga wants to merge 2 commits into
base: dev
Choose a base branch
from
Draft

return delay files in geo coords if geocoded aoi #408

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e0a3e69
return delay files in geo coords if geocoded aoi
bbuzzanga Nov 23, 2022
4ba25a1
support geocoded file in any coordinates
bbuzzanga Nov 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion tools/RAiDER/cli/validators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Copy link
Collaborator

@jlmaurer jlmaurer Nov 28, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we add a try/except here so either will work? Or is it always [0]?

dem_bounds = rio_extents(rio_profile(args.dem)[0])
lats = dem_bounds[:2]
lons = dem_bounds[2:]
if isOutside(
Expand Down
4 changes: 2 additions & 2 deletions tools/RAiDER/delay.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Copy link
Collaborator

@jlmaurer jlmaurer Nov 28, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same thing here as above, and below as well.

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)
Expand Down Expand Up @@ -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)
Expand Down
10 changes: 5 additions & 5 deletions tools/RAiDER/llreader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down
2 changes: 1 addition & 1 deletion tools/RAiDER/models/hrrr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
57 changes: 52 additions & 5 deletions tools/RAiDER/utilFcns.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@
from RAiDER.logger import logger


pbar = None


def projectDelays(delay, inc):
'''Project zenith delays to LOS'''
return delay / cosd(inc)
Expand Down Expand Up @@ -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)
Copy link
Collaborator

@jlmaurer jlmaurer Nov 28, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this the part that wasn't working to correctly reproject the outputs?

with rasterio.open(fname_wgs84) as src:
profile = src.profile
fname = fname_wgs84

return profile, fname


def rio_extents(profile):
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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:
Expand All @@ -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,
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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