File:Pressure of tunguska event if 10mt 1.png

   dem_file = 'tunguska_30m_dem.tif'
   
   if not os.path.exists(dem_file):
       command = [
           'eio', 'clip',
           '-o', dem_file,
           '--bounds', str(bounds[0]), str(bounds[1]), str(bounds[2]), str(bounds[3])
       ]
       
       try:
           subprocess.run(command, check=True)
           print("Korkeusdata ladattu ja leikattu onnistuneesti.")
       except subprocess.CalledProcessError as e:
           print(f"Virhe: {e}")
           return

   # Lataa korkeusdata
   with rasterio.open(dem_file) as src:
       elevation = src.read(1)
       transform = src.transform

   # Coordinates of the Tunguska event
   explosion_coords = (101.909722, 60.903056)
   explosion_height = 7000   
   # Calculate distances for the Pythagorean formula
   lon, lat = np.meshgrid(
       np.linspace(bounds[0], bounds[2], elevation.shape[1]),
       np.linspace(bounds[3], bounds[1], elevation.shape[0])
   )
   radius = 6371000  # Radius of Earth in meters

   # Convert lat/lon to distances in meters
   d_lat = np.radians(lat - explosion_coords[1])
   d_lon = np.radians(lon - explosion_coords[0])
   a = np.sin(d_lat / 2)**2 + np.cos(np.radians(explosion_coords[1])) * np.cos(np.radians(lat)) * np.sin(d_lon / 2)**2
   c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
   distance_from_explosion = radius * c
   distance_from_explosion_alt=np.sqrt(distance_from_explosion*distance_from_explosion+explosion_height*explosion_height)
   # Paineen laskeminen
   base_energy = 1  # 1 psi at 42.6 km if 10 Mt
   distance_threshold = 42600  # 40 km in meters
   distance_threshold2 = 400000  # 400 km in meters
   heat_energy = base_energy / (distance_from_explosion_alt / distance_threshold)**2
   heat_energy[distance_from_explosion_alt > distance_threshold2] = 0  # Nollaa etäisyydet yli 65 km

   # Rajaa maksimimäärä 25 cal/cm²
   heat_energy = np.clip(heat_energy, 0, 1000)

   # Custom color map
   colors = [(0, "darkblue"), (0.25, "cyan"), (0.5, "yellow"), (0.75, "orange"), (1, "red")]
   cmap = mcolors.LinearSegmentedColormap.from_list("custom_hot", colors, N=256)

   # Calculate hillshade using LightSource
   ls = LightSource(azdeg=280, altdeg=45)
   hillshade_raster = ls.shade(elevation, cmap=plt.cm.gray, vert_exag=1, blend_mode='overlay')

   # Visualisointi Cartopy-kirjastolla
   fig, ax = plt.subplots(figsize=(15, 10), subplot_kw={'projection': ccrs.PlateCarree()})

   # Plot hillshade with transparency
   img_hs = ax.imshow(hillshade_raster, extent=(bounds[0], bounds[2], bounds[1], bounds[3]), origin='upper', transform=ccrs.PlateCarree(), alpha=0.6)

   # Plot heat energy raster with custom color map
   img_heat = ax.imshow(heat_energy, cmap=cmap, extent=(bounds[0], bounds[2], bounds[1], bounds[3]), origin='upper', transform=ccrs.PlateCarree(), alpha=0.4)

   plt.colorbar(img_heat, ax=ax, label='Pressure estimation 40km=1 PSI', extend='both')

   # Plot features: Tunguska and Vanavara
   tunguska_event = Circle((101.909722, 60.903056), 0.01, color='red', alpha=0.7, transform=ccrs.PlateCarree())
   vanavara = Rectangle((102.2848 - 0.005, 60.3457 - 0.005), 0.03, 0.03, color='blue', alpha=0.7, transform=ccrs.PlateCarree())

   ax.add_patch(tunguska_event)
   ax.add_patch(vanavara)

   # Add circles for distances
   circle_7km = Circle((101.909722, 60.903056), 7 / 111.32, color='black', alpha=0.7, fill=False, linestyle='-', transform=ccrs.PlateCarree())
   circle_15km = Circle((101.909722, 60.903056), 15 / 111.32, color='red', alpha=0.7, fill=False, linestyle='-', transform=ccrs.PlateCarree())
   circle_35km = Circle((101.909722, 60.903056), 35 / 111.32, color='black', alpha=0.7, fill=False, linestyle='--', transform=ccrs.PlateCarree())

   #ax.add_patch(circle_7km)
   #ax.add_patch(circle_15km)
   #ax.add_patch(circle_35km)

   # Load and plot hydro data
   hydro_data = gpd.read_file('./hydro/HydroRIVERS_v10_si.shp')
   hydro_data = hydro_data.to_crs(ccrs.PlateCarree().proj4_init)  # Reproject to PlateCarree
   hydro_data.plot(ax=ax, color='blue', linewidth=1)

   # Aseta akselien rajat
   ax.set_xlim(bounds[0], bounds[2])
   ax.set_ylim(bounds[1], bounds[3])

   # Aseta suhteellinen laajuus
   ax.set_aspect('equal')

   # Add labels
   ax.text(101.909722+0.5, 60.903056+0.04, 'Tunguska event', transform=ccrs.PlateCarree(),
           horizontalalignment='right', fontsize=32, color='black')
   ax.text(102.2848, 60.3457, 'Vanavara', transform=ccrs.PlateCarree(),
           horizontalalignment='right', fontsize=32, color='black')

   ax.set_title('Pressure in Tunguska, PSI if 10 Mt', fontsize=16)

   # Kontuurit lämpöenergiasta
   energy_levels = [0.03, 0.1, 0.2,0.5,1,2,3,5,10,20,100,300,1000]  # Määrittele kontuurit
   contours = plt.contour(lon, lat, heat_energy, levels=energy_levels, colors='black', linewidths=1, alpha=0.7)

   # Lisää kontuurien labelit
   plt.clabel(contours, inline=True, fontsize=10, fmt='%.1f')

   plt.show()

1. Määritä koordinaatit (lon_min, lat_min, lon_max, lat_max) ja räjähdyksen korkeus

bounds = (101.5, 60.2, 102.5, 61.4) explosion_height = 7000 # Korkeus metrinä load_and_process_elevation(bounds, explosion_height)