Converta a pilha de séries temporais do raster GTiff em um único NetCDF

12

Movendo-se da lista de discussão gdal-dev:

Em segunda-feira, 2 de setembro de 2013 às 19:09, David Shean escreveu:

Olá, estou tentando empacotar séries temporais de rasters GTiff com projeção / extensão / resolução idênticas como um único arquivo NetCDF para distribuição. Passei a última hora consultando o documento on-line e jogando com gdal_translate, gdalbuildvrt e gdalwarp sem sucesso.

Existe uma maneira fácil de fazer isso usando os utilitários de linha de comando gdal existentes? Pensei em perguntar antes de recorrer a uma solução personalizada usando a API NetCDF Python.

Obrigado. -David

Em terça-feira, 3 de setembro de 2013 às 10:15, Etienne Tourigny escreveu:

o que você quer provavelmente está fora do escopo da gdal. Seria necessário um gerenciamento inteligente de metadados para que o gdal_translate os coloque em um único arquivo ...

Aconselho que você os converta todos em netcdf usando gdal_translate e, em seguida, use python-netcdf4 (não o de numpy / scipy) para empilhá-los na dimensão temporal.

Em terça-feira, 3 de setembro de 2013, às 07:55, "Signell, Richard" escreveu:

David, Se você postar sua pergunta no grupo GIS stackexchange /gis// , fornecerei um código de exemplo que deve ser útil.

-Rico

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

Atualização 9/3/13 17:04 PDT

Aqui está a saída gdalinfo para um dos meus conjuntos de dados de entrada:


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

Seguindo a abordagem sugerida por Luke.

A geração vrt funciona bem:

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

Mas quando tento traduzir para nc, recebo o seguinte erro:


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)

Assim, após uma inspeção mais detalhada, parece que gdal está descontente com a projeção estereográfica polar que estou usando (EPSG: 3413). Veja as linhas 1570-1582 de netcdfdataset.cpp:

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

Minha projeção tem latitude_of_origin especificada, mas nenhum paralelo padrão conforme o esperado pelo driver netcdf.

python gdal netcdf David Shean
fonte
1
Qual versão do GDAL? Houve várias alterações no driver NetCDF no GDAL> = 1.9.0. Essa página menciona especificamente alterações no manuseio da projeção estereográfica polar. Você pode solucionar o problema substituindo a projeção pelo parâmetro gdal_translate -a_srs e especificando uma sequência de projeção válida, mas equivalente. Veja também ( trac.osgeo.org/gdal/wiki/NetCDF_ProjectionTestingStatus )
user2856
gdalinfo --version GDAL 1.11dev, lançado em 13/04/2013
David Shean
1
Agradecemos a Rich e Luke pela contribuição útil. Preciso atualizar para a versão mais recente do GDAL, avaliar a funcionalidade estereográfica polar do driver netcdf e acompanhar o gdal-dev sobre quaisquer problemas remanescentes. Embora ambas as respostas funcionem, gosto da receita de Rich e adotarei para meus próprios propósitos. Sei que outras pessoas acharão útil essa discussão - feliz por ter sido arquivada no SE.
David Shean

Respostas:

21

Aqui está um código python que faz o que você deseja, lendo arquivos GDAL que representam dados em horários específicos e gravando em um único arquivo NetCDF compatível com 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 e NetCDF4 Python podem ser um pouco trabalhosos, mas a boa notícia é que eles fazem parte da maioria das distribuições científicas de python (Python (x, y), Enthought Python Distribution, Anaconda, ...)

Atualização: ainda não fiz o estereográfico polar no NetCDF compatível com CF, mas devo parecer algo assim. Aqui eu assumi que central_meridiane latitude_of_originno GDAL são os mesmos que straight_vertical_longitude_from_polee latitude_of_projection_originno 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()
Rich Signell
fonte
Ótimo código Rich! Isso é muito útil, e usarei isso no futuro. Parece que sua projeção de entrada é considerada geográfica com unidades de lat / lon (EPSG: 4326). Estou trabalhando com dados de alta resolução em latitudes polares, então isso não é o ideal, mas tentarei converter para WGS84.
David Shean
lat / lon foi apenas um exemplo. Você pode usar o que quiser. Quais aplicativos você está direcionando? ArcGIS, apenas para arquivamento ou o quê?
Rico Signell
Bem, tenho muitas séries temporais como esta e estou avaliando opções para armazenamento e análise eficientes. No momento, estou empacotando dados para ingestão por modelos de fluxo. A comunidade de modelagem, pelo menos modelagem de fluxo de gelo, parece gostar de netcdf.
David Shean
Existe um URL onde possamos encontrar uma amostra desses dados?
Rico Signell
Infelizmente, não posso distribuir neste momento, mas há planos para arquivar no futuro.
David Shean
2

É fácil colocá-los em um único NetCDF com utilitários GDAL, exemplo abaixo. Mas você não recebe a dimensão temporal / outros metadados da resposta do @ RichSignell. Os tiffs são despejados em subdatasets.

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
fonte
Eu tentei essa abordagem e ela falhou nos meus dados de entrada - postarei a saída acima.
David Shean
Como teste, usei o gdalwarp para reprojetar o vrt multi-banda EPSG: 3413 para o EPSG: 4326 e, em seguida, usei o gdal_translate para converter para netcdf4. Como Lucas sugere, isso funciona sem problemas. Como Etienne sugeriu no segmento gdal-dev original, há um controle limitado sobre os metadados para essa abordagem.
David Shean