Przejście z listy mailingowej gdal-dev:
W poniedziałek, 2 września 2013 o 19:09, David Shean napisał:
Cześć lista, próbuję spakować rastry GTiff o identycznej projekcji / zasięgu / rozdzielczości jak pojedynczy plik NetCDF do dystrybucji. Spędziłem ostatnią godzinę konsultując się z internetowym dokumentem i grając bez żadnych sukcesów z gdal_translate, gdalbuildvrt i gdalwarp.
Czy istnieje prosty sposób, aby to zrobić za pomocą istniejących narzędzi wiersza poleceń gdal? Pomyślałem, że zapytam przed skorzystaniem z niestandardowego rozwiązania korzystającego z interfejsu API NetCDF Python.
Dzięki. -David
W wtorek, 3 września 2013 o 10:15, Etienne Tourigny napisała:
to, czego chcesz, prawdopodobnie wykracza poza zakres gdal. Wymagałoby to sprytnego zarządzania metadanymi, aby gdal_translate umieścił je w jednym pliku ...
Radzę przekonwertować je wszystkie na netcdf za pomocą gdal_translate, a następnie użyć python-netcdf4 (nie ten z numpy / scipy), aby ułożyć je w stos czasowy.
W wtorek, 3 września 2013 r., O godzinie 07:55, „Signell, Richard” napisał:
David, jeśli opublikujesz swoje pytanie w grupie GIS stackexchange /gis//, podam przykładowy kod, który powinien być pomocny.
-Bogaty
====================
Aktualizacja 9/3/13 17:04 PDT
Oto dane wyjściowe gdalinfo dla jednego z moich wejściowych zestawów danych:
gdalinfo 20120901T2024_align_x+22.19_y+3.68_z+14.97_warp.tif
Driver: GTiff/GeoTIFF
Files: 20120901T2024_align_x+22.19_y+3.68_z+14.97_warp.tif
Size is 10666, 13387
Coordinate System is:
PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",70],
PARAMETER["central_meridian",-45],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (-211346.063781524338992,-2245136.291794800199568)
Pixel Size = (5.000000000000000,-5.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -211346.064,-2245136.292) ( 50d22'39.70"W, 69d23'55.59"N)
Lower Left ( -211346.064,-2312071.292) ( 50d13'22.38"W, 68d48'10.75"N)
Upper Right ( -158016.064,-2245136.292) ( 49d 1'33.33"W, 69d26'16.42"N)
Lower Right ( -158016.064,-2312071.292) ( 48d54'35.06"W, 68d50'27.28"N)
Center ( -184681.064,-2278603.792) ( 49d38' 1.32"W, 69d 7'17.04"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
NoData Value=-32767
Kontynuacja sugerowanego podejścia Luke'a.
Generacja VRR działa dobrze:
gdalbuildvrt -separate newtest.vrt *warp.tif
<VRTDataset rasterXSize="10666" rasterYSize="13387">
<SRS>PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]</SRS>
<GeoTransform> -2.1134606378152434e+05, 5.0000000000000000e+00, 0.0000000000000000e+00, -2.2451362917948002e+06, 0.0000000000000000e+00, -5.0000000000000000e+00</GeoTransform>
<VRTRasterBand dataType="Float32" band="1">
<NoDataValue>-3.27670000000000E+04</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="1">20110619T2024_align_x+15.51_y+1.15_z+12.10_warp.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="10666" RasterYSize="13387" DataType="Float32" BlockXSize="256" BlockYSize="256" />
<SrcRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<DstRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<NODATA>-32767</NODATA>
</ComplexSource>
</VRTRasterBand>
<VRTRasterBand dataType="Float32" band="2">
<NoDataValue>-3.27670000000000E+04</NoDataValue>
<ComplexSource>
<SourceFilename relativeToVRT="1">20110802T2024_align_x+16.33_y+2.14_z+12.02_warp.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="10666" RasterYSize="13387" DataType="Float32" BlockXSize="256" BlockYSize="256" />
<SrcRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<DstRect xOff="0" yOff="0" xSize="10666" ySize="13387" />
<NODATA>-32767</NODATA>
</ComplexSource>
</VRTRasterBand>
...
Ale kiedy próbuję przetłumaczyć na nc, pojawia się następujący błąd:
gdal_translate -of netcdf newtest.vrt newtest.nc
Input file size is 10666, 13387
Warning 1: Variable has 0 dimension(s) - not supported.
0...10...20...30...40...50ERROR 1: netcdf error #-62 : NetCDF: One or more variable sizes violate format constraints .
at (netcdfdataset.cpp,SetDefineMode,1574)
ERROR 1: netcdf error #-39 : NetCDF: Operation not allowed in define mode .
at (netcdfdataset.cpp,IWriteBlock,1435)
ERROR 1: netCDF scanline write failed: NetCDF: Operation not allowed in define mode
ERROR 1: An error occured while writing a dirty block
...ERROR 1: netcdf error #-39 : NetCDF: Operation not allowed in define mode .
at (netcdfdataset.cpp,IWriteBlock,1435)
ERROR 1: netCDF scanline write failed: NetCDF: Operation not allowed in define mode
ERROR 1: netcdf error #-62 : NetCDF: One or more variable sizes violate format constraints .
at (netcdfdataset.cpp,~netCDFDataset,1548)
Po bliższym przyjrzeniu się wydaje się, że gdal nie jest zadowolona z polarnej stereograficznej projekcji, której używam (EPSG: 3413). Zobacz wiersze 1570-1582 pliku netcdfdataset.cpp:
Moja projekcja ma określoną szerokość geograficzną, ale nie ma żadnych standardowych podobieństw, zgodnie z oczekiwaniami sterownika netcdf.
Odpowiedzi:
Oto kod Pythona, który robi to, co chcesz, odczytując pliki GDAL, które reprezentują dane w określonym czasie i zapisując w jednym pliku NetCDF, który jest zgodny z CF
Python GDAL i NetCDF4 może być trochę trudny do zbudowania, ale dobrą wiadomością jest to, że są częścią większości naukowych dystrybucji python (Python (x, y), Enthought Python Distribution, Anaconda, ...)
Aktualizacja: Nie robiłem jeszcze stereografii biegunowej w NetCDF zgodnym z CF, ale powinienem wyglądać mniej więcej tak. Tutaj założyłem, że
central_meridian
iwlatitude_of_origin
GDAL są takie same jakstraight_vertical_longitude_from_pole
iwlatitude_of_projection_origin
CF:źródło
Łatwo jest umieścić je w jednym NetCDF z narzędziami GDAL, przykład poniżej. Ale nie otrzymujesz wymiaru czasowego / innych metadanych odpowiedzi @ RichSignell. Tiffy zostają po prostu wrzucone do subdanych.
źródło