Skip to content

Commit

Permalink
remove spatial interpolation method to separate package geostatista
Browse files Browse the repository at this point in the history
  • Loading branch information
@MAfarrag
MAfarrag committed Jun 6, 2022
1 parent f21d44c commit 9027c58
Showing 1 changed file with 0 additions and 213 deletions.
213 changes: 0 additions & 213 deletions statista/tools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,219 +15,6 @@ class Tools:
def __init__(self):
pass

@staticmethod
def IDW(raster, coordinates, data, No_data_cells=False):
"""IDW.
this function generates distributred values from reading at stations
using inverse distance weighting method
Parameters
raster:
GDAL calss represent the GIS raster
coordinates:
dict {'x':[],'y':[]} with two columns contains x coordinates and y
coordinates of the stations
data:
numpy array contains values of the timeseries at each gauge in the
same order as the coordinates in the coordinates lists (x,y)
No_data_cells:
boolen value (True or False) if the user want to calculate the
values in the cells that has no data value (cropped) No_data_cells
equal True if not No_data_cells equals False (default is false )
Returns
-------
array:
array with the same dimension of the raster
"""
# get the shaoe of the raster
shape_base_dem = raster.ReadAsArray().shape
# get the coordinates of the top left corner and cell size [x,dx,y,dy]
geo_trans = raster.GetGeoTransform()
# get the no_value
no_val = np.float32(
raster.GetRasterBand(1).GetNoDataValue()
) # get the value stores in novalue cells
# read the raster
raster_array = raster.ReadAsArray()
# calculate the coordinates of the center of each cell
# X_coordinate= upperleft corner x+ index* cell size+celsize/2
coox = np.ones(shape_base_dem)
cooy = np.ones(shape_base_dem)
# calculate the coordinates of the cells
if (
not No_data_cells
): # if user decide not to calculate values in the o data cells
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if raster_array[i, j] != no_val:
coox[i, j] = (
geo_trans[0] + geo_trans[1] / 2 + j * geo_trans[1]
) # calculate x
cooy[i, j] = (
geo_trans[3] + geo_trans[5] / 2 + i * geo_trans[5]
) # calculate y
else:
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
coox[i, j] = (
geo_trans[0] + geo_trans[1] / 2 + j * geo_trans[1]
) # calculate x
cooy[i, j] = (
geo_trans[3] + geo_trans[5] / 2 + i * geo_trans[5]
) # calculate y

# inverse the distance from the cell to each station
inverseDist = np.ones(
(shape_base_dem[0], shape_base_dem[1], len(coordinates["x"]))
)
# denominator of the equation (sum(1/d))
denominator = np.ones((shape_base_dem[0], shape_base_dem[1]))

for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if not np.isnan(coox[i, j]): # calculate only if the coox is not nan
for st in range(
len(coordinates["x"])
): # iteration by station [0]: #
# distance= sqrt((xstation-xcell)**2+ystation-ycell)**2)
inverseDist[i, j, st] = 1 / (
np.sqrt(
np.power((coordinates["x"][st] - coox[i, j]), 2)
+ np.power((coordinates["y"][st] - cooy[i, j]), 2)
)
)
denominator = np.sum(inverseDist, axis=2)

sp_dist = (
np.ones((len(raster[:, 0]), shape_base_dem[0], shape_base_dem[1])) * np.nan
)

for t in range(len(raster[:, 0])): # iteration by time step
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if not np.isnan(
coox[i, j]
): # calculate only if the coox is not nan
sp_dist[i, j, t] = (
np.sum(inverseDist[i, j, :] * raster[t, :])
/ denominator[i, j]
)
# change the type to float 32
sp_dist = sp_dist.astype(np.float32)

return sp_dist


@staticmethod
def ISDW(raster, coordinates, data, No_data_cells=False):
"""ISDW.
this function generates distributred values from reading at stations using
inverse squared distance weighting method
Parameters
----------
raster: [dataset]
GDAL calss represent the GIS raster
coordinates: [dict]
dict {'x':[],'y':[]} with two columns contains x coordinates and y
coordinates of the stations
data: [array]
numpy array contains values of the timeseries at each gauge in the
same order as the coordinates in the coordinates lists (x,y)
No_data_cells: []
boolen value (True or False) if the user want to calculate the
values in the cells that has no data value (cropped) No_data_cells
equal True if not No_data_cells equals False
(default is false )
Returns
-------
sp_dist: [array]
array with the same dimension of the raster
"""
shape_base_dem = raster.ReadAsArray().shape
# get the coordinates of the top left corner and cell size [x,dx,y,dy]
geo_trans = raster.GetGeoTransform()
# get the no_value
no_val = np.float32(
raster.GetRasterBand(1).GetNoDataValue()
) # get the value stores in novalue cells
# read the raster
raster_array = raster.ReadAsArray()
# calculate the coordinates of the center of each cell
# X_coordinate= upperleft corner x+ index* cell size+celsize/2
coox = np.ones(shape_base_dem) * np.nan
cooy = np.ones(shape_base_dem) * np.nan
# calculate the coordinates of the cells
if (
not No_data_cells
): # if user decide not to calculate values in the o data cells
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if raster_array[i, j] != no_val:
coox[i, j] = (
geo_trans[0] + geo_trans[1] / 2 + j * geo_trans[1]
) # calculate x
cooy[i, j] = (
geo_trans[3] + geo_trans[5] / 2 + i * geo_trans[5]
) # calculate y
else:
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
coox[i, j] = (
geo_trans[0] + geo_trans[1] / 2 + j * geo_trans[1]
) # calculate x
cooy[i, j] = (
geo_trans[3] + geo_trans[5] / 2 + i * geo_trans[5]
) # calculate y

print("step 1 finished")
# inverse the distance from the cell to each station
inverseDist = np.ones(
(shape_base_dem[0], shape_base_dem[1], len(coordinates["x"]))
)
# denominator of the equation (sum(1/d))
denominator = np.ones((shape_base_dem[0], shape_base_dem[1]))

for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if not np.isnan(coox[i, j]): # calculate only if the coox is not nan
for st in range(
len(coordinates["x"])
): # iteration by station [0]: #
inverseDist[i, j, st] = (
1
/ (
np.sqrt(
np.power((coordinates["x"][st] - coox[i, j]), 2)
+ np.power((coordinates["y"][st] - cooy[i, j]), 2)
)
)
** 2
)
print("step 2 finished")
denominator = np.sum(inverseDist, axis=2)
sp_dist = (
np.ones((shape_base_dem[0], shape_base_dem[1], len(data[:, 0]))) * np.nan
)
for t in range(len(data[:, 0])): # iteration by time step
for i in range(shape_base_dem[0]): # iteration by row
for j in range(shape_base_dem[1]): # iteration by column
if not np.isnan(
coox[i, j]
): # calculate only if the coox is not nan
sp_dist[i, j, t] = (
np.sum(inverseDist[i, j, :] * data[t, :])
/ denominator[i, j]
)
# change the type to float 32
sp_dist = sp_dist.astype(np.float32)

return sp_dist

@staticmethod
def normalize(x):
"""Normalizer
Expand Down

0 comments on commit 9027c58

Please sign in to comment.