Wyodrębnij dane z NetCDF za pomocą wielokąta kształtu pliku w Pythonie [duplikat]

Muszę podgrupować dane z NetCDF przy użyciu określonego pliku shapefile. Dane dotyczą temperatury powierzchni morza i koloru oceanu przy rozdzielczości 1/4 stopnia. Mam 4 wielokąty opisujące USA. Duży ekosystem morski północno-wschodniego szelfu kontynentalnego i jego podskładniki, których muszę użyć do wyodrębnienia danych. Pracuję z miesięcznymi plikami złożonymi z lat 1982-2014, więc ta procedura ekstrakcji danych musi zostać zautomatyzowana. Pliki są już podzielone na przybliżoną siatkę obszaru roboczego wynoszącą [35, 45, -80, -60].

Wcześniej konwertujemy pliki danych HDF5 na rastry w R i przetwarzamy je w ten sposób, ale ta metoda jest naprawdę nieefektywna i jestem pewien, że istnieje lepsze rozwiązanie w Pythonie przy użyciu bieżących plików NetCDF.

Do tej pory korzystałem z GDAL i Fiony do odczytu plików shapefiles i NetCDF4 do ładowania plików danych. Nie jestem pewien, jak przejść do podzbioru danych. Znalazłem to:

GDAL dla Pythona: wyodrębnianie poddomen z pliku NetCDF?

Ale nie mam mglistego pomysłu na temat tego, jak podgrupować plik NetCDF przy użyciu czegokolwiek innego niż prosta obwiednia, którym z pewnością nie są te wielokąty.

Punkt w procedurach wielokąta prawdopodobnie zajmie wieczność do pracy, ale może mógłbym podgrupować dane za pomocą mniejszej ramki granicznej, która jest obracana, aby dopasować te kształty jako początkowy punkt początkowy, a następnie przeprowadzić wyszukiwanie punktowe:

Podzbiór krzywoliniowego pliku netCDF (wyjście modelu ROMS) za pomocą ramki ograniczającej lon / lat.

Jakieś pomysły?

EDYCJA 1:

Właśnie natrafiłem na pakiet OpenClimateGIS, który wydaje się, że może idealnie pasować do rachunku ... Spróbuję to sprawdzić, czy uda mi się go uruchomić: http://ncpp.github.io/ocgis/examples. HTML # Advanced-Subsetting

Może to być dostosowane do twoich potrzeb.

Jeśli nie przeszkadza ci wywoływanie linii poleceń z Pythona, możesz zrobić coś takiego gdalwarp -cutline clip.shp -cl clip -crop_to_cutline input_raster output_raster_clipped.tif. -cwherei -csqlmogą być bardziej odpowiednie opcje gdalwarp do wybierania jednego z czterech wielokątów do obcinania.

Spójrz na to: /programming/34585582/how-to-mask-the-specific-array-data-based-on-the-shapefile

Pamiętaj, że po załadowaniu NetCDF pracujesz z tablicą NumPy.

Co próbujesz wyprodukować? Statystyki podsumowujące oparte na obszarach wielokąta?

Tak czy inaczej, oto co bym zrobił:

1. Załaduj plik shapefile i ustaw obszary w obsługiwanym formacie (celowanie w proces maskowania matplotlib w powyższym linku brzmi dobrze)
2. Załaduj plik NetCDF i umieść dane w jednej tablicy X, Y, T
3. Zamaskuj tę tablicę za pomocą wielokątów (pojedynczo?)
4. Eksportuj statystyki podsumowujące.

Możesz użyć rioxarray. Oto przykład: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray
import geopandas

geodf = geopandas.read_file(...)
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)