Конвертируйте файлы NETCDF в файлы TIF ⇐ Python

[Anonymous]

Я создаю программу на Python, которая преобразует файлы NetCDF в файлы GeoTIFF, позволяя пользователю выбирать переменную и временной шаг с помощью интерфейса Tkinter.
Преобразование работает, но у меня возникает проблема с визуализацией:
все значения, равные 0,0, превращаются в прозрачные пиксели, чего я и хочу, но небольшие значения, такие как 0,00001 или 0,1, также исчезают (расцениваются как прозрачные или NoData).
Я хочу, чтобы значения строго 0,0 были прозрачными (NaN), и любое положительное или отрицательное ненулевое значение, даже очень маленькое (например, 0,00001), должно оставаться видимым.
Вот основная часть моего кода:

Код: Выделить всё

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()

Однако, когда я проверяю выходные данные в QGIS, пиксели со значениями вроде 0,0001 или 0,000001 также отсутствуют — кажется, где-то в процессе они обрабатываются как NaN.
Что я пробовал:

• Использование np.isclose(data, 0.0, atol=1e-12) → та же проблема.
• Установка nodata=None → все пиксели становятся непрозрачными (даже за пределами сетки).
• Проверка данных перед интерполяцией показывает, что небольшие значения действительно существуют.

Ожидаемое поведение:

• 0.0 → прозрачный (NaN)
• 0,0 или < 0,0 → видимый (цветной в QGIS)

Среда:

• Python 3.10
• Библиотеки: xarray, rasterio, scipy, shapely, tkinter

Подробнее здесь: https://stackoverflow.com/questions/797 ... -tif-files