<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Plane-DEM intersections
First dated version
Step1:  Test case 1
The first test case is illustrated in the image below. We have a horizontal topographic surface, at a height of 0, with 100 x 100 cells with a cell size of 1. The geological plane dips 45° towards East. The source point for the plane is located at (0, 50, 50).
The locations of the expected intersection points are (50, *, 0).
First, a horizontal plane was created with Saga GIS and saved in pygsf/example_data/horiz_plane.asc.
Step2:  We read the data source with success. So we may unpack the result.
Step3:  Hmmm, there is no projection info. In fact, there shouldn't..
Step4:  A dictionary, as suspected. Try to see the content..
Step5:  A very horizontal surface, we agree..
Step6:  Given these data, we store them into a GeoArray
Step7:  There is a single band provided in the geoarray, and represented by the data array.
The signature of the plane-DEM intersection function is
Step8:  The source point is located at (0, 50, 50)
Step9:  Now we try calculating the intersection
Step10:  As expected, all the intersection points lie at (50, *, 0)
Plotting with Bokeh..
Step11:  Test case 2
Now we consider a horizontal plane at z = 0 as topographic surface (same as case 1) and another horizontal surface at z = 1 as geological plane. We should get no intersection.
Step12:  The horizontal geological plane definition
Step13:  The source point located at (0, 50, 1)
Step14:  Ok, list is empty, as expected.
Test case 3
Now we consider a horizontal plane at z = 0 as topographic surface (same as case 1) and another horizontal surface at z = 0 as geological plane. We should get all grid points as intersections.
The variables are the same as Case 2, apart from the point definition
Step15:  They seem correct, just quite numerous..
We visualize them with Bokeh.
<ASSISTANT_TASK:>
Python Code:
from pygsf.io.gdal.raster import try_read_raster_band
Explanation: Plane-DEM intersections
First dated version: 2019-06-11
Current version: 2021-04-24
Last run: 2021-04-24
A few simulated topographic surfaces were used to validate the routine for calculating the plane-DEM intersection.
Loading the dataset can be made with the following function:
End of explanation
source_data = "/home/mauro/Documents/projects/gsf/example_data/others/horiz_plane.asc"
success, cntnt = try_read_raster_band(raster_source=source_data)
print(success)
Explanation: Test case 1
The first test case is illustrated in the image below. We have a horizontal topographic surface, at a height of 0, with 100 x 100 cells with a cell size of 1. The geological plane dips 45° towards East. The source point for the plane is located at (0, 50, 50).
The locations of the expected intersection points are (50, *, 0).
First, a horizontal plane was created with Saga GIS and saved in pygsf/example_data/horiz_plane.asc.
End of explanation
geotransform, projection, band_params, data = cntnt
type(geotransform)
print(geotransform)
type(projection)
print(projection)
Explanation: We read the data source with success. So we may unpack the result.
End of explanation
type(band_params)
Explanation: Hmmm, there is no projection info. In fact, there shouldn't..
End of explanation
print(band_params)
Explanation: A dictionary, as suspected. Try to see the content..
End of explanation
type(data)
data.shape
data.min()
data.max()
Explanation: A very horizontal surface, we agree..
End of explanation
from pygsf.georeferenced.rasters import GeoArray
ga = GeoArray(inGeotransform=geotransform, inLevels=[data])
Explanation: Given these data, we store them into a GeoArray:
End of explanation
from pygsf.orientations.orientations import Plane
gplane = Plane(azim=90.0, dip_ang=45.0)
print(gplane)
Explanation: There is a single band provided in the geoarray, and represented by the data array.
The signature of the plane-DEM intersection function is:
plane_dem_intersection (srcPlaneAttitude: Plane, srcPt: Point, geo_array: GeoArray, level_ndx: int=0) 
-> Tuple[List[Point], List[Point]]:
We already have the geoarray, we need to define the source plane attitue and the source point.
The geoplane is East-dipping with a dip angle of 45°:
End of explanation
from pygsf.geometries.shapes.space3d import Point3D
pt = Point3D(0, 50, 50)
Explanation: The source point is located at (0, 50, 50)
End of explanation
from pygsf.georeferenced.rasters import plane_dem_intersection
inters_pts = plane_dem_intersection(
        srcPlaneAttitude=gplane,
        srcPt=pt,
        geo_array=ga)
print(inters_pts)
Explanation: Now we try calculating the intersection:
End of explanation
from bokeh.plotting import figure, output_notebook, show
x = list(map(lambda pt: pt.x, inters_pts))
y = list(map(lambda pt: pt.y, inters_pts))
output_notebook()
p = figure()
p.circle(x, y, size=2, color="navy", alpha=0.5)
show(p)
Explanation: As expected, all the intersection points lie at (50, *, 0)
Plotting with Bokeh..
End of explanation
source_data = "/home/mauro/Documents/projects/gsf/example_data/others/horiz_plane.asc"
success, cntnt = try_read_raster_band(raster_source=source_data)
print(success)
geotransform, projection, band_params, data = cntnt
ga = GeoArray(inGeotransform=geotransform, inLevels=[data])
Explanation: Test case 2
Now we consider a horizontal plane at z = 0 as topographic surface (same as case 1) and another horizontal surface at z = 1 as geological plane. We should get no intersection.
End of explanation
from pygsf.orientations.orientations import Plane
gplane = Plane(azim=90.0, dip_ang=0.0)
Explanation: The horizontal geological plane definition:
End of explanation
pt = Point3D(0, 50, 1)
inters_pts = plane_dem_intersection(
        srcPlaneAttitude=gplane,
        srcPt=pt,
        geo_array=ga)
print(inters_pts)
Explanation: The source point located at (0, 50, 1)
End of explanation
pt = Point3D(0, 50, 0)
inters_pts = plane_dem_intersection(
        srcPlaneAttitude=gplane,
        srcPt=pt,
        geo_array=ga)
print(inters_pts)
Explanation: Ok, list is empty, as expected.
Test case 3
Now we consider a horizontal plane at z = 0 as topographic surface (same as case 1) and another horizontal surface at z = 0 as geological plane. We should get all grid points as intersections.
The variables are the same as Case 2, apart from the point definition:
End of explanation
from bokeh.plotting import figure, output_notebook, show
x = list(map(lambda pt: pt.x, inters_pts))
y = list(map(lambda pt: pt.y, inters_pts))
output_notebook()
p = figure()
p.circle(x, y, size=2, color="navy", alpha=0.5)
show(p)
Explanation: They seem correct, just quite numerous..
We visualize them with Bokeh.
End of explanation
<END_TASK>