import warnings
warnings.filterwarnings('ignore')

import pyvista as pv
import rioxarray as riox
import numpy as np

# Read the data from a DEM file

path_file = r'C:/Users/m.burhanuddin/OneDrive - dservice GmbH/Desktop/PRESENTATION/MATERI SURVEY BURHAN\DEMNAS/DEMNAS_2516-34_v1.0.tif'
data = riox.open_rasterio(path_file)
data = data[0]

# Save the raster data as an array
values = np.asarray(data)
 
# Create a mesh grid
x, y = np.meshgrid(data['x'], data['y'])
 
# Set the z values and create a StructuredGrid
z = np.zeros_like(x)
mesh = pv.StructuredGrid(x, y, z)
 
# Assign Elevation Values
mesh["Elevation"] = values.ravel(order='F')
 
# Warp the mesh by scalar
topo = mesh.warp_by_scalar(scalars="Elevation", factor=0.0000015)

# Plot the elevation map
p = pv.Plotter()
p.add_mesh(mesh=topo, scalars=topo["Elevation"], cmap='terrain')
p.show_grid(color='black')
p.set_background(color='white')
p.show(cpos="xy")

import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt

# Path to your DEM TIFF file
dem_file = path_file 

# Open the DEM file using rasterio
with rasterio.open(dem_file) as dem:
    # Read the DEM data
    dem_data = dem.read(1)  # Read the first band
    
    # Plot the DEM data
    plt.figure(figsize=(10, 6))
    show(dem, title='Digital Elevation Model (DEM)')
    plt.show()

import numpy as np

# Open the DEM file using rasterio
with rasterio.open(dem_file) as dem:
    # Read the DEM data
    dem_data = dem.read(1)  # Read the first band
    
    # Mask no-data values (usually indicated by a specific value, e.g., -9999)
    dem_data = np.ma.masked_equal(dem_data, dem.nodata)

    # Create a figure and axis
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # Plot the DEM data with a colormap and colorbar
    cax = ax.imshow(dem_data, cmap='terrain')
    cbar = fig.colorbar(cax, ax=ax, orientation='vertical')
    cbar.set_label('Elevation (meters)')
    
    # Set plot title and labels
    ax.set_title('Digital Elevation Model (DEM)')
    ax.set_xlabel('Longitude')
    ax.set_ylabel('Latitude')
    
    # Show the plot
    plt.show()

from osgeo import gdal
import numpy as np
import matplotlib 
import matplotlib.pyplot as plt
import elevation 

gdal_data = gdal.Open(path_file)
gdal_band = gdal_data.GetRasterBand(1)
nodataval = gdal_band.GetNoDataValue()

# convert to a numpy array
data_array = gdal_data.ReadAsArray().astype(np.float64)
data_array

# replace missing values if necessary
if np.any(data_array == nodataval):
    data_array[data_array == nodataval] = np.nan

#Plot out data with Matplotlib's 'contour'
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contour(data_array, cmap = "viridis", 
            levels = list(range(0, 5000, 100)))
plt.title("Elevation Contours of Selat Malaka")
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

#Plot our data with Matplotlib's 'contourf'
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contourf(data_array, cmap = "viridis", 
            levels = list(range(0, 5000, 500)))
plt.title("Elevation Contours of Selat Malaka")
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

import os
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# Open the DEM file
dem_file = r'C:/Users/m.burhanuddin/OneDrive - dservice GmbH/Desktop/PRESENTATION/MATERI SURVEY BURHAN\DEMNAS/BATNAS_100E-105E_000-05N_MSL_v1.1.tif'
dataset = gdal.Open(dem_file)

# Get the first band
band = dataset.GetRasterBand(1)

# Get NoData value
nodataval = band.GetNoDataValue()

# Read the data as an array and convert to float
data_array = band.ReadAsArray().astype(np.float64)  # Use np.float64 or float

# Mask the no-data values
data_array = np.ma.masked_equal(data_array, nodataval)

# Plot the DEM data
plt.figure(figsize=(10, 6))
cax = plt.imshow(data_array, cmap='terrain')
cbar = plt.colorbar(cax, orientation='vertical')
cbar.set_label('Elevation (meters)')
plt.title('Digital Elevation Model (DEM)')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()