11# std lib
2+ from typing import Optional
23from functools import partial
34import math
45
@@ -21,7 +22,6 @@ class cupy(object):
2122from xrspatial .utils import has_cuda
2223from xrspatial .utils import is_cupy_backed
2324
24- from typing import Optional
2525
2626def _run_numpy (data , azimuth = 225 , angle_altitude = 25 ):
2727 azimuth = 360.0 - azimuth
@@ -47,14 +47,23 @@ def _run_dask_numpy(data, azimuth, angle_altitude):
4747 meta = np .array (()))
4848 return out
4949
50+
5051def _run_dask_cupy (data , azimuth , angle_altitude ):
51- msg = 'Upstream bug in dask prevents cupy backed arrays '
52+ msg = 'Not implemented. '
5253 raise NotImplementedError (msg )
5354
55+
5456@cuda .jit
55- def _gpu_calc_numba (data , output , sin_altituderad , cos_altituderad , azimuthrad ):
57+ def _gpu_calc_numba (
58+ data ,
59+ output ,
60+ sin_altituderad ,
61+ cos_altituderad ,
62+ azimuthrad
63+ ):
64+
5665 i , j = cuda .grid (2 )
57- if i > 0 and i < data .shape [0 ]- 1 and j > 0 and j < data .shape [1 ]- 1 :
66+ if i > 0 and i < data .shape [0 ]- 1 and j > 0 and j < data .shape [1 ] - 1 :
5867 x = (data [i + 1 , j ]- data [i - 1 , j ])/ 2
5968 y = (data [i , j + 1 ]- data [i , j - 1 ])/ 2
6069
@@ -64,21 +73,14 @@ def _gpu_calc_numba(data, output, sin_altituderad, cos_altituderad, azimuthrad):
6473
6574 sin_slope = math .sin (slope )
6675 sin_part = sin_altituderad * sin_slope
67-
76+
6877 cos_aspect = math .cos (aspect )
6978 cos_slope = math .cos (slope )
7079 cos_part = cos_altituderad * cos_slope * cos_aspect
7180
7281 res = sin_part + cos_part
7382 output [i , j ] = (res + 1 ) * 0.5
7483
75- def calc_dims (shape ):
76- threadsperblock = (32 ,32 )
77- blockspergrid = (
78- (shape [0 ] + (threadsperblock [0 ] - 1 )) // threadsperblock [0 ],
79- (shape [1 ] + (threadsperblock [1 ] - 1 )) // threadsperblock [1 ]
80- )
81- return blockspergrid , threadsperblock
8284
8385def _run_cupy (d_data , azimuth , angle_altitude ):
8486 # Precompute constant values shared between all threads
@@ -87,19 +89,22 @@ def _run_cupy(d_data, azimuth, angle_altitude):
8789 cos_altituderad = np .cos (altituderad )
8890 azimuthrad = (360.0 - azimuth ) * np .pi / 180.
8991
90- # Allocate output buffer and launch kernel with appropriate dimensions
92+ # Allocate output buffer and launch kernel with appropriate dimensions
9193 output = cupy .empty (d_data .shape , np .float32 )
92- griddim , blockdim = calc_dims (d_data .shape )
93- _gpu_calc_numba [griddim , blockdim ](d_data , output , sin_altituderad , cos_altituderad , azimuthrad )
94+ griddim , blockdim = cuda_args (d_data .shape )
95+ _gpu_calc_numba [griddim , blockdim ](
96+ d_data , output , sin_altituderad , cos_altituderad , azimuthrad
97+ )
9498
9599 # Fill borders with nans.
96- output [ 0 , :] = cupy .nan
100+ output [0 , :] = cupy .nan
97101 output [- 1 , :] = cupy .nan
98102 output [:, 0 ] = cupy .nan
99103 output [:, - 1 ] = cupy .nan
100104
101105 return output
102106
107+
103108def hillshade (agg : xr .DataArray ,
104109 azimuth : int = 225 ,
105110 angle_altitude : int = 25 ,
@@ -114,7 +119,7 @@ def hillshade(agg: xr.DataArray,
114119 agg : xarray.DataArray
115120 2D NumPy, CuPy, NumPy-backed Dask, or Cupy-backed Dask array
116121 of elevation values.
117- altitude : int, default=25
122+ angle_altitude : int, default=25
118123 Altitude angle of the sun specified in degrees.
119124 azimuth : int, default=225
120125 The angle between the north vector and the perpendicular
@@ -136,83 +141,38 @@ def hillshade(agg: xr.DataArray,
136141 --------
137142 .. plot::
138143 :include-source:
139-
140- import datashader as ds
141- import matplotlib.pyplot as plt
142- from xrspatial import generate_terrain, hillshade
143-
144- # Create Canvas
145- W = 500
146- H = 300
147- cvs = ds.Canvas(plot_width = W,
148- plot_height = H,
149- x_range = (-20e6, 20e6),
150- y_range = (-20e6, 20e6))
151-
152- # Generate Example Terrain
153- terrain_agg = generate_terrain(canvas = cvs)
154-
155- # Edit Attributes
156- terrain_agg = terrain_agg.assign_attrs(
157- {
158- 'Description': 'Example Terrain',
159- 'units': 'km',
160- 'Max Elevation': '4000',
161- }
162- )
163-
164- terrain_agg = terrain_agg.rename({'x': 'lon', 'y': 'lat'})
165- terrain_agg = terrain_agg.rename('Elevation')
166-
167- # Create Hillshade Aggregate Array
168- hillshade_agg = hillshade(agg = terrain_agg, name = 'Illumination')
169-
170- # Edit Attributes
171- hillshade_agg = hillshade_agg.assign_attrs({'Description': 'Example Hillshade',
172- 'units': ''})
173-
174- # Plot Terrain
175- terrain_agg.plot(cmap = 'terrain', aspect = 2, size = 4)
176- plt.title("Terrain")
177- plt.ylabel("latitude")
178- plt.xlabel("longitude")
179-
180- # Plot Terrain
181- hillshade_agg.plot(cmap = 'Greys', aspect = 2, size = 4)
182- plt.title("Hillshade")
183- plt.ylabel("latitude")
184- plt.xlabel("longitude")
144+ import numpy as np
145+ import xarray as xr
146+ from xrspatial import hillshade
147+
148+ data = np.array([
149+ [0., 0., 0., 0., 0.],
150+ [0., 1., 0., 2., 0.],
151+ [0., 0., 3., 0., 0.],
152+ [0., 0., 0., 0., 0.],
153+ [0., 0., 0., 0., 0.]
154+ ])
155+ n, m = data.shape
156+ raster = xr.DataArray(data, dims=['y', 'x'], name='raster')
157+ raster['y'] = np.arange(n)[::-1]
158+ raster['x'] = np.arange(m)
159+
160+ hillshade_agg = hillshade(raster)
185161
186162 .. sourcecode:: python
187163
188- >>> print(terrain_agg[200:203, 200:202])
189- <xarray.DataArray 'Elevation' (lat: 3, lon: 2)>
190- array([[1264.02249454, 1261.94748873],
191- [1285.37061171, 1282.48046696],
192- [1306.02305679, 1303.40657515]])
193- Coordinates:
194- * lon (lon) float64 -3.96e+06 -3.88e+06
195- * lat (lat) float64 6.733e+06 6.867e+06 7e+06
196- Attributes:
197- res: 1
198- Description: Example Terrain
199- units: km
200- Max Elevation: 4000
201-
202- >>> print(hillshade_agg[200:203, 200:202])
203- <xarray.DataArray 'Illumination' (lat: 3, lon: 2)>
204- array([[1264.02249454, 1261.94748873],
205- [1285.37061171, 1282.48046696],
206- [1306.02305679, 1303.40657515]])
164+ >>> print(hillshade_agg)
165+ <xarray.DataArray 'hillshade' (y: 5, x: 5)>
166+ array([[ nan, nan, nan, nan, nan],
167+ [ nan, 0.71130913, 0.44167341, 0.71130913, nan],
168+ [ nan, 0.95550163, 0.71130913, 0.52478473, nan],
169+ [ nan, 0.71130913, 0.88382559, 0.71130913, nan],
170+ [ nan, nan, nan, nan, nan]])
207171 Coordinates:
208- * lon (lon) float64 -3.96e+06 -3.88e+06
209- * lat (lat) float64 6.733e+06 6.867e+06 7e+06
210- Attributes:
211- res: 1
212- Description: Example Hillshade
213- units:
214- Max Elevation: 4000
172+ * y (y) int32 4 3 2 1 0
173+ * x (x) int32 0 1 2 3 4
215174 """
175+
216176 # numpy case
217177 if isinstance (agg .data , np .ndarray ):
218178 out = _run_numpy (agg .data , azimuth , angle_altitude )
0 commit comments