Правильный подход к регриддингу в GDAL

Ответить
Аватара пользователя
Jasen
Активный участник
Сообщения: 100
Зарегистрирован: 27 янв 2006, 18:33
Репутация: 4
Ваше звание: Специалист
Откуда: Москва
Контактная информация:

Правильный подход к регриддингу в GDAL

Сообщение Jasen » 18 май 2020, 16:21

Привет. Комплексный вопрос, потому что уже все перерыл. Задача проста: из src сделать dst в сетке и охвате mask.
В GDAL пришел к этому:

Код: Выделить всё

mask=gdal.Open(path1,gdal.GA_ReadOnly)
src=gdal.Open(path2,gdal.GA_ReadOnly)
dst=gdal.GetDriverByName('NetCDF').Create(pathdst,xsize=mask.RasterXSize,ysize=mask.RasterYSize)
dst.SetGeoTransform(mask.GetGeoTransform())
dst.SetProjection(mask.GetProjection())
optt=gdal.WarpOptions(format='NetCDF',xRes=mask.RasterXSize,yRes=mask.RasterYSize)
dst=gdal.Warp(dst,src,options=optt)
dst=None
src=None
mask=None
но оно генерирует файл пустого типа, короче - неправильно.
Добился результата с помощью Reprojectimage, но там обрабатывается только 1 слой, а я хочу как-то для всех сразу сделать warp. С Translate та же история (это логично).
Итак, вопрос: как с помощью Python и GDAL сделать перепроецирование растра в сетку другого растра для всех каналов (слоев времени в данном случае) сразу?
Если кто-то найдет в интернете ссылку по этому вопросу, которую я ещё не видел - тому мой огромный респект. Но я пересмотрел уже все что только можно, в том числе документацию. Оно не работает.
Что написал - то написал!

Аватара пользователя
Эдуард Казаков
Гуру
Сообщения: 546
Зарегистрирован: 23 апр 2014, 17:11
Репутация: 532
Откуда: Planet Earth
Контактная информация:

Re: Правильный подход к регриддингу в GDAL

Сообщение Эдуард Казаков » 18 май 2020, 18:45

Здравствуйте, может поможет
https://github.com/eduard-kazakov/RasterAdjuster
Задача немного другая - быстро привести растр1 в пространственный домен растра 2, но разница не большая, можно код посмотреть.

Из общих вещей:
- при запуске warp желательно указывать параметры, целиком описывающие пространственный домен, то есть разрешение+экстент или размер в пикселях + экстент. Чтобы не возникло неопределенностей в общем.
- Зачем создавать датасет предварительно? Вы можете запускать так
dst=gdal.Warp(<путь до целевого файла>,src,options=opt)
А в opt собственно описать всё: целевую проекцию, xRes, yRes, экстент, формат.

Аватара пользователя
Jasen
Активный участник
Сообщения: 100
Зарегистрирован: 27 янв 2006, 18:33
Репутация: 4
Ваше звание: Специалист
Откуда: Москва
Контактная информация:

Re: Правильный подход к регриддингу в GDAL

Сообщение Jasen » 18 май 2020, 20:38

Так в том и проблема, что у меня нет проекции в файлах, есть только GetGeoTransform, а GetProjection пустой.

Код: Выделить всё

mask=gdal.Open(path1,gdal.GA_ReadOnly)#.GetRasterBand(1)
print(type(mask))
src=gdal.Open(path2,gdal.GA_ReadOnly).GetRasterBand(1)
bnds=(mask.GetGeoTransform()[0],mask.GetGeoTransform()[3]-mask.RasterYSize,mask.GetGeoTransform()[0]+mask.RasterXSize,mask.GetGeoTransform()[3])
optt=gdal.WarpOptions(format='NetCDF',xRes=mask.RasterXSize,yRes=mask.RasterYSize,outputBounds=bnds,width=mask.GetGeoTransform()[1])
dst=gdal.Warp(pathdst,src,options=optt)
dst=None
src=None
mask=None
Если для привязки использую Dataset - ошибки нет, но файла нет тоже (dst - NoneType). Если вычленяю из файла mask первый слой, то есть использую RasterBand, по нему не хочет почему-то определяться GetGeoTransform().
В моем src файле 16000 слоев (моментов времени), да, но я полагал, что GDAL достаточно низкоуровневая штука, чтобы справиться, не?
Последний раз редактировалось Jasen 18 май 2020, 20:48, всего редактировалось 1 раз.
Что написал - то написал!

