import sys import os sys.path.insert(0, "/var/www/html/web/scripts/") import requests, h5py, urllib3, json, math from bs4 import BeautifulSoup import numpy as np from ftplib import FTP from datetime import datetime, timedelta # Vypnutie varovaní urllib3.disable_warnings(urllib3.exceptions.InsecureRequestWarning) # --- PÔVODNÁ KONFIGURÁCIA (ZACHOVANÁ) --- BASE_PATH = "/tmp/rhi_work" LOCAL_DIR = os.path.join(BASE_PATH, "downloads_volume") OUTPUT_GEOJSON = f"{BASE_PATH}/rhi_3d_latest.json" FTP_HOST, FTP_USER, FTP_PASS = "ftp.meteopocasie.sk", "ncikala.meteopocasie.sk", "Supercela150" FTP_BASE_DIR = "web/meteoradar/rhi" # --- NOVINKA: CESTA PRE TVOJ ARCHÍV NA VPS --- ARCHIVE_DIR = "/var/www/html/web/modely/ecmwf/europe/" RADAR_LAT, RADAR_LON = 48.2556, 17.1524 MAX_DIST, MAX_ALT = 160, 22 DIST_STEP_KM = 1.5 ALT_STEP_KM = 0.5 WALL_THICKNESS_KM = 1.0 # Vytvorenie lokálnych priečinkov os.makedirs(LOCAL_DIR, exist_ok=True) os.makedirs(ARCHIVE_DIR, exist_ok=True) os.chmod(BASE_PATH, 0o777) def get_latest_volume_url(): base_url = "https://opendata.shmu.sk/meteorology/weather/radar/volume/skjav/dBZ/" try: r = requests.get(base_url, verify=False, timeout=15) soup = BeautifulSoup(r.text, 'html.parser') dirs = sorted([a['href'] for a in soup.find_all('a') if a['href'].startswith('202')]) if not dirs: return None latest_dir = base_url + dirs[-1] r = requests.get(latest_dir, verify=False, timeout=15) soup = BeautifulSoup(r.text, 'html.parser') files = sorted([a['href'] for a in soup.find_all('a') if a['href'].endswith('.hdf')]) return latest_dir + files[-1] if files else None except: return None def find_storm_azimuth(hdf_path): with h5py.File(hdf_path, 'r') as f: data = f['dataset1/data1/data'][:] valid_data = np.where(data == 255, 0, data) idx = np.unravel_index(np.argmax(valid_data, axis=None), valid_data.shape) return int(idx[0]) def calculate_end_point(lat, lon, az, dist_km): R = 6371.0 brng = np.radians(az) d = dist_km lat1, lon1 = np.radians(lat), np.radians(lon) lat2 = np.arcsin(np.sin(lat1)*np.cos(d/R) + np.cos(lat1)*np.sin(d/R)*np.cos(brng)) lon2 = lon1 + np.arctan2(np.sin(brng)*np.sin(d/R)*np.cos(lat1), np.cos(d/R)-np.sin(lat1)*np.sin(lat2)) return [np.degrees(lon2), np.degrees(lat2)] def get_oriented_polygon(center_lon, center_lat, azimuth_deg, length_km, width_km): lat_to_km = 111.32 lon_to_km = 111.32 * math.cos(math.radians(center_lat)) az_rad = math.radians(azimuth_deg) u_x, u_y = math.sin(az_rad), math.cos(az_rad) v_x, v_y = math.cos(az_rad), -math.sin(az_rad) L, W = length_km / 2.0, width_km / 2.0 corners = [ ( L*u_x - W*v_x, L*u_y - W*v_y), ( L*u_x + W*v_x, L*u_y + W*v_y), (-L*u_x + W*v_x, -L*u_y + W*v_y), (-L*u_x - W*v_x, -L*u_y - W*v_y) ] poly_coords = [] for dx, dy in corners: plon = center_lon + (dx / lon_to_km) plat = center_lat + (dy / lat_to_km) poly_coords.append([round(plon, 5), round(plat, 5)]) poly_coords.append(poly_coords[0]) return [poly_coords] def get_color_for_dbz(val): if val < 15: return "#00a08a" if val < 25: return "#00ff00" if val < 35: return "#ffff00" if val < 45: return "#ff8000" if val < 55: return "#ff0000" return "#ff00ff" def generate_rhi_3d(hdf_path, azimuth): with h5py.File(hdf_path, 'r') as f: raw_time = f['what'].attrs.get('time').decode('utf-8') raw_date = f['what'].attrs.get('date').decode('utf-8') dt_local = datetime.strptime(f"{raw_date}{raw_time}", "%Y%m%d%H%M%S") + timedelta(hours=1) datasets = sorted([k for k in f.keys() if k.startswith('dataset')], key=lambda x: int(x.replace('dataset', ''))) dist_grid = np.arange(0, MAX_DIST, DIST_STEP_KM) alt_grid = np.arange(0, MAX_ALT, ALT_STEP_KM) rhi_matrix = np.full((len(alt_grid), len(dist_grid)), np.nan) for ds_name in datasets: ds = f[ds_name] el = ds['where'].attrs.get('elangle') gain = ds['data1/what'].attrs.get('gain', 0.5) offset = ds['data1/what'].attrs.get('offset', -32.0) nodata = ds['data1/what'].attrs.get('nodata', 255) rscale = ds['where'].attrs.get('rscale', 1000) / 1000.0 raw_vals = ds['data1/data'][int(azimuth % 360), :].astype(float) dbz_vals = np.where(raw_vals != nodata, (raw_vals * gain) + offset, np.nan) for i, d in enumerate(dist_grid): Re = 6371.0 * (4/3) h_min = np.sqrt(d**2 + Re**2 + 2*d*Re*np.sin(np.radians(el - 0.5))) - Re h_max = np.sqrt(d**2 + Re**2 + 2*d*Re*np.sin(np.radians(el + 0.5))) - Re bin_idx = int(d / rscale) if bin_idx < len(dbz_vals) and not np.isnan(dbz_vals[bin_idx]): val = dbz_vals[bin_idx] if val > 8: alt_mask = (alt_grid >= h_min) & (alt_grid <= h_max) rhi_matrix[alt_mask, i] = np.fmax(rhi_matrix[alt_mask, i], val) features = [] for i, d in enumerate(dist_grid): if d == 0: continue center_lon, center_lat = calculate_end_point(RADAR_LAT, RADAR_LON, azimuth, d) polygon_coords = get_oriented_polygon(center_lon, center_lat, azimuth, DIST_STEP_KM, WALL_THICKNESS_KM) for j, h in enumerate(alt_grid): val = rhi_matrix[j, i] if not np.isnan(val) and val >= 10: features.append({ "type": "Feature", "geometry": {"type": "Polygon", "coordinates": polygon_coords}, "properties": { "dbz": float(val), "color": get_color_for_dbz(val), "base": int(h * 1000), "height": int((h + ALT_STEP_KM) * 1000) } }) geojson_data = { "type": "FeatureCollection", "features": features, "metadata": { "timestamp": dt_local.strftime("%H:%M"), "azimuth": azimuth, "total_blocks": len(features) } } # Uloženie lokálne do /tmp/ with open(OUTPUT_GEOJSON, 'w') as f: json.dump(geojson_data, f) # Uloženie kópie do archívu /modely/ecmwf/europe/ archive_path = os.path.join(ARCHIVE_DIR, f"rhi_3d_{dt_local.strftime('%H%M')}.json") with open(archive_path, 'w') as f: json.dump(geojson_data, f) return dt_local, raw_date + raw_time[:4] def upload_rhi_3d(dt, rtime, az): try: with FTP(FTP_HOST) as ftp: ftp.login(FTP_USER, FTP_PASS) ftp.cwd('/') for d in FTP_BASE_DIR.split('/'): if d: ftp.cwd(d) with open(OUTPUT_GEOJSON, 'rb') as f: ftp.storbinary("STOR rhi_3d_latest.json", f) return True except Exception as e: return str(e) if __name__ == "__main__": azimuth_to_use = None if len(sys.argv) > 1: try: azimuth_to_use = int(sys.argv[1]) except: pass url = get_latest_volume_url() if url: fname = url.split('/')[-1] hdf_p = os.path.join(LOCAL_DIR, fname) r = requests.get(url, stream=True, verify=False) with open(hdf_p, 'wb') as f: for c in r.iter_content(8192): f.write(c) if azimuth_to_use is None: azimuth_to_use = find_storm_azimuth(hdf_p) dt, rtime = generate_rhi_3d(hdf_p, azimuth_to_use) status = upload_rhi_3d(dt, rtime, azimuth_to_use) if os.path.exists(hdf_p): os.remove(hdf_p) # FINÁLNY VÝSTUP PRE PHP print(json.dumps({ "status": "success", "azimuth": azimuth_to_use, "ftp_upload": status, "archive_saved": True }))