0

我正在尝试将 .grib 文件转换为 GeoTIFF 以在 GIS 中使用（尤其是 ArcGIS），但无法正确投影图像。我已经能够在 Python 中使用 GDAL 创建一个 GeoTIFF，它可以显示数据，但在带入 ArcGIS 时未显示在正确的位置。生成的图像如下。

将数据放入 ArcGIS 时的样子

我正在使用的数据可以从以下位置下载：https ://gimms.gsfc.nasa.gov/SMOS/SMAP/L05/

我正在尝试将数据投影到 WGS84 Web Mercator (Auxiliary Sphere), EPSG: 3857

注意：我尝试通过创建一个可以处理 .grib 数据的栅格马赛克来通过 ArcMap 引入数据，但我没有任何运气。

更新：我也尝试过使用 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

我不太精通在 Python 中使用 GDAL 或 GDAL，因此将不胜感激任何帮助或提示。

4

2 回答 2

0

尝试使用 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))

过去我在使用选项 2 创建多波段 GeoTiff 时遇到过麻烦，因此我建议尽可能使用选项 1。

于 2020-07-22T17:46:40.810 回答
-1

像这样的东西应该将您的原生坐标转换为您想要的投影。这还没有经过测试。（可以按纬度而不是纬度）。

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)
于 2019-08-04T13:16:18.007 回答