Аватара пользователя
Jasen
Активный участник
Сообщения: 100
Зарегистрирован: 27 янв 2006, 18:33
Репутация: 4
Ваше звание: Специалист
Откуда: Москва
Контактная информация:

Re: Правильный подход к регриддингу в GDAL

Сообщение Jasen » 18 май 2020, 20:40

Или мне сделать через ReprojectImage в цикле 16000 слоев? Почему-то мне кажется, что это не быстрое будет дело :)
Что написал - то написал!

Аватара пользователя
antonv
Активный участник
Сообщения: 229
Зарегистрирован: 29 ноя 2016, 10:44
Репутация: 114
Откуда: Санкт-Петербург

Re: Правильный подход к регриддингу в GDAL

Сообщение antonv » 19 май 2020, 11:52

Jasen,
Если вычленяю из файла mask первый слой, то есть использую RasterBand, по нему не хочет почему-то определяться GetGeoTransform().
GetGeoTransform - это ведь метод самого датасета, а не вложенных bands, так и должно быть.
А у вас на вход точно не пустой src подаётся? После того, как вы сделали src=gdal.Open(path2,gdal.GA_ReadOnly).GetRasterBand(1), можете удостовериться при помощи src.ReadAsArray(), что внутри не пусто?
Какого формата файл на входе? Если в src=gdal.Open(...) вы открываете NetCDF файл, то всё немного хитрее, у этого формата есть subdatasets. Открывать надо конкретную переменную командой gdal.Open('NETCDF:"'+datafile+'":variable').

Аватара пользователя
Jasen
Активный участник
Сообщения: 100
Зарегистрирован: 27 янв 2006, 18:33
Репутация: 4
Ваше звание: Специалист
Откуда: Москва
Контактная информация:

Re: Правильный подход к регриддингу в GDAL

Сообщение Jasen » 19 май 2020, 12:17

Не знаю. Сегодня попробовал в QGis с теми же файлами сделать то же самое. Затем перенес все параметры командной строки в Python. Результат тот же. В QGis получилась прекрасная привязка. Значит, дело не в файлах, а в Python.
Командная строка из окошка QGis:

Код: Выделить всё

gdalwarp -ot Float32 -t_srs EPSG:4326 -r bilinear -of GTiff -tr 1 1 -te 19.5 " 40.0" " 70.5" " 75.0" -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -wo OPTIMIZE_SIZE=TRUE C:\Users\jasinskiy\Documents\WORK\mask.nc
Новый код:

Код: Выделить всё

mask=gdal.Open(path1,gdal.GA_ReadOnly)
src=gdal.Open(os.getcwd()+'/test1/mask.nc',gdal.GA_ReadOnly)
bnds=(mask.GetGeoTransform()[0],mask.GetGeoTransform()[3]-mask.RasterYSize,mask.GetGeoTransform()[0]+mask.RasterXSize,mask.GetGeoTransform()[3])
optt=gdal.WarpOptions(format='netCDF',resampleAlg='bilinear',warpOptions=('OPTIMIZE_SIZE=TRUE'),creationOptions=('COMPRESS=DEFLATE','PREDICTOR=1','ZLEVEL=6'),outputType=gdalconst.GDT_Float32,srcSRS='EPSG:4326',dstSRS='EPSG:4326',xRes=mask.RasterXSize,yRes=mask.RasterYSize,outputBounds=bnds,width=mask.GetGeoTransform()[1])
print(mask.RasterXSize,mask.RasterYSize,mask.GetGeoTransform()[1])
dst=gdal.Warp(pathdst,src,options=optt)
Что написал - то написал!

Аватара пользователя
Jasen
Активный участник
Сообщения: 100
Зарегистрирован: 27 янв 2006, 18:33
Репутация: 4
Ваше звание: Специалист
Откуда: Москва
Контактная информация:

Re: Правильный подход к регриддингу в GDAL

Сообщение Jasen » 19 май 2020, 12:21

antonv писал(а):
19 май 2020, 11:52
Какого формата файл на входе? Если в src=gdal.Open(...) вы открываете NetCDF файл, то всё немного хитрее, у этого формата есть subdatasets. Открывать надо конкретную переменную командой gdal.Open('NETCDF:"'+datafile+'":variable').
Спасибо за замечание, про band я действительно не сразу понял что метод не годится. Но по сути GetGeoTransform из целого netCDF файла у меня работает, он ведь печатает массив. Значит данные нет проблем получить. А на выходе мне нужно весь большой nc файл перепроецировать, что я благополучно сделал в QGis. То есть subdatasets тут получать наверное не нужно?
Что написал - то написал!

Ответить

Вернуться в «GDAL/OGR»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 3 гостя