перепроектировать netcdf в синусоидальную проекцию - PullRequest
Новая
       6

перепроектировать netcdf в синусоидальную проекцию

0 голосов
/ 07 февраля 2019

Я хочу, чтобы глобальный файл был разбит на несколько файлов netcdf, соответствующих плиткам modis.

В настоящее время у меня есть каталог глобальных файлов netcdf со временем.

ПРИМЕЧАНИЕ. Временные данные ERA включенырегулярная сетка широта / долгота.

Долготы варьируются от 0 до 360, что эквивалентно от -180 до +180 в географических системах координат. LINK

Вот так выглядит мой каталог: 

ERAIN_SFC00_6H_10mWind.nc  ERAIN_SFC00_6H_2mtemperature.nc  ERAIN_SFC00_6H_TCW-O3.nc

Мне нужен набор файлов в sinusoidal проекция, один файл для каждой плитки MODIS.Вот так: 

ERAIN_SFC00_6H_TCW-O3_h07v06.nc  ERAIN_SFC00_6H_TCW-O3_h24v03.nc         ERAIN_SFC00_6H_10mWind_H15V05.nc  ERAIN_SFC00_6H_10mWind_H31V06.nc
ERAIN_SFC00_6H_TCW-O3_h08v03.nc  ERAIN_SFC00_6H_TCW-O3_h24v04.nc         ERAIN_SFC00_6H_10mWind_H15V07.nc  ERAIN_SFC00_6H_10mWind_H31V07.nc
ERAIN_SFC00_6H_TCW-O3_h08v04.nc  ERAIN_SFC00_6H_TCW-O3_h24v05.nc         ERAIN_SFC00_6H_10mWind_H16V05.nc  ERAIN_SFC00_6H_10mWind_H31V08.nc
ERAIN_SFC00_6H_TCW-O3_h08v05.nc  ERAIN_SFC00_6H_TCW-O3_h24v06.nc         ERAIN_SFC00_6H_10mWind_H16V06.nc  ERAIN_SFC00_6H_10mWind_H31V09.nc
ERAIN_SFC00_6H_TCW-O3_h08v06.nc  ERAIN_SFC00_6H_TCW-O3_h24v07.nc         ERAIN_SFC00_6H_10mWind_H16V07.nc  ERAIN_SFC00_6H_10mWind_H31V10.nc
ERAIN_SFC00_6H_TCW-O3_h08v07.nc  ERAIN_SFC00_6H_TCW-O3_h25v02.nc         ERAIN_SFC00_6H_10mWind_H16V08.nc  ERAIN_SFC00_6H_10mWind_H31V11.nc
ERAIN_SFC00_6H_TCW-O3_h09v02.nc  ERAIN_SFC00_6H_TCW-O3_h25v03.nc         ERAIN_SFC00_6H_10mWind_H17V02.nc  ERAIN_SFC00_6H_10mWind_H31V12.nc

Псевдокод 

1. Convert the files to sinusoidal projection
2. determine the extent of the MODIS tiles
3. chop the MODIS tile extents and save as new netCDF

Я абсолютно не представляю, а) как преобразовать в синусоидальную проекцию б) как нарезать плитки модиса.

Я сделал эти экстенты плитки из уже существующих файлов modis netcdf.Я не знаю, точны ли они (есть ли официальный источник, определяющий плитки MODIS?), Но у них есть значения.

Данные

ПРИМЕРНЫЕ ДАННЫЕ ЗДЕСЬ 

import xarray as xr

ds = xr.open_dataset('ERA_test.nc')
ds

Out[]:
<xarray.Dataset>
Dimensions:        (lat: 256, lon: 512, time: 5)
Coordinates:
  * lon            (lon) float64 0.0 0.7031 1.406 ... -2.109 -1.406 -0.7031
  * lat            (lat) float64 89.46 88.77 88.07 ... -88.07 -88.77 -89.46
  * time           (time) datetime64[ns] 2008-01-02 ... 2008-01-03
Data variables:
    temperature2m  (time, lat, lon) float32 ...
