Преобразование стека временных рядов растра GTiff в один NetCDF

12

Переезд из списка рассылки gdal-dev:

В понедельник, 2 сентября 2013 года, в 19:09 Дэвид Шин написал:

Привет список, я пытаюсь упаковать временные ряды растров GTiff с идентичной проекцией / экстентом / разрешением как один файл NetCDF для распространения. Я провел последний час, консультируясь с онлайн-документом и безуспешно играя с gdal_translate, gdalbuildvrt и gdalwarp.

Есть ли простой способ сделать это, используя существующие утилиты командной строки gdal? Я решил спросить, прежде чем прибегнуть к индивидуальному решению с использованием NetCDF Python API.

Благодарю. -Давид

В четверг, 3 сентября 2013 г., в 10:15 Этьен Туриньи написал:

то, что вы хотите, вероятно, выходит за рамки gdal. Это потребовало бы умного управления метаданными, чтобы gdal_translate поместил их в один файл ...

Я бы посоветовал вам преобразовать их все в netcdf, используя gdal_translate, а затем использовать python-netcdf4 (не тот, что из numpy / scipy), чтобы сложить их во временное измерение.

Во вторник, 3 сентября 2013 г., в 7:55 утра «Сигнелл, Ричард» написал:

Дэвид, если вы разместите свой вопрос в группе обмена ГИС-стеками /gis//, я приведу пример кода, который должен быть полезным.

-Богатый

====================

Обновление 03.09.13 17:04 PDT

Вот вывод gdalinfo для одного из моих входных наборов данных:


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

Вслед за предложенным Люком подходом.

Генерация vrt работает нормально:

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>
...

Но когда я пытаюсь перевести на nc, я получаю следующую ошибку:


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)

Таким образом, при ближайшем рассмотрении кажется, что gdal недоволен полярной стереографической проекцией, которую я использую (EPSG: 3413). Смотрите строки 1570-1582 из netcdfdataset.cpp:

https://code.vpac.org/gitorious/gdal-netcdf-testing/gdal-netcdf-driver/blobs/8fa3582669969ad4d55e461f5846b3ed33727f63/gdal/frmts/netcdf/netcdfdataset.cpp

В моей проекции указан latitude_of_origin, но нет стандартных параллелей, как ожидается драйвером netcdf.

python gdal netcdf Дэвид Шин
источник
1
Какая версия GDAL? В драйвере NetCDF произошел ряд изменений в GDAL> = 1.9.0. На этой странице конкретно упоминаются изменения в обработке полярной стереографической проекции. Вы можете обойти эту проблему, переопределив проекцию с помощью параметра gdal_translate -a_srs и указав допустимую, но эквивалентную строку проекции. См. Также ( trac.osgeo.org/gdal/wiki/NetCDF_ProjectionTestingStatus )
user2856
gdalinfo - версия GDAL 1.11dev, выпущена 2013/04/13
Дэвид Шин
1
Спасибо и Ричу, и Люку за полезный вклад. Мне нужно обновиться до последней версии GDAL, оценить новейшую полярную стереографическую функциональность драйвера netcdf и проконсультироваться с gdal-dev по любым нерешенным вопросам. Хотя оба ответа будут работать, мне нравится рецепт Рича, и я приму его в своих целях. Я знаю, что другие найдут эту дискуссию полезной - рад, что она заархивирована на SE.
Дэвид Шин

Ответы:

22

Вот некоторый код Python, который делает то, что вы хотите, читая файлы GDAL, которые представляют данные в определенное время, и записывая в один файл NetCDF, который является CF-совместимым

#!/usr/bin/env python
'''
Convert a bunch of GDAL readable grids to a NetCDF Time Series.
Here we read a bunch of files that have names like:
/usgs/data0/prism/1890-1899/us_tmin_1895.01
/usgs/data0/prism/1890-1899/us_tmin_1895.02
...
/usgs/data0/prism/1890-1899/us_tmin_1895.12
'''

import numpy as np
import datetime as dt
import os
import gdal
import netCDF4
import re

