Hammer坐标系转换到WGS1984
Hammer坐标系转换到WGS1984,以处理FY3D/MERSI NVI数据为例。FY3D/MERSI NVI空间分辨率250m,全球 10°×10°分幅。
01读取数据
import gdal
import h5py
infile = r"D:\微信公众号\FY3D_MERSI_5090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.HDF"
outfile = r"D:\微信公众号\FY3D_MERSI_5090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.tif"
hdf_ds = h5py.File(infile, "r")
left_x = hdf_ds.attrs['Left-Top X'][0]
left_y = hdf_ds.attrs["Left-Top Y"][0]
res_x = hdf_ds.attrs['Resolution X'][0]
res_y = hdf_ds.attrs['Resolution Y'][0]
ndviname = list(hdf_ds.keys())[6]
ndvi_ds = hdf_ds[ndviname]
rows = ndvi_ds.shape[0]
cols = ndvi_ds.shape[1]
data = ndvi_ds[()]
driver = gdal.GetDriverByName("GTiff")
outds = driver.Create(outfile, cols, rows, 1, gdal.GDT_Int16)
outds.SetGeoTransform( (float(left_x) * 1000 / 4.0, float(res_x) * 1000, 0, float(left_y) * 1000 / 4.0, 0, -1 * float(res_y) * 1000))
proj = 'PROJCS["World_Hammer",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not_specified_based_on_custom_spheroid",SPHEROID["Custom spheroid",6363961,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Hammer_Aitoff"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(data)
FY3D_MERSI_4090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS以Hammer投影坐标系写出。
FY3D_MERSI_5090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS以Hammer投影坐标系写出。
02数据拼接
import gdal
filelist = [r"D:\微信公众号\FY3D_MERSI_4090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.tif",
r"D:\微信公众号\FY3D_MERSI_5090_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.tif"]
outfile = r"D:\微信公众号\FY3D_MERSI_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.tif"
ds = gdal.Open(filelist[0])
cols = ds.RasterXSize
rows = ds.RasterYSize
ingeo = ds.GetGeoTransform()
proj = ds.GetProjection()
minx = ingeo[0]
maxy = ingeo[3]
maxx = ingeo[0] + ingeo[1] * cols
miny = ingeo[3] + ingeo[5] * rows
ds = None
for file in filelist[1:]:
ds = gdal.Open(file)
cols = ds.RasterXSize
rows = ds.RasterYSize
geo = ds.GetGeoTransform()
minx_ = geo[0]
maxy_ = geo[3]
maxx_ = geo[0] + geo[1] * cols
miny_ = geo[3] + geo[5] * rows
minx = min(minx, minx_)
maxy = max(maxy, maxy_)
maxx = max(maxx, maxx_)
miny = min(miny, miny_)
geo = None
ds = None
newcols = int((maxx - minx) / abs(ingeo[1]))
newrows = int((maxy - miny) / abs(ingeo[5]))
driver = gdal.GetDriverByName("GTiff")
outds = driver.Create(outfile, newcols, newrows, 1, gdal.GDT_Int16)
outgeo = (minx, ingeo[1], 0, maxy, 0, ingeo[5])
outds.SetGeoTransform(outgeo)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
for file in filelist:
ds = gdal.Open(file)
data = ds.ReadAsArray()
geo = ds.GetGeoTransform()
x = int(abs((geo[0] - minx) / ingeo[1]))
y = int(abs((geo[3] - maxy) / ingeo[5]))
outband.WriteArray(data, x, y)
ds = None
outband.FlushCache()
上边两幅影像的拼接结果（坐标系仍然为Hammer投影）。
03坐标系转换
import gdal
import numpy as np
import math
import osr
infile = r"D:\微信公众号\FY3D_MERSI_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS.tif"
outfile = r"D:\微信公众号\FY3D_MERSI_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS_WGS1984.tif"
ds = gdal.Open(infile)
ingeo = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
or_x = ingeo[0]
or_y = ingeo[3]
end_x = ingeo[0] + cols * ingeo[1]
end_y = ingeo[3] + rows * ingeo[5]
# X方向分块
xblocksize = int((cols + 1) / 5)
# Y方向分块
yblocksize = int((rows + 1) / 5)
lon_max = -360
lon_min = 360
lat_max = -90
lat_min = 90
for i in range(0, rows + 1, yblocksize):
if i + yblocksize < rows + 1:
numrows = yblocksize
else:
numrows = rows + 1 - i
for j in range(0, cols + 1, xblocksize):
if j + xblocksize < cols + 1:
numcols = xblocksize
else:
numcols = cols + 1 - j
# 计算所有点的Hammer坐标系下X方向坐标数组
x = ingeo[0] + j * ingeo[1]
y = ingeo[3] + i * ingeo[5]
xgrid, ygrid = np.meshgrid(np.linspace(x, x + numcols * ingeo[1], num=numcols),
np.linspace(y, y + numrows * ingeo[5], num=numrows))
# 将hammer坐标转化为经纬度坐标
# 首先将Hammer转化为-1到1
xgrid = np.where(xgrid > (18000.0 * 1000.0), (18000.0 * 1000.0) - xgrid, xgrid)
xgrid = xgrid / (18000.0 * 1000.0)
ygrid = np.where(ygrid > (9000.0 * 1000.0), (9000.0 * 1000.0) - ygrid, ygrid)
ygrid = ygrid / (9000.0 * 1000.0)
z = np.sqrt(1 - np.square(xgrid) / 2.0 - np.square(ygrid) / 2.0)
lon = 2 * np.arctan(np.sqrt(2) * xgrid * z / (2.0 * (np.square(z)) - 1))
xgrid = None
lat = np.arcsin(np.sqrt(2) * ygrid * z)
ygrid = None
z = None
lon = lon / math.pi * 180.0
lat = lat / math.pi * 180.0
lon[lon < 0] = lon[lon < 0] + 360.0
# lat[lat<0]=lat[lat<0]+180
lon_max = max(lon_max, np.max(lon))
lon_min = min(lon_min, np.min(lon))
lon = None
lat_max = max(lat_max, np.max(lat))
lat_min = min(lat_min, np.min(lat))
lat = None
newcols = math.ceil((lon_max - lon_min) / 0.0025)
newrows = math.ceil((lat_max - lat_min) / 0.0025)
driver = gdal.GetDriverByName("GTiff")
outds = driver.Create(outfile, newcols, newrows, 1, gdal.GDT_Int16)
geo2 = (lon_min, 0.0025, 0, lat_max, 0, -1 * 0.0025)
oproj_srs = osr.SpatialReference()
proj_4 = "+proj=longlat +datum=WGS84 +no_defs"
oproj_srs.ImportFromProj4(proj_4)
outds.SetGeoTransform(geo2)
outds.SetProjection(oproj_srs.ExportToWkt())
outband = outds.GetRasterBand(1)
datav = ds.ReadAsArray()
data = np.full((datav.shape[0] + 1, datav.shape[1] + 1), -32750, dtype=int)
data[0:datav.shape[0], 0:datav.shape[1]] = datav
xblocksize = int(newcols / 5)
yblocksize = int(newrows / 5)
for i in range(0, newrows, yblocksize):
if i + yblocksize < newrows:
numrows = yblocksize
else:
numrows = newrows - i
for j in range(0, newcols, xblocksize):
if j + xblocksize < newcols:
numcols = xblocksize
else:
numcols = newcols - j
x = lon_min + j * 0.0025 + 0.0025 / 2.0
y = lat_max + i * (-1 * 0.0025) - 0.0025 / 2.0
newxgrid, newygrid = np.meshgrid(np.linspace(x, x + numcols * 0.0025, num=numcols),
np.linspace(y, y + numrows * (-1 * 0.0025), num=numrows))
# 将经纬度坐标转化为Hammer坐标
newxgrid = np.where(newxgrid > 180.0, newxgrid - 360.0, newxgrid)
newxgrid = newxgrid / 180.0 * math.pi
newygrid = newygrid / 180.0 * math.pi
newz = np.sqrt(1 + np.cos(newygrid) * np.cos(newxgrid / 2.0))
x = np.cos(newygrid) * np.sin(newxgrid / 2.0) / newz
newxgrid = None
y = np.sin(newygrid) / newz
newygrid = None
newz = None
x = x * (18000.0 * 1000.0)
y = y * (9000.0 * 1000.0)
x_index = (np.floor((x - or_x) / ingeo[1])).astype(int)
x_index = np.where(x_index < 0, data.shape[1] - 1, x_index)
x_index = np.where(x_index >= data.shape[1], data.shape[1] - 1, x_index)
y_index = (np.floor((y - or_y) / ingeo[5])).astype(int)
y_index = np.where(y_index < 0, data.shape[0] - 1, y_index)
y_index = np.where(y_index >= data.shape[0], data.shape[0] - 1, y_index)
newdata = data[y_index, x_index]
outband.WriteArray(newdata, j, i)
outband.SetNoDataValue(-32750)
outband.FlushCache()
投影到WGS1984坐标系下。
04数据裁剪
import gdal
infile = r"D:\微信公众号\FY3D_MERSI_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS_WGS1984.tif"
shapefile = r"D:\我们自己的全国县级矢量\shp\province_shp\china_北京市_WGS1984.shp"
outfile = r"D:\微信公众号\FY3D_MERSI_L3_NVI_MLT_HAM_20210620_AOTD_0250M_MS_WGS1984_clip.tif"
warp_parameters = gdal.WarpOptions(format='GTiff',
cutlineDSName = shapefile ,
cropToCutline = True)
gdal.Warp(outfile, infile, options = warp_parameters)
裁剪出北京市的结果。
关注我个人wx_gzh:小Rser
