def upscale (src_file_path, dest_file_path, upscale_factor=0.5):
""" This function upscales a dataset.
In the current context, upscaling refers to resampling to lower resolution/larger cell sizes.
Note that upscaling can also be referred to as downsampling.
src_file_path: path of the source dataset which will be upscaled
dest_file_path: pile of the destintation where the upscaled dataset will be saved
upscale_factor: multipicative factor by which the dataset size will be reduced
Resampling based on https://rasterio.readthedocs.io/en/latest/topics/resampling.html
"""
with rio.open(src_file_path) as dataset:
# resample data to target shape
data = dataset.read(
out_shape=(
dataset.count,
int(dataset.height * upscale_factor),
int(dataset.width * upscale_factor)
),
resampling=Resampling.bilinear
)
### scale image transform
transform = dataset.transform * dataset.transform.scale(
(dataset.width / data.shape[-1]),
(dataset.height / data.shape[-2])
)
out_meta = dataset.meta.copy()
### Update the metadata to write the upscaled dataset
out_meta.update({"driver": "GTiff",
"height": data.shape[1],
"width": data.shape[2],
"transform": transform
}
)
with rio.open(dest_file_path, "w", **out_meta) as dest:
dest.write(data)
def get_Points_gdf(src_file):
"""This function reads in a raster, and obtains Long/Lat values of the center
of each pixel. It returns a geodataframe."""
raster = rio.open(src_file)
arr = raster.read(1)
indices = np.argwhere(~np.isnan(arr))
coords = raster.xy(indices[:,0], indices[:,1])
Points_gdf= gpd.GeoDataFrame({'Longitude':coords[0],'Latitude': coords[1]},
geometry=[shapely.geometry.Point(x, y) for x, y in zip(coords[0], coords[1])])
return Points_gdf
def get_nearest_neighbors (this_date_string):
"""Take in a string value (this_date_string) which represents the date.
Return nearest neighbor values (w_1 to w_10) for that date"""
n_stations = top10_stations_df.shape[1]
n_pixels = top10_stations_df.shape[0]
date_index = Dates[Dates == this_date_string].index
col_names = ['Station_' + str(index) for index in np.arange(n_stations)]
nearest_neighbors = pd.DataFrame(index=np.arange(n_pixels), columns=col_names).astype('float32')
for ind in np.arange(n_pixels):
nearest_neighbors.iloc[ind,:] = Point_precip.loc[date_index,top10_stations_df.loc[ind,:]].values
return nearest_neighbors
def unflatten_index(flat_ind, ncols):
"""Take in the index of a flattened array (row-wise) and return the row and column index of an unflattened array"""
row_ind = flat_ind//ncols ## quotient
col_ind = flat_ind%ncols ## remainder
return (row_ind, col_ind)
def flatten_index (row_ind, col_ind, ncols):
"""Take in the row and column index and return the index of a flattened array that goes row-wise."""
return row_ind*ncols + col_ind
def get_800_rowcols (scale, row_ind_sub800, col_ind_sub800):
"""Take in (row,col) index corresponding to a sub-800m (here, 400 m) grid, and return (row,col) index for 800m-res grid.
The smaller grid is assumed to fit perfectly within the larger grid."""
row_ind_800 = row_ind_sub800//scale
col_ind_800 = col_ind_sub800//scale
return (row_ind_800, col_ind_800)
def get_unflattened_800res_index (flattened_index_sub800, ncols_sub800, ncols_800, scale):
## Get the indices of row,col for the regular unflattened shape
unflattened_ind_sub800 = unflatten_index(flattened_index_sub800, ncols_sub800)
## Get the indices of row,col for the 800m-res input
unflattened_ind_800 = get_800_rowcols(scale, *unflattened_ind_sub800)
## Get the index corresponding to the flattened 800m-res file
flattened_index_800 = flatten_index(*unflattened_ind_800, ncols_800)
return flattened_index_800
def get_downscaling_features(Points_sub800_gdf, Features_topography_sub800_df, PRISM_800_upscaled_file, precip_corr):
with rio.open(PRISM_800_upscaled_file) as src:
ncols_800 = src.read(1).shape[1]
ncols_sub800 = int(ncols_800*downscale_factor)
parent_indices = get_unflattened_800res_index (Points_sub800_gdf.index, ncols_sub800, ncols_800, downscale_factor)
precip_corr_sub800 = precip_corr.loc[parent_indices]
precip_corr_sub800.reset_index(drop=True, inplace=True)
Features_sub800_df = pd.merge(Features_topography_sub800_df, precip_corr_sub800, left_index=True, right_index=True)
return Features_sub800_df
def get_smooth_value (row_ind,col_ind, noisy_PRISM, height, width):
x_cell = ds_grid_x[row_ind,col_ind]
y_cell = ds_grid_y[row_ind,col_ind]
value_cell = noisy_PRISM[row_ind,col_ind]
bottomrow_ind = min(row_ind + max_cells + 1, height)
toprow_ind = max(row_ind - max_cells, 0)
leftcol_ind = max(col_ind - max_cells, 0)
rightcol_ind = min(col_ind + max_cells + 1, width)
x_grid_trunc = ds_grid_x[toprow_ind:bottomrow_ind, leftcol_ind:rightcol_ind]
y_grid_trunc = ds_grid_y[toprow_ind:bottomrow_ind, leftcol_ind:rightcol_ind]
noisy_PRISM_trunc = noisy_PRISM[toprow_ind:bottomrow_ind, leftcol_ind:rightcol_ind]
distances = euclidean_distances(np.array([x_cell, y_cell])[np.newaxis], np.array([x_grid_trunc.flatten(), y_grid_trunc.flatten()]).T).flatten()*downscale_res ##distance in kilometers
distances = np.reshape(distances, x_grid_trunc.shape)
neighbor_values = noisy_PRISM_trunc[np.logical_and(distances <= max_dist, distances > 0)]
neighbor_distances = distances[np.logical_and(distances <= max_dist, distances > 0)]
mean_abs_diff = np.mean(np.abs(value_cell - neighbor_values))
mean_abs_diff_max = value_cell*0.04
a = a_max if mean_abs_diff >= mean_abs_diff_max else a_max*mean_abs_diff/mean_abs_diff_max
value_cell_smoothed = np.sum(neighbor_values/(neighbor_distances**a))/np.sum(1.0/(neighbor_distances**a))
return value_cell_smoothed
