在处理SMOS数据时,遇到了读取nc数据并存为tif后,影像投影无法改变,因此全球数据无法重叠。源数据的投影为'EPSG:6933',希望转为'EPSG:4326'。
解决代码。
```python
import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal, gdalconst
# *************************导入netCDF4
file_path = r'I:\\'
input_name = file_path + 'SM_SCIE_MIR_L4AGBB_20110101T000000_20211231T235959_100_001_8.nc'
out_path = r'I:\out\\'
res = 0.25 #目标数据的分辨率
top = 57 #目标数据的top坐标
geo_transform = [-179, res, 0.00, top, 0.00, -res] #投影转换
os.makedirs(out_path, exist_ok=True)
with nc.Dataset(input_name) as DS:
data = DS.variables['AGB'] #nc文件数据
data_arr1 = np.asarray(data) #转为数组
src_proj = DS.proj4text #nc原始文件投影信息
for j in range(0, 11):
year = str(2011 + j)
data_arr3 = np.where(data_arr1[j] < 0, 0, data_arr1[j])
data_arr4=np.flipud(data_arr3) #影像上下变换
# 获取数组的行数和列数
rows, cols = data_arr4.shape
# 创建一个新的Dataset对象
driver = gdal.GetDriverByName('MEM')
dataset = driver.Create('', cols, rows, 1, gdal.GDT_Float32)
# 将数组的数据写入Dataset对象
band = dataset.GetRasterBand(1)
band.WriteArray(data_arr4)
# 设置地理变换信息
dataset.SetGeoTransform(geo_transform) # 这里的参数需要设置一下,尽管这个参数在warp后被抹去。
# 设置投影信息
dataset.SetProjection(src_proj) # 源文件的投影。
option = gdal.WarpOptions(format='GTiff', srcSRS=src_proj, dstSRS='EPSG:4326', outputType=gdal.GDT_Float32, resampleAlg=gdalconst.GRA_NearestNeighbour)
out_name = out_path + '**' + year + '.tif'
result_data=gdal.Warp(out_name, dataset, options=option, xRes=res, yRes=res, creationOptions=["COMPRESS=LZW"], dstNodata=0)
result_data.SetGeoTransform(geo_transform) #warp后geo_transform信息被抹去,这里需要再次设定。
result_data.FlushCache()
del result_data
print(out_name)
'''