Konvertieren Sie NETCDF-Dateien in TIF-Dateien ⇐ Python

[Anonymous]

Ich erstelle ein Python-Programm, das NetCDF-Dateien in GeoTIFFs konvertiert, sodass der Benutzer die Variable und den Zeitschritt über eine Tkinter-Schnittstelle auswählen kann.
Die Konvertierung funktioniert, aber ich habe ein Visualisierungsproblem:
Alle Werte gleich 0,0 werden in transparente Pixel umgewandelt, was ich möchte – aber auch kleine Werte wie 0,00001 oder 0,1 verschwinden (als transparent oder behandelt). NoData).
Ich möchte, dass nur 0,0-Werte transparent sind (NaN), und jeder positive oder negative Wert ungleich Null, auch sehr klein (z. B. 0,00001), sollte sichtbar bleiben.
Hier ist der Kernteil meines Codes:

Code: Select all

import tkinter as tk
from tkinter import filedialog, messagebox
import xarray as xr
import numpy as np
import rasterio
from rasterio.transform import from_origin
from scipy.interpolate import griddata
import pandas as pd

# ---------------- Conversion: real values ---------------- #
from shapely.geometry import Polygon, Point
from shapely.ops import unary_union

def convert_nc_to_tiff(nc_file, variable, timesteps, output_file_pattern, resolution=1.0):
"""Convert selected timesteps of a NetCDF variable to GeoTIFF(s) with float values,
clipped to the exact mesh domain (union of face polygons)."""

ds = xr.open_dataset(nc_file)

if variable not in ds.variables:
raise ValueError(f"Variable {variable} not found in file")

# Build mesh domain polygon
if "Mesh2d_face_x_bnd" in ds and "Mesh2d_face_y_bnd" in ds:
xb = ds["Mesh2d_face_x_bnd"].values
yb = ds["Mesh2d_face_y_bnd"].values

polygons = []
for i in range(xb.shape[0]):
coords = [(xb[i, j], yb[i, j]) for j in range(xb.shape[1]) if not np.isnan(xb[i, j])]
if len(coords) >= 3:
polygons.append(Polygon(coords))

domain = unary_union(polygons)

# Use centroids for interpolation points
x = np.nanmean(xb, axis=1)
y = np.nanmean(yb, axis=1)

elif "Mesh2d_face_x" in ds and "Mesh2d_face_y"  in ds:
x = ds["Mesh2d_face_x"].values
y = ds["Mesh2d_face_y"].values
domain = Polygon(np.column_stack((x, y))).convex_hull
else:
raise ValueError("Mesh face coordinates not found in NetCDF")

var = ds[variable]

# Grid extents
x_min, x_max = x.min(), x.max()
y_min, y_max = y.min(), y.max()

width = int(np.ceil((x_max - x_min) / resolution)) + 1
height = int(np.ceil((y_max - y_min) / resolution)) + 1

xi = np.linspace(x_min, x_max, width)
yi = np.linspace(y_min, y_max, height)
grid_x, grid_y = np.meshgrid(xi, yi)

# Precompute mask (pixels inside mesh domain)
mask_in = np.zeros(grid_x.shape, dtype=bool)
for i in range(grid_x.shape[0]):
for j in range(grid_x.shape[1]):
if domain.contains(Point(grid_x[i, j], grid_y[i, j])):
mask_in[i, j] = True

for t in timesteps:
data = var.isel(time=t).values
data = np.where(data == 0, np.nan, data)

# Interpolation
points = np.column_stack((x, y))
grid_data = griddata(points, data.flatten(), (grid_x, grid_y), method="linear")

# Fill gaps inside mesh
mask_nan = np.isnan(grid_data)
if np.any(mask_nan):
grid_data[mask_nan] = griddata(points, data.flatten(), (grid_x, grid_y), method="nearest")[mask_nan]

# Apply domain mask
grid_data[~mask_in] = np.nan

# Flip Y so north is up
grid_data = np.flipud(grid_data)

transform = from_origin(x_min, y_max, resolution, resolution)
out_path = output_file_pattern.format(timestep=t)