ds = gdal.Open('/usgs/data0/prism/1890-1899/us_tmin_1895.01')
a = ds.ReadAsArray()
nlat,nlon = np.shape(a)

b = ds.GetGeoTransform() #bbox, interval
lon = np.arange(nlon)*b[1]+b[0]
lat = np.arange(nlat)*b[5]+b[3]


basedate = dt.datetime(1858,11,17,0,0,0)

# create NetCDF file
nco = netCDF4.Dataset('time_series.nc','w',clobber=True)

# chunking is optional, but can improve access a lot: 
# (see: http://www.unidata.ucar.edu/blogs/developer/entry/chunking_data_choosing_shapes)
chunk_lon=16
chunk_lat=16
chunk_time=12

# create dimensions, variables and attributes:
nco.createDimension('lon',nlon)
nco.createDimension('lat',nlat)
nco.createDimension('time',None)
timeo = nco.createVariable('time','f4',('time'))
timeo.units = 'days since 1858-11-17 00:00:00'
timeo.standard_name = 'time'

lono = nco.createVariable('lon','f4',('lon'))
lono.units = 'degrees_east'
lono.standard_name = 'longitude'

lato = nco.createVariable('lat','f4',('lat'))
lato.units = 'degrees_north'
lato.standard_name = 'latitude'

# create container variable for CRS: lon/lat WGS84 datum
crso = nco.createVariable('crs','i4')
csro.long_name = 'Lon/Lat Coords in WGS84'
crso.grid_mapping_name='latitude_longitude'
crso.longitude_of_prime_meridian = 0.0
crso.semi_major_axis = 6378137.0
crso.inverse_flattening = 298.257223563

# create short integer variable for temperature data, with chunking
tmno = nco.createVariable('tmn', 'i2',  ('time', 'lat', 'lon'), 
   zlib=True,chunksizes=[chunk_time,chunk_lat,chunk_lon],fill_value=-9999)
tmno.units = 'degC'
tmno.scale_factor = 0.01
tmno.add_offset = 0.00
tmno.long_name = 'minimum monthly temperature'
tmno.standard_name = 'air_temperature'
tmno.grid_mapping = 'crs'
tmno.set_auto_maskandscale(False)

nco.Conventions='CF-1.6'

#write lon,lat
lono[:]=lon
lato[:]=lat

pat = re.compile('us_tmin_[0-9]{4}\.[0-9]{2}')
itime=0

#step through data, writing time and data to NetCDF
for root, dirs, files in os.walk('/usgs/data0/prism/1890-1899/'):
    dirs.sort()
    files.sort()
    for f in files:
        if re.match(pat,f):
            # read the time values by parsing the filename
            year=int(f[8:12])
            mon=int(f[13:15])
            date=dt.datetime(year,mon,1,0,0,0)
            print(date)
            dtime=(date-basedate).total_seconds()/86400.
            timeo[itime]=dtime
           # min temp
            tmn_path = os.path.join(root,f)
            print(tmn_path)
            tmn=gdal.Open(tmn_path)
            a=tmn.ReadAsArray()  #data
            tmno[itime,:,:]=a
            itime=itime+1

nco.close()

GDAL и NetCDF4 Python могут быть немного трудны для сборки, но хорошая новость заключается в том, что они являются частью большинства научных дистрибутивов Python (Python (x, y), Enthought Python Distribution, Anaconda, ...)

Обновление: я еще не делал полярную стереографию в CFCD-совместимом NetCDF, но я должен выглядеть примерно так. Здесь я предположил, что central_meridianи latitude_of_originв GDAL такие же, как straight_vertical_longitude_from_poleи latitude_of_projection_originв CF:

#!/usr/bin/env python
'''
Convert a bunch of GDAL readable grids to a NetCDF Time Series.
Here we read a bunch of files that have names like:
/usgs/data0/prism/1890-1899/us_tmin_1895.01
/usgs/data0/prism/1890-1899/us_tmin_1895.02
...
/usgs/data0/prism/1890-1899/us_tmin_1895.12
'''

