2017-08-31 3 views
1

私はSMAPデータ衛星、とりわけ湿気と土壌の比重のために働いています。EASE-2グリッド製品SMAPから地理座標への再投影を修正するにはどうすればよいですか?

私が使用GDALのアイデアはすべてを解決し、これに似たようなコードとテスト・修正Link to first approach to download SMAP data

に発表された作り、次のとおりです、だから、

import os 
import h5py 
import numpy as np 
from osgeo import gdal, gdal_array, osr 


    # the file to download 

https://n5eil01u.ecs.nsidc.org/SMAP/SPL4SMAU.003/2017.08.01/SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5

path = "/path/to/data" 
    h5File = h5py.File(path + "SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5", 'r') 
data = h5File.get('Analysis_Data/sm_rootzone_analysis') 
lat = h5File.get("cell_lat") 
lon = h5File.get("cell_lon") 

np_data = np.array(data) 
np_lat = np.array(lat) 
np_lon = np.array(lon) 

num_cols = float(np_data.shape[1]) 
num_rows = float(np_data.shape[0]) 

xmin = np_lon.min() 
xmax = np_lon.max() 
ymin = np_lat.min() 
ymax = np_lat.max() 
xres = (xmax - xmin)/num_cols 
yres = (ymax - ymin)/num_rows 

nrows, ncols = np_data.shape 
xres = (xmax - xmin)/float(ncols) 
yres = (ymax - ymin)/float(nrows) 
geotransform = (xmin, xres, 0, ymax, 0, -xres) 

dataFileOutput = path + "sm_rootzone_analysis.tif" 
output_raster = gdal.GetDriverByName('GTiff').Create(dataFileOutput, ncols, nrows, 1, gdal.GDT_Float32) # Open the file 
output_raster.SetGeoTransform(geotransform) 
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326) 

output_raster.SetProjection(srs.ExportToWkt()) 
output_raster.GetRasterBand(1).WriteArray(np_data) # Writes my array to the raster 

del output_raster 

をこのアプローチを使用すると、その結果は投影の多くの問題を伴うグローバルマップであり、例えば以下の画像はpyによって生成される上のコード。 Map produced using the code above

正しいデータと比較するために、HEG nasaソフトウェアを使用して、同じ画像をh5から抽出しました。 Map correct

答えて

0

実際にデータがEASE2グローバルグリッドにある場合は、ジオトランスフォームの緯度/経度の座標系としてEPSG:4326を割り当てないでください。

あなたは緯度/経度は9キロでEASE2グリッドに座標変換する場合は、あなたのgeotransformのようなものでなければなりません:

geotransform = (-17367530.44516138, 9000, 0, 7314540.79258289, 0, -9000.0)

およびSRS:

srs.ImportFromEPSG(6933)

関連する問題