with rasterio.open(
out_path,
"w",
driver="GTiff",
height=grid_data.shape[0],
width=grid_data.shape[1],
count=1,
dtype="float32",
crs=None,
transform=transform,
nodata=np.nan,
) as dst:
dst.write(grid_data.astype("float32"), 1)

ds.close()

# ---------------- GUI ---------------- #
class NC2TIFFApp:
def __init__(self, root):
self.root = root
self.root.title("NetCDF to GeoTIFF Converter (Values)")
self.root.geometry("600x700")  # visible window size

print("NC2TIFF GUI launched")  # debug

self.nc_file = None
self.ds = None

tk.Button(root, text="Select .nc File", command=self.load_file).pack(pady=5)

tk.Label(root, text="Select variable:").pack()
self.var_listbox = tk.Listbox(root, width=50, height=10, exportselection=False)
self.var_listbox.pack(pady=5)

tk.Label(root, text="Select time:").pack()
self.time_listbox = tk.Listbox(
root, width=50, height=10, selectmode=tk.MULTIPLE, exportselection=False
)
self.time_listbox.pack(pady=5)

tk.Label(root, text="Output file pattern:").pack()
self.output_entry = tk.Entry(root, width=50)
self.output_entry.insert(0, "output_t{timestep}.tif")
self.output_entry.pack(pady=5)

tk.Button(root, text="Convert to GeoTIFF(s)", command=self.convert).pack(pady=10)

def load_file(self):
file = filedialog.askopenfilename(filetypes=[("NetCDF files", "*.nc")])
if file:
if self.ds is not None:
self.ds.close()
self.nc_file = file
self.ds = xr.open_dataset(file)

# Variables
self.var_listbox.delete(0, tk.END)
for v in self.ds.variables:
self.var_listbox.insert(tk.END, v)

# Timesteps (works even if 'time' is coordinate)
self.time_listbox.delete(0, tk.END)
if "time"  in self.ds:
for i, t in enumerate(self.ds["time"].values):
self.time_listbox.insert(tk.END, f"{i}: {pd.to_datetime(str(t))}")
else:
self.time_listbox.insert(tk.END, "0")

messagebox.showinfo("File Loaded", f"Loaded {file} and populated variables/timesteps.")

def convert(self):
if not self.nc_file:
messagebox.showerror("Error", "Please select a .nc file first")
return

sel_var = self.var_listbox.curselection()
if not sel_var:
messagebox.showerror("Error", "Please select a variable")
return
variable = self.var_listbox.get(sel_var[0])

sel_times = self.time_listbox.curselection()
if not sel_times:
messagebox.showerror("Error", "Please select at least one timestep")
return

pattern = self.output_entry.get().strip()
if not pattern:
messagebox.showerror("Error", "Please specify output file pattern")
return

try:
convert_nc_to_tiff(self.nc_file, variable, sel_times, pattern)
messagebox.showinfo("Success", f"Saved {len(sel_times)} GeoTIFF file(s).")
except Exception as e:
messagebox.showerror("Error", str(e))

if __name__ == "__main__":
root = tk.Tk()
app = NC2TIFFApp(root)
root.mainloop()

Wenn ich jedoch die Ausgabe in QGIS überprüfe, fehlen auch Pixel mit Werten wie 0,0001 oder 0,000001 – es scheint, dass sie irgendwo im Prozess als NaN behandelt werden.
Was ich versucht habe:

• Mit np.isclose(data, 0.0, atol=1e-12) → dasselbe Problem.
• Einstellung von nodata=None → alle Pixel werden undurchsichtig (auch außerhalb des Netzes).
• Die Überprüfung der Daten vor der Interpolation zeigt, dass kleine Werte vorhanden sind.

Erwartetes Verhalten:

• 0,0 → transparent (NaN)
• 0,0 oder < 0,0 → sichtbar (farbig in QGIS)

Umgebung:

• Python 3.10
• Bibliotheken: xarray, rasterio, scipy, shapely, tkinter