- 积分
- 284
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-6-15
- 最后登录
- 1970-1-1
|
发表于 2023-3-3 10:15:05
|
显示全部楼层
这是一个插入到等经纬度网格的代码。
#-------------读辅助数据---------------------------
dim = 2748
data = np.fromfile(rawfile,dtype=float,count=dim*dim*2)
latlon = np.reshape(data,(dim, dim, 2))
lat = latlon[:,:,0]
lon = latlon[:,:,1]
#--Define partial domain of latitude and longitude for interpolation -----------------------
MaxLon = 150.0
MinLon = 70.0
MaxLat = 50.0
MinLat = 10.0
Resolution = 0.05
Center_Lon = (MaxLon+ MinLon)/2.0
#利用np.argmax()函数获取最大值所在索引,amax()获取最大值的值。
#LowLeftIJ = np.argmin(Distance_LowLeft ) return 1-dimensional array
#Indices of the minimum elements of a N-dimensional array:
Distance_LowLeft = abs(lat -MinLat) + abs(lon -MinLon) ######2748*2748 纬向差+经向差
value1 = np.amin(Distance_LowLeft) #####最小值
LowLeftIJ = np.unravel_index(np.argmin(Distance_LowLeft, axis=None), Distance_LowLeft.shape) ######最小值得索引 (967,1248) 距离左下边界最近的一个点
Distance_UpRight = abs(lon-MaxLon) + abs(lat-MaxLat)
UpRightIJ = np.unravel_index(np.argmin(Distance_UpRight, axis=None), Distance_UpRight.shape) ######最小值得索引 (349,2108) 距离右上边界最近的一个点
Distance_UpLeft = abs(lon-MinLon) + abs(lat-MaxLat)
UpLeftIJ = np.unravel_index(np.argmin(Distance_UpLeft, axis=None), Distance_UpLeft.shape) ######最小值得索引 (***,***) 距离左上边界最近的一个点
Distance_UpCenter = abs(lon-Center_Lon) + abs(lat-MaxLat)
UpCenterIJ = np.unravel_index(np.argmin(Distance_UpCenter, axis=None), Distance_UpCenter.shape) ######最小值得索引 (***,***) 距离上边界中间最近的一个点
Distance_LowRight = abs(lon-MaxLon) + abs(lat-MinLat)
LowRightIJ = np.unravel_index(np.argmin(Distance_LowRight, axis=None), Distance_LowRight.shape) ######最小值得索引 (***,***) 距离右下边界最近的一个点
StartCol = min([LowLeftIJ[1], UpLeftIJ[1]]) #####开始的列
EndCol = max([LowRightIJ[1], UpRightIJ[1]]) #####结束的列
StartLine = min([UpLeftIJ[0], UpRightIJ[0], UpCenterIJ[0]]) ###开始的行
EndLine = max([LowLeftIJ[0], LowRightIJ[0]]) #####结束的行
Start = [StartLine, StartCol] #####起始行列号
Count = [EndLine - StartLine + 1,
EndCol - StartCol + 1, ] ##行列数量 count need not -1
Lon_cut = lon[StartLine:EndLine, StartCol:EndCol] #### 截取的经度
Lat_cut = lat[StartLine:EndLine, StartCol:EndCol] #### 截取的纬度
Sub_dataset1 = data2[StartLine:EndLine,StartCol:EndCol] ####截取所选范围数据
###-------------------------------------插值到等经纬度格点 ----------------
Lon_cut_1D = np.array(Lon_cut).flatten() ####将经度数据转换为1维 777988
Lat_cut_1D = np.array(Lat_cut).flatten() ####将纬度数据转换为1维
Sub_dataset1_1D = np.array(Sub_dataset1).flatten() ####将卫星数据转换为1维
points = np.array([Lon_cut_1D,Lat_cut_1D]) #####建立经度*纬度的新数组 2*777988
pointsT = points.transpose(1,0) # [2,*] 转置为 [*,2] 777988*2 0为经度,1为纬度
lat = np.arange(MinLat, MaxLat + 0.01, Resolution) # +0.01 为了包含右侧短点
lon = np.arange(MinLon, MaxLon + 0.01, Resolution)
lon_mesh, lat_mesh = np.meshgrid(lon, lat) ### 601*1001
grid_Sub_dataset1 = griddata(pointsT, Sub_dataset1_1D, (lon_mesh, lat_mesh), method='nearest')
# grid_Sub_dataset1 = np.flipud(grid_Sub_dataset1) # 行互换 ,上下翻转 ######shape 601*1001
print(grid_Sub_dataset1.shape) |
|