import numpy as np
import datetime as dt
import os
import gdal
import netCDF4
import re

ds = gdal.Open('/usgs/data0/prism/1890-1899/us_tmin_1895.01')
a = ds.ReadAsArray()
ny,nx = np.shape(a)

b = ds.GetGeoTransform() #bbox, interval
x = np.arange(nx)*b[1]+b[0]
y = np.arange(ny)*b[5]+b[3]


basedate = dt.datetime(1858,11,17,0,0,0)

# create NetCDF file
nco = netCDF4.Dataset('time_series.nc','w',clobber=True)

# chunking is optional, but can improve access a lot: 
# (see: http://www.unidata.ucar.edu/blogs/developer/entry/chunking_data_choosing_shapes)
chunk_x=16
chunk_y=16
chunk_time=12

# create dimensions, variables and attributes:
nco.createDimension('x',nx)
nco.createDimension('y',ny)
nco.createDimension('time',None)
timeo = nco.createVariable('time','f4',('time'))
timeo.units = 'days since 1858-11-17 00:00:00'
timeo.standard_name = 'time'

xo = nco.createVariable('x','f4',('x'))
xo.units = 'm'
xo.standard_name = 'projection_x_coordinate'

yo = nco.createVariable('y','f4',('y'))
yo.units = 'm'
yo.standard_name = 'projection_y_coordinate'

# create container variable for CRS: x/y WGS84 datum
crso = nco.createVariable('crs','i4')
crso.grid_mapping_name='polar_stereographic'
crso.straight_vertical_longitude_from_pole = -45.
crso.latitude_of_projection_origin = 70.
crso.scale_factor_at_projection_origin = 1.0
crso.false_easting = 0.0
crso.false_northing = 0.0
crso.semi_major_axis = 6378137.0
crso.inverse_flattening = 298.257223563

# create short integer variable for temperature data, with chunking
tmno = nco.createVariable('tmn', 'i2',  ('time', 'y', 'x'), 
   zlib=True,chunksizes=[chunk_time,chunk_y,chunk_x],fill_value=-9999)
tmno.units = 'degC'
tmno.scale_factor = 0.01
tmno.add_offset = 0.00
tmno.long_name = 'minimum monthly temperature'
tmno.standard_name = 'air_temperature'
tmno.grid_mapping = 'crs'
tmno.set_auto_maskandscale(False)

nco.Conventions='CF-1.6'

#write x,y
xo[:]=x
yo[:]=y

pat = re.compile('us_tmin_[0-9]{4}\.[0-9]{2}')
itime=0

#step through data, writing time and data to NetCDF
for root, dirs, files in os.walk('/usgs/data0/prism/1890-1899/'):
    dirs.sort()
    files.sort()
    for f in files:
        if re.match(pat,f):
            # read the time values by parsing the filename
            year=int(f[8:12])
            mon=int(f[13:15])
            date=dt.datetime(year,mon,1,0,0,0)
            print(date)
            dtime=(date-basedate).total_seconds()/86400.
            timeo[itime]=dtime
           # min temp
            tmn_path = os.path.join(root,f)
            print(tmn_path)
            tmn=gdal.Open(tmn_path)
            a=tmn.ReadAsArray()  #data
            tmno[itime,:,:]=a
            itime=itime+1

nco.close()
Рич Синьелл
источник
Отличный код Рич! Это очень полезно, и я буду использовать это в будущем. Похоже, ваша входная проекция считается географической с единицами широты / долготы (EPSG: 4326). Я работаю с данными высокого разрешения в полярных широтах, так что это не идеально, но я постараюсь перейти на WGS84.
Дэвид Шин
широта / долгота была только примером. Вы можете использовать все, что пожелаете. На какие приложения вы ориентируетесь? ArcGIS, просто для архивации или как?
Рич Сигнелл
Ну, у меня есть много подобных временных рядов, и я оцениваю варианты для эффективного хранения и анализа. Но в данный момент я упаковываю данные для приема по моделям потока. Модельное сообщество, по крайней мере, моделирование ледяного потока, похоже, любит netcdf.
Дэвид Шин
Есть ли URL, где мы могли бы найти образец этих данных?
Рич Сигнелл
К сожалению, я не могу распространять в данный момент, но есть планы архивировать в будущем.
Дэвид Шин
2

