Как преобразовать файл.grib в GeoTIFF с правильной проекцией, используя Python, GDAL, ArcPy

Question

Я пытаюсь преобразовать файл.grib в GeoTIFF для использования в ГИС (в частности, в ArcGIS), но у меня возникают проблемы с правильным отображением изображения. Я смог создать GeoTIFF, используя GDAL в Python, который показывает данные, но не отображается в правильном месте при переносе в ArcGIS. Полученное изображение ниже.

Данные, с которыми я работаю, можно загрузить по адресу: https://gimms.gsfc.nasa.gov/SMOS/SMAP/L05/

Я пытаюсь проецировать данные в WGS84 Web Mercator (Вспомогательная сфера), EPSG: 3857

Примечание. Я попытался ввести данные через ArcMap, создав растровую мозаику, которая должна работать с данными.grib, но мне не повезло.

Обновление: я также пытался использовать инструмент Project Raster, но ArcGIS не нравится проекция по умолчанию, которая берется из файла.grib и выдает ошибку.

Код, который я использую:

import gdal

src_filename = r"C:\att\project\arcshare\public\disaster_response\nrt_products\smap\20150402_20150404_anom1.grib"
dst_filename = r"C:\att\project\arcshare\public\disaster_response\nrt_products\smap\smap_py_test1.tif"

#Open existing dataset
src_ds = gdal.Open(src_filename)

#Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

#Output to new format
dst_file = driver.CreateCopy( dst_filename, src_ds, 0 )

#Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

Я не очень хорошо разбираюсь в использовании GDAL или GDAL в Python, поэтому любая помощь или советы будут с благодарностью.

Answer 1

Попробуйте использовать gdal.Translate (в Python) или gdal_translate (из командной строки). Вот два примера того, как я использовал каждый подход в прошлом:

Вариант 1: подход Python

from osgeo import gdal

# Open existing dataset
src_ds = gdal.Open(src_filename)

# Ensure number of bands in GeoTiff will be same as in GRIB file. 
bands = [] # Set up array for gdal.Translate(). 
if src_ds is not None:
    bandNum = src_ds.RasterCount # Get band count
for i in range(bandNum+1): # Update array based on band count
    if (i==0): #gdal starts band counts at 1, not 0 like the Python for loop does.
        pass
    else:
        bands.append(i)

# Open output format driver
out_form= "GTiff"

# Output to new format using gdal.Translate. See https://gdal.org/python/ for osgeo.gdal.Translate options.
dst_ds = gdal.Translate(dst_filename, src_ds, format=out_form, bandList=bands)

# Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

Вариант 2: подход с использованием командной строки gdal_translate (вызывается из Python)

import os

# Open output format driver, see gdal_translate --formats for list
out_form = "GTiff"

# Pull out specific band of interest
band=3

# Convert from GRIB to GeoTIFF using system gdal_translate
src_ds = src_filename
dst_ds = dst_filename
os.system("gdal_translate -b %s -of %s %s %s" %(str(band), out_form, src_ds, dst_ds))

В прошлом у меня были проблемы с созданием многодиапазонного GeoTiff с использованием варианта 2, поэтому я рекомендую использовать вариант 1, когда это возможно.

Answer 2

Нечто подобное должно преобразовать ваши родные координаты в желаемую проекцию. Это еще не проверено. (Может широтой вместо широт).

from cfgrib import xarray_store
from pyproj import Proj, transform

grib_data = xarray_store.open_dataset('your_grib_file.grib')
lat = grib_data.latitudes.value
lon = grib_data.longitudes.value
lon_transformed, lat_transformed = transform (Proj(init='init_projection'), 
 Proj(init='target_projection', lon, lat)