Деформированный VRT, не трансформирующийся - PullRequest
Новая
       5

Деформированный VRT, не трансформирующийся

0 голосов
/ 31 октября 2018

Я работаю с данными nlcd и cdl в системе HPC. Я хочу деформировать данные nlcd в растр cdl. Я надеялся использовать файлы vrt, поскольку хранилище ограничено, но у меня есть доступ к большому количеству вычислительной мощности.

У меня есть данные cdl, хранящиеся в файле zarr с информацией о преобразовании и crs: 

dat_cdl = xr.open_zarr('/project/cper_neon_aop/regionalization/data/cdl/cdl_data.zarr/').sel(year=2011)#.persist()

<xarray.Dataset>
Dimensions:   (x: 153811, y: 96523)
Coordinates:
  * x         (x) float64 -2.356e+06 -2.356e+06 ... 2.258e+06 2.258e+06
  * y         (y) float64 3.173e+06 3.173e+06 3.173e+06 ... 2.77e+05 2.769e+05
    year      int64 2011
Data variables:
    CDL_Data  (y, x) uint8 dask.array<shape=(96523, 153811), chunksize=(2500, 2500)>

CRS, преобразование и форма этих данных: 

print(rasterio.Affine(*dat_cdl['CDL_Data'].attrs['transform']))
print(dat_cdl['CDL_Data'].attrs['crs'])
print(dat_cdl['CDL_Data'].shape[0])
print(dat_cdl['CDL_Data'].shape[1])

| 30.00, 0.00,-2356095.00|
| 0.00,-30.00, 3172605.00|
| 0.00, 0.00, 1.00|
+ellps=GRS80 +lat_0=23 +lat_1=29.5 +lat_2=45.5 +lon_0=-96 +no_defs +proj=aea +towgs84=0,0,0,-0,-0,-0,0 +units=m +x_0=0 +y_0=0
96523
153811

И я пытаюсь преобразовать данные nlcd в cdl shape / transform / crs через: 

f_dir = '/project/cper_neon_aop/regionalization/data/nlcd/nlcd_2011_landcover_2011_edition_2014_10_10/'
nlcd_f = 'nlcd_2011_landcover_2011_edition_2014_10_10.img'
from rasterio.vrt import WarpedVRT as rio_WarpedVRT
with rasterio.open(f_dir+nlcd_f) as src:
    print(src.profile)
    with rio_WarpedVRT(src,transfrom=rasterio.Affine(*dat_cdl['CDL_Data'].attrs['transform']),crs=dat_cdl['CDL_Data'].attrs['crs'],height=dat_cdl['CDL_Data'].shape[0],width=dat_cdl['CDL_Data'].shape[1],NUM_THREADS=2,warp_mem_limit=1024,) as vrt:    
        rio_shutil.copy(vrt, f_dir+nlcd_f[0:-4]+'_cdl_grid.vrt', driver='VRT')
        print(vrt.profile)
        print(vrt.shape)

{'driver': 'HFA', 'dtype': 'uint8', 'nodata': None, 'width': 161190, 'height': 104424, 'count': 1, 'crs': CRS({'proj': 'aea', 'lat_1': 29.5, 'lat_2': 45.5, 'lat_0': 23, 'lon_0': -96, 'x_0': 0, 'y_0': 0, 'ellps': 'GRS80', 'towgs84': '0,0,0,-0,-0,-0,0', 'units': 'm', 'no_defs': True}), 'transform': Affine(30.0, 0.0, -2493045.0,
       0.0, -30.0, 3310005.0), 'blockxsize': 64, 'blockysize': 64, 'tiled': True}
{'driver': 'VRT', 'dtype': 'uint8', 'nodata': None, 'width': 161190, 'height': 104424, 'count': 1, 'crs': CRS({'ellps': 'GRS80', 'lat_0': 23, 'lat_1': 29.5, 'lat_2': 45.5, 'lon_0': -96, 'no_defs': True, 'proj': 'aea', 'towgs84': '0,0,0,-0,-0,-0,0', 'units': 'm', 'x_0': 0, 'y_0': 0}), 'transform': Affine(30.0, 0.0, -2493045.0,
       0.0, -30.0, 3310005.0), 'blockxsize': 512, 'blockysize': 128, 'tiled': True}
(104424, 161190)

Похоже, это не меняет форму или трансформацию ... Я неправильно интерпретирую функцию WarpedVRT?

Соответствующие версии пакета:
растерио: 1.0.8
гдал: 2.3.1

Спасибо!

...