Их легко поместить в один NetCDF с утилитами GDAL, пример ниже. Но вы не получите временное измерение / другие метаданные ответа @ RichSignell. TIFF просто сбрасывается в наборы данных.

C:\remotesensing\testdata>dir /b ndvi*.tif
ndvi1.tif
ndvi2.tif
ndvi3.tif

C:\remotesensing\testdata>gdalbuildvrt -separate ndvi.vrt ndvi*.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\remotesensing\testdata>gdal_translate -of netcdf ndvi.vrt ndvi.nc
Input file size is 96, 88
0...10...20...30...40...50...60...70...80...90...100 - done.

C:\remotesensing\testdata>gdalinfo ndvi.nc
Driver: netCDF/Network Common Data Format
Files: ndvi.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 1.10.0, released 2013/04/24
  NC_GLOBAL#history=Wed Sep 04 09:49:11 2013: GDAL CreateCopy( ndvi.nc, ... )
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"ndvi.nc":Band1
  SUBDATASET_1_DESC=[88x96] Band1 (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"ndvi.nc":Band2
  SUBDATASET_2_DESC=[88x96] Band2 (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"ndvi.nc":Band3
  SUBDATASET_3_DESC=[88x96] Band3 (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

C:\remotesensing\testdata>gdalinfo NETCDF:"ndvi.nc":Band1
Driver: netCDF/Network Common Data Format
Files: ndvi.nc
Size is 96, 88
Coordinate System is:
GEOGCS["GCS_GDA_1994",
    DATUM["Geocentric_Datum_of_Australia_1994",
        SPHEROID["GRS 1980",6378137,298.2572221010002,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6283"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (115.810500000000000,-32.260249999999999)
Pixel Size = (0.000250000000000,-0.000250000000000)
Metadata:
  Band1#_FillValue=0
  Band1#grid_mapping=crs
  Band1#long_name=GDAL Band Number 1
  crs#GeoTransform=115.8105 0.00025 0 -32.26025 0 -0.00025
  crs#grid_mapping_name=latitude_longitude
  crs#inverse_flattening=298.2572221010002
  crs#longitude_of_prime_meridian=0
  crs#semi_major_axis=6378137
  crs#spatial_ref=GEOGCS["GCS_GDA_1994",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 1.10.0, released 2013/04/24
  NC_GLOBAL#history=Wed Sep 04 09:49:11 2013: GDAL CreateCopy( ndvi.nc, ... )
Corner Coordinates:
Upper Left  ( 115.8105000, -32.2602500) (115d48'37.80"E, 32d15'36.90"S)
Lower Left  ( 115.8105000, -32.2822500) (115d48'37.80"E, 32d16'56.10"S)
Upper Right ( 115.8345000, -32.2602500) (115d50' 4.20"E, 32d15'36.90"S)
Lower Right ( 115.8345000, -32.2822500) (115d50' 4.20"E, 32d16'56.10"S)
Center      ( 115.8225000, -32.2712500) (115d49'21.00"E, 32d16'16.50"S)
Band 1 Block=96x1 Type=Float32, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    _FillValue=0
    grid_mapping=crs
    long_name=GDAL Band Number 1
    NETCDF_VARNAME=Band1
user2856
источник
Я попробовал этот подход, и он не удался для моих входных данных - я опубликую вывод выше.
Дэвид Шин
В качестве теста я использовал gdalwarp для перепроектирования многоканального vrt EPSG: 3413 в EPSG: 4326, а затем использовал gdal_translate для преобразования в netcdf4. Как предполагает Люк, это работает без проблем. Как предложил Этьен в оригинальной ветке gdal-dev, для этого подхода ограничен контроль над метаданными.
Дэвид Шин