Attributes:
    content:  ERA INTERIM data of var167; converted from grb to netCDF with C...
    history:

Размер плитки 

import pickle

tiles = pickle.load(open('tile_extents.pkl','rb'))
tiles['h16v06']

Out[]:
{'lon_min': -23.09401076758503,
 'lon_max': -10.641777724759121,
 'lat_min': 20.0,
 'lat_max': 30.0}

Преобразование в синусоидальную проекцию

Это то, чего я действительно не знаю, как сделать.Я думаю, что нашел проекцию MODIS от pyproj. 

import pyproj as Proj

p_modis_grid = Proj('+proj=sinu +R=6371007.181 +nadgrids=@null +wktext')

Но как мне использовать это для перепроектирования данных netCDF?Рад узнать, что я использую не тот инструмент.rasterio подойдет лучше?Или NCO / CDO?

Спасибо за вашу помощь.

Текущее наилучшее предположение.

Так что это, похоже, меняет проекцию, но также портит все измерения.Я думаю, что это отправная точка, но откуда?

gdalwarp -t_srs "+proj=sinu +R=6371007.181 +nadgrids=@null +wktext" ERA_test.nc out.nc

$ ncdump -h out.nc

netcdf out {
dimensions:
    x = 512 ;
    y = 256 ;
variables:
    short Band1(y, x) ;
        Band1:long_name = "GDAL Band Number 1" ;
        Band1:_FillValue = -32767s ;
        Band1:grid_mapping = "sinusoidal" ;
        Band1:units = "K" ;
    short Band2(y, x) ;
        Band2:long_name = "GDAL Band Number 2" ;
        Band2:_FillValue = -32767s ;
        Band2:grid_mapping = "sinusoidal" ;
        Band2:units = "K" ;
    short Band3(y, x) ;
        Band3:long_name = "GDAL Band Number 3" ;
        Band3:_FillValue = -32767s ;
        Band3:grid_mapping = "sinusoidal" ;
        Band3:units = "K" ;
    short Band4(y, x) ;
        Band4:long_name = "GDAL Band Number 4" ;
        Band4:_FillValue = -32767s ;
        Band4:grid_mapping = "sinusoidal" ;
        Band4:units = "K" ;
    short Band5(y, x) ;
        Band5:long_name = "GDAL Band Number 5" ;
        Band5:_FillValue = -32767s ;
        Band5:grid_mapping = "sinusoidal" ;
        Band5:units = "K" ;
    char sinusoidal ;
        sinusoidal:grid_mapping_name = "sinusoidal" ;
        sinusoidal:false_easting = 0. ;
        sinusoidal:false_northing = 0. ;
        sinusoidal:longitude_of_central_meridian = 0. ;
        sinusoidal:long_name = "CRS definition" ;
        sinusoidal:longitude_of_prime_meridian = 0. ;
        sinusoidal:semi_major_axis = 6371007.181 ;
        sinusoidal:inverse_flattening = 0. ;
        sinusoidal:spatial_ref = "PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6371007.181,0],EXTENSION[\"PROJ4_GRIDS\",\"@null\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Sinusoidal\"],PARAMETER[\"longitude_of_center\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],EXTENSION[\"PROJ4\",\"+proj=sinu +R=6371007.181 +nadgrids=@null +wktext\"]]" ;
        sinusoidal:GeoTransform = "-360.3515625 0.7028340789935529 0 89.81365616296401 0 -0.7028340789935529 " ;
    double x(x) ;
        x:standard_name = "projection_x_coordinate" ;
        x:long_name = "x coordinate of projection" ;
        x:units = "m" ;
    double y(y) ;
        y:standard_name = "projection_y_coordinate" ;
        y:long_name = "y coordinate of projection" ;
        y:units = "m" ;

// global attributes:
        :Conventions = "CF-1.5" ;
        :GDAL = "GDAL 2.3.3, released 2018/12/14" ;
        :history = "Thu Feb 07 14:27:39 2019: GDAL Create( out.nc, ... )" ;
}
...