Display a DEM as a 3D surface
Kaufmann, 2024.
[1]:
import os
os.environ["TRAME_DISABLE_V3_WARNING"] = "1" # disable trame v3 warning
import pyvista as pv
from geometron.utils.geoservices import esri_rest_server_geotiff
from geometron.vtk import geotiff_to_grid
Get the DEM from WalOnMap (not possible because these are images?)
[2]:
host = 'www.imagetest.geomatics.gov.nt.ca'
folder = 'GNWT_Basemaps'
service = 'Aster_DEM_Basemap'
service_type = 'MapServer'
[3]:
xmin, xmax, ymin, ymax = 620400.0, 620600.0, 626700.0, 627000.0 # map extent in epsg:31370
size_width, size_height = 768, 768
epsg = 3857
[4]:
# Parameters
params = {'bbox': f'{xmin},{ymin},{xmax},{ymax}',
'size': f'{size_width}, {size_height}',
'bboxSR':{"wkid" : epsg},
'imageSR':{"wkid" : epsg},
'format': 'image/geotiff',
'transparent': 'false',
'layerDefs': '',
'layers' : '',
'dpi' : '100',
'time':'',
'f':'image'
}
[5]:
image_file='./data/dem.tif'
[6]:
esri_rest_server_geotiff(image_file=image_file, host=host,folder=folder, service=service, service_type=service_type, parameter=params, verbose=True)
url:
https://www.imagetest.geomatics.gov.nt.ca/arcgis/rest/services/GNWT_Basemaps/Aster_DEM_Basemap/MapServer/export?bbox=620400.0%2C626700.0%2C620600.0%2C627000.0&size=768%2C+768&bboxSR=%7B%27wkid%27%3A+3857%7D&imageSR=%7B%27wkid%27%3A+3857%7D&format=image%2Fgeotiff&transparent=false&layerDefs=&layers=&dpi=100&time=&f=pjson
url:
https://www.imagetest.geomatics.gov.nt.ca/arcgis/rest/services/GNWT_Basemaps/Aster_DEM_Basemap/MapServer/export?bbox=620400.0%2C626700.0%2C620600.0%2C627000.0&size=768%2C+768&bboxSR=%7B%27wkid%27%3A+3857%7D&imageSR=%7B%27wkid%27%3A+3857%7D&format=image%2Fgeotiff&transparent=false&layerDefs=&layers=&dpi=100&time=&f=image
/home/su530201/.local/share/virtualenvs/geometron-VhFx4lLH/lib/python3.10/site-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
Load DEM as a grid
[7]:
dem = geotiff_to_grid('./data/dem.tif', warp=True, factor=5.)
ERROR:root:Cannot read tile : 0,0 from file
/home/su530201/.local/share/virtualenvs/geometron-VhFx4lLH/lib/python3.10/site-packages/pyvista/core/utilities/fileio.py:203: UserWarning: The VTK reader `vtkTIFFReader` in pyvista reader `TIFFReader('/home/su530201/PycharmProjects/geometron/docs/source/examples/geological/data/dem.tif')` raised an errorwhile reading the file.
"Cannot read tile : 0,0 from file"
warnings.warn(
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[7], line 1
----> 1 dem = geotiff_to_grid('./data/dem.tif', warp=True, factor=5.)
File ~/PycharmProjects/geometron/geometron/vtk/image_to_vtk.py:25, in geotiff_to_grid(filename, warp, scalars_name, factor)
23 transform = shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src.read_transform()), force_3d=True)
24 grid = grid.transform(transform)
---> 25 no_data = [n for n, i in enumerate(grid.point_data['Tiff Scalars']) if i in src.nodatavals]
26 grid.point_data['Tiff Scalars'][no_data] = np.nan
27 grid.set_active_scalars('Tiff Scalars')
File ~/PycharmProjects/geometron/geometron/vtk/image_to_vtk.py:25, in <listcomp>(.0)
23 transform = shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src.read_transform()), force_3d=True)
24 grid = grid.transform(transform)
---> 25 no_data = [n for n, i in enumerate(grid.point_data['Tiff Scalars']) if i in src.nodatavals]
26 grid.point_data['Tiff Scalars'][no_data] = np.nan
27 grid.set_active_scalars('Tiff Scalars')
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
[8]:
import rasterio as rio
from geometron.geometries.transforms import shapely_matrix_to_vtk_matrix, rasterio_matrix_to_shapely_matrix
[9]:
filename = './data/dem.tif'
[10]:
with rio.open(filename, mode='r') as src:
img = src.read()
img
[10]:
array([[[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
...,
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253]],
[[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
...,
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253]],
[[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
...,
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253],
[253, 253, 253, ..., 253, 253, 253]],
[[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255],
...,
[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255],
[255, 255, 255, ..., 255, 255, 255]]], dtype=uint8)
[11]:
grid = img.cast_to_structured_grid()
src = rio.open(filename)
transform = shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src.read_transform()), force_3d=True)
grid = grid.transform(transform)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[11], line 1
----> 1 grid = img.cast_to_structured_grid()
2 src = rio.open(filename)
3 transform = shapely_matrix_to_vtk_matrix(rasterio_matrix_to_shapely_matrix(src.read_transform()), force_3d=True)
AttributeError: 'numpy.ndarray' object has no attribute 'cast_to_structured_grid'
[12]:
no_data = [n for n, i in enumerate(grid.point_data['Tiff Scalars']) if i in src.nodatavals]
grid.point_data['Tiff Scalars'][no_data] = np.nan
grid.set_active_scalars('Tiff Scalars')
if warp:
grid = grid.warp_by_scalar('Tiff Scalars', factor=factor)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 no_data = [n for n, i in enumerate(grid.point_data['Tiff Scalars']) if i in src.nodatavals]
2 grid.point_data['Tiff Scalars'][no_data] = np.nan
3 grid.set_active_scalars('Tiff Scalars')
NameError: name 'grid' is not defined
[13]:
grid.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 grid.plot()
NameError: name 'grid' is not defined
[14]:
grid
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 grid
NameError: name 'grid' is not defined
[15]:
grid.set_active_scalars('Tiff Scalars')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 grid.set_active_scalars('Tiff Scalars')
NameError: name 'grid' is not defined
[16]:
grid.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 grid.plot()
NameError: name 'grid' is not defined
[ ]: