from tool.algorithm.transforml1a import SAR_GEO as SAR_GEO from tool.algorithm.image import ImageHandle import numpy as np import scipy from scipy.interpolate import griddata, RegularGridInterpolator import logging import pyresample as pr # 插值模块 from pyresample.bilinear import NumpyBilinearResampler from pyresample import geometry from pyresample.geometry import AreaDefinition from osgeo import osr import os # os.environ['PROJ_LIB'] = r"D:\Anaconda\envs\micro\Lib\site-packages\osgeo\data\proj" logger = logging.getLogger("mylog") ############## # 多项式回归组件 ############## # def griddata_geo(points, data, lon_grid, lat_grid, method, i, end_i): grid_data = griddata(points, data, (lon_grid, lat_grid), method=method, ) grid_data = grid_data[:, :, 0] return [i, end_i, grid_data] def griddataBlock(start_x, len_x, start_y, len_y, grid_data_input, grid_x, grid_y, method): grid_x = grid_x.reshape(-1) grid_y = grid_y.reshape(-1) grid_data_input = grid_data_input.reshape(-1) x_list = np.array(list(range(len_x))) + start_x y_list = np.array(list(range(len_y))) + start_y x_grid, y_grid = np.meshgrid(x_list, y_list) idx = np.argsort(grid_x) grid_x = grid_x[idx].reshape(-1) grid_y = grid_y[idx].reshape(-1) grid_data_input = grid_data_input[idx].reshape(-1) interp_func = RegularGridInterpolator((grid_x.reshape(-1), grid_y.reshape(-1)), grid_data_input.reshape(-1), method='slinear', bounds_error=False, fill_value=np.nan) grid_data = interp_func((x_grid, y_grid)) # grid_data = griddata(p, grid_data_input, (x_grid, y_grid), method=method) grid_data = grid_data[:, :, 0] return (x_grid, y_grid, grid_data) class polyfit2d_U: def __init__(self, x, y, z) -> None: # 定义参数 X = np.ones((x.shape[0], 10)) X[:, 0] = 1 X[:, 1] = x X[:, 2] = y X[:, 3] = x * y X[:, 4] = x ** 2 X[:, 5] = y ** 2 X[:, 6] = x * X[:, 5] X[:, 7] = y * X[:, 4] X[:, 8] = x ** 3 X[:, 9] = y ** 3 Y = z.reshape(-1, 1) A = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), Y) self.paras_fit = A def fit(self, x, y): X = np.ones((x.shape[0], 10)) X[:, 0] = 1 X[:, 1] = x X[:, 2] = y X[:, 3] = x * y X[:, 4] = x ** 2 X[:, 5] = y ** 2 X[:, 6] = x * X[:, 5] X[:, 7] = y * X[:, 4] X[:, 8] = x ** 3 X[:, 9] = y ** 3 z = np.matmul(X, self.paras_fit) return np.sum(z) class TransImgL1A: def __init__(self, ori_sim_path, roi): self._begin_r, self._begin_c, self._end_r, self._end_c = 0, 0, 0, 0 self._mask = None self._min_lon, self._max_lon, self._min_lat, self._max_lat = 0, 0, 0, 0 self.init_trans_para(ori_sim_path, roi) def get_roi_points(self): rowcol = np.where(self._mask == 1) data = [(self._begin_r + row, self._begin_c + col) for (row, col) in zip(rowcol[0], rowcol[1])] return data ###################### # 插值方法 ###################### def init_trans_para(self, ori_sim_path, roi): """裁剪L1a_img Args: src_img_path (_type_): 原始L1A影像 cuted_img_path (_type_): 待裁剪对象 roi (_type_): 裁剪roi """ ori2geo_img = ImageHandle.ImageHandler.get_data(ori_sim_path) point_list = np.array(roi) min_lon = np.nanmin(point_list[:, 0]) max_lon = np.nanmax(point_list[:, 0]) min_lat = np.nanmin(point_list[:, 1]) max_lat = np.nanmax(point_list[:, 1]) self._min_lon, self._max_lon, self._min_lat, self._max_lat = min_lon, max_lon, min_lat, max_lat r_c_list = np.where( (ori2geo_img[0, :, :] >= min_lon) & (ori2geo_img[0, :, :] <= max_lon) & (ori2geo_img[1, :, :] >= min_lat) & (ori2geo_img[1, :, :] <= max_lat)) # if len(r_c_list) == 0 or r_c_list[0] == [] or r_c_list[1] == [] or np.array(r_c_list).size == 0: msg = 'csv_roi:' + str(roi) + 'not in box,please revise csv data！' print(msg) else: # print("csv_roi:") # print(roi) r_min = np.nanmin(r_c_list[0]) r_max = np.nanmax(r_c_list[0]) c_min = np.nanmin(r_c_list[1]) c_max = np.nanmax(r_c_list[1]) ori2geo_img = ori2geo_img[:, r_min:r_max + 1, c_min:c_max + 1] # 开始调用组件 计算 mask = SAR_GEO.cut_L1A_img(ori2geo_img.astype(np.float64), point_list) self._begin_r = r_min self._end_r = r_max self._begin_c = c_min self._end_c = c_max self._mask = mask def cut_L1A(self, in_path, out_path): img = ImageHandle.ImageHandler.get_data(in_path) if len(img.shape) == 3: cut_img = img[:, self._begin_r:self._end_r + 1, self._begin_c:self._end_c + 1] cut_img[0, :, :] = cut_img[0, :, :] * self._mask cut_img[1, :, :] = cut_img[1, :, :] * self._mask ImageHandle.ImageHandler.write_img(out_path, '', [0, 0, 0, 0, 0, 0], cut_img) else: cut_img = img[self._begin_r:self._end_r + 1, self._begin_c:self._end_c + 1] cut_img[:, :] = cut_img[:, :] * self._mask cut_img[:, :] = cut_img[:, :] * self._mask ImageHandle.ImageHandler.write_img(out_path, '', [0, 0, 0, 0, 0, 0], cut_img) def grid_interp_to_station(self, all_data, station_lon, station_lat, method='linear'): ''' func: 将等经纬度网格值 插值到 离散站点。使用griddata进行插值 inputs: all_data,形式为：[grid_lon,grid_lat,data] 即[经度网格，纬度网格，数值网格] station_lon: 站点经度 station_lat: 站点纬度。可以是 单个点，列表或者一维数组 method: 插值方法,默认使用 linear ''' station_lon = np.array(station_lon).reshape(-1, 1) station_lat = np.array(station_lat).reshape(-1, 1) lon = all_data[0].reshape(-1, 1) lat = all_data[1].reshape(-1, 1) data = all_data[2].reshape(-1, 1) points = np.concatenate([lon, lat], axis=1) station_value = griddata(points, data, (station_lon, station_lat), method=method) station_value = station_value[:, :, 0] return station_value ##################### # 当存在 ori2geo.tif ##################### @staticmethod def cut_L1a_img(src_img_path, cuted_img_path, roi): """裁剪L1a_img Args: src_img_path (_type_): 原始L1A影像 cuted_img_path (_type_): 待裁剪对象 roi (_type_): 裁剪roi """ ori2geo_img = ImageHandle.ImageHandler.get_data(src_img_path) point_list = np.array(roi) # 开始调用组件 计算 mask = SAR_GEO.cut_L1A_img(ori2geo_img.astype(np.float64), point_list) # ori2geo_img[0, :, :] = ori2geo_img[0, :, :] * mask ori2geo_img[1, :, :] = ori2geo_img[1, :, :] * mask ImageHandle.ImageHandler.write_img(cuted_img_path, '', [0, 0, 0, 0, 0, 0], ori2geo_img) return ori2geo_img # 保存成影像 def tran_geo_to_l1a(self, geo_img_path, out_l1a_img_path, ori_sim_img_path, is_class=False): """裁剪后的有投影信息的影像(cover、ndvi)转换到L1A裁剪影像的尺寸 Args: geo_img_path (_type_): _description_ out_l1a_img_path (_type_): _description_ ori_sim_img_path (_type_): _description_ geo_img_path:地理影像路径 out_l1a_img_path：转换L1A坐标系图像路径 ori_sim_img_path：裁剪后模拟影像路径 is_clss: 是否是 定性类产品 """ inverse_gt = ImageHandle.ImageHandler.get_invgeotransform(geo_img_path) ori2geo_tif = ImageHandle.ImageHandler.get_data(ori_sim_img_path) height = ImageHandle.ImageHandler.get_img_height(geo_img_path) width = ImageHandle.ImageHandler.get_img_width(geo_img_path) # 计算投影 x = ori2geo_tif[0, :, :] # lon lat x,y y = ori2geo_tif[1, :, :] ori2geo_tif[0, :, :] = inverse_gt[0] + inverse_gt[1] * x + inverse_gt[2] * y # x ori2geo_tif[1, :, :] = inverse_gt[3] + inverse_gt[4] * x + inverse_gt[5] * y # y del x, y geo_tif = ImageHandle.ImageHandler.get_data(geo_img_path) # 获取目标影像 ori2geo_tif_shape = ori2geo_tif.shape # height,width if is_class: ori2geo_tif = np.round(ori2geo_tif).astype(np.int32) mask = (ori2geo_tif[0, :, :] >= 0) & (ori2geo_tif[0, :, :] < width) & (ori2geo_tif[1, :, :] >= 0) & ( ori2geo_tif[1, :, :] < height) ori2geo_tif[0, :, :] = ori2geo_tif[0, :, :] * mask ori2geo_tif[1, :, :] = ori2geo_tif[1, :, :] * mask geo_tif_shape = geo_tif.shape geo_tif_l1a = geo_tif[ori2geo_tif[1, :, :].reshape(-1), ori2geo_tif[0, :, :].reshape(-1)].reshape( ori2geo_tif.shape[1], ori2geo_tif.shape[2]).astype(np.float32) del ori2geo_tif, geo_tif one_ids = np.where(mask == False) geo_tif_l1a[one_ids[0], one_ids[1]] = np.nan ImageHandle.ImageHandler.write_img(out_l1a_img_path, '', [0, 0, 0, 0, 0, 0], geo_tif_l1a) # save_temp_L1A(out_l1a_img_path,geo_tif_l1a) return geo_tif_l1a else: # 数值性插值 mask = (ori2geo_tif[0, :, :] > 0) & (ori2geo_tif[0, :, :] < width - 1) & (ori2geo_tif[1, :, :] > 0) & ( ori2geo_tif[1, :, :] < height - 1) one_ids = np.where(mask == 1) x, y = np.meshgrid(np.arange(0, width), np.arange(0, height)) result_data = self.grid_interp_to_station([y.reshape(-1), x.reshape(-1), geo_tif.reshape(-1)], ori2geo_tif[1, one_ids[0], one_ids[1]].reshape(-1), ori2geo_tif[0, one_ids[0], one_ids[1]].reshape(-1), method='linear').reshape(-1) mask = mask.reshape(-1) result_data_result = np.zeros((ori2geo_tif.shape[1], ori2geo_tif.shape[2])) result_data_result[:, :] = np.nan result_data_result = SAR_GEO.insert_data(result_data_result, one_ids[0].astype(np.int32), one_ids[1].astype(np.int32), result_data) ImageHandle.ImageHandler.write_img(out_l1a_img_path, '', [0, 0, 0, 0, 0, 0], result_data_result) # save_temp_L1A(out_l1a_img_path,result_data_result) return result_data_result def tran_lonlats_to_rowcols(self, lonlats, ori_sim_img_path): """ 功能：输入经纬度坐标，输出图像行列号 函数名称：tran_lonlats_to_rowcols(lonlats,out_rowcols,ori_sim_img_path) Lonlats:经纬度坐标，示例[[120.53, 31.5], [120.61, 31.5], [120.53, 31.45], [120.61, 31.45]] out_rowcols:图像行列号，示例[[0, 0], [7000, 0], [7000, 8000], [0, 8000]] ori_sim_img_path：裁剪后模拟影像路径 """ ori2geo_tif = ImageHandle.ImageHandler.get_data(ori_sim_img_path) min_lon = np.nanmin(ori2geo_tif[0, :, :]) max_lon = np.nanmax(ori2geo_tif[0, :, :]) min_lat = np.nanmin(ori2geo_tif[1, :, :]) max_lat = np.nanmax(ori2geo_tif[1, :, :]) result = [] for i in range(len(lonlats)): p = lonlats[i] if min_lon > p[0] or max_lon < p[0] or min_lat > p[1] or max_lat < p[1]: result.append([-1, -1]) continue temp_x = np.square(ori2geo_tif[0, :, :] - p[0]) + np.square(ori2geo_tif[1, :, :] - p[1]) r_c_list = [] r_c = np.argmin(temp_x) r_c = [r_c // temp_x.shape[1], r_c % temp_x.shape[1]] r_c_list.append([r_c[0], r_c[1], ori2geo_tif[0, r_c[0], r_c[1]], ori2geo_tif[1, r_c[0], r_c[1]]]) # 插值 for i in range(r_c[0] - 3, r_c[0] + 3): if i < 0 or i > temp_x.shape[0] - 1: continue for j in range(r_c[1] - 3, r_c[1] + 3): if j < 0 or j > temp_x.shape[1] - 1: continue r_c_list.append([i, j, ori2geo_tif[0, i, j], ori2geo_tif[1, i, j]]) r_c_list = np.array(r_c_list) points = r_c_list[:, 2:] f_r = scipy.interpolate.interp2d(r_c_list[:, 2], r_c_list[:, 3], r_c_list[:, 0], kind='linear') f_c = scipy.interpolate.interp2d(r_c_list[:, 2], r_c_list[:, 3], r_c_list[:, 1], kind='linear') tar_get_r = f_r(p[0], p[1])[0] tar_get_c = f_c(p[0], p[1])[0] if tar_get_r < ori2geo_tif.shape[1] and tar_get_c < ori2geo_tif.shape[ 2] and tar_get_r >= 0 and tar_get_c >= 0: lon_temp = ori2geo_tif[0, int(round(tar_get_r)), int(round(tar_get_c))] lon_lat = ori2geo_tif[1, int(round(tar_get_r)), int(round(tar_get_c))] # 增加条件筛选 result.append([tar_get_r, tar_get_c]) else: result.append([-1, -1]) return result def tran_lonlats_to_L1A_rowcols(self, meas_data, ori_sim_path): lonlats = [] data_roi = [] for data in meas_data: lon = float(data[1]) lat = float(data[2]) if (lon > self._min_lon and lon < self._max_lon and lat > self._min_lat and lat < self._max_lat): lonlats.append([lon, lat]) data_roi.append(data) rowcols = self.tran_lonlats_to_rowcols(lonlats, ori_sim_path) measdata_list = [] for data, rowcol in zip(data_roi, rowcols): if (rowcol[0] != -1 and rowcol[1] != -1): measdata_list.append( [round(rowcol[0]) - self._begin_r, round(rowcol[1]) - self._begin_c, float(data[3])]) return measdata_list @staticmethod def get_radius_of_influence(lalo_step, src_meta='radar2geo', ratio=3): EARTH_RADIUS = 6378122.65 # m """Get radius of influence based on the lookup table resolution in lat/lon direction""" if src_meta == "geo2radar": # geo2radar radius = 100e3 else: # radar2geo step_deg = max(np.abs(lalo_step)) step_m = step_deg * np.pi / 180.0 * EARTH_RADIUS radius = step_m * ratio return radius def interp2d_station_to_grid(self, lon, lat, data, loc_range=[18, 54, 73, 135], det_grid=1, method='linear', projCode=4326): # 参考链接 https://blog.csdn.net/weixin_43718675/article/details/103497930 ''' func : 将站点数据插值到等经纬度格点 inputs: lon: 站点的经度 lat: 站点的纬度 data: 对应经纬度站点的 气象要素值 loc_range: [lat_min,lat_max,lon_min,lon_max]。站点数据插值到loc_range这个范围 det_grid: 插值形成的网格空间分辨率 method: 所选插值方法，默认 0.125 return: [lon_grid,lat_grid,data_grid] ''' # step1: 先将 lon,lat,data转换成 n*1 的array数组 lon = np.array(lon).reshape(-1, 1) lat = np.array(lat).reshape(-1, 1) data = np.array(data).reshape(-1, 1) # step2:确定插值区域的经纬度网格 lat_min = loc_range[0] # y lat_max = loc_range[1] # y lon_min = loc_range[2] # x lon_max = loc_range[3] # x gt = [0, 0, 0, 0, 0, 0] gt[0] = lon_min # x gt[1] = det_grid gt[3] = lat_max # y gt[5] = -det_grid lat_count = int((lat_max - lat_min) / det_grid + 1) # y lon_count = int((lon_max - lon_min) / det_grid + 1) # x # 替换为pyresample 插值方法 proj_osr = osr.SpatialReference() proj_osr.ImportFromEPSG(projCode) projection = proj_osr.ExportToPROJJSON() # （lower_left_x、lower_left_y、upper_right_x、upper_right_y） target_def = AreaDefinition("id1", "WGS84", "proj_id", projection, lon_count, lat_count, [lon_min, lat_min, lon_max, lat_max]) source_def = geometry.SwathDefinition(lons=lon, lats=lat) lalo_step = [det_grid, -det_grid] radius = TransImgL1A.get_radius_of_influence(lalo_step, src_meta='radar2geo') result = pr.bilinear.resample_bilinear(data, source_def, target_def, radius=radius, neighbours=32, nprocs=8, fill_value=np.nan, epsilon=0) # return result def geocoding(self, ori_geo_tif, produc_arr, pixel_delta=1, method='linear'): # 参考链接 https://blog.csdn.net/weixin_43718675/article/details/103497930 ori_geo_tif[np.isnan(ori_geo_tif)] = -1 lon_data = ori_geo_tif[0, :, :].reshape(-1) lat_data = ori_geo_tif[1, :, :].reshape(-1) idx = np.where(lat_data != -1) lat_data = lat_data[idx] lon_data = lon_data[idx] idx = np.where(lon_data != -1) lat_data = lat_data[idx] lon_data = lon_data[idx] # ########################################### result = self.interp2d_station_to_grid(lon_data, lat_data, produc_arr, [self._min_lat, self._max_lat, self._min_lon, self._max_lon], det_grid=pixel_delta, method=method) return result def l1a_2_geo(self, ori_geo_path, l1a_produc_path, geo_produc_path, method='linear'): ori_geo_tif = ImageHandle.ImageHandler.get_data(ori_geo_path) # l1a_produc = ImageHandle.ImageHandler.get_data(l1a_produc_path) l1a_produc = ImageHandle.ImageHandler.get_band_array(l1a_produc_path, 1) pixel_delta_y = (self._max_lat - self._min_lat) / (self._end_r - self._begin_r) # 0.001 pixel_delta_x = (self._max_lon - self._min_lon) / (self._end_c - self._begin_c) lon_data = ori_geo_tif[0, :, :].reshape(-1) lat_data = ori_geo_tif[1, :, :].reshape(-1) l1a_produc = l1a_produc.reshape(-1) idx = np.logical_not(np.isnan(lon_data)) lat_data = lat_data[idx] lon_data = lon_data[idx] l1a_produc = l1a_produc[idx] idx = np.logical_not(np.isnan(lat_data)) lat_data = lat_data[idx] lon_data = lon_data[idx] l1a_produc = l1a_produc[idx] gt = [self._min_lon, pixel_delta_x, 0.0, self._max_lat, 0.0, -pixel_delta_y] [lat_min, lat_max, lon_min, lon_max] = [self._min_lat, self._max_lat, self._min_lon, self._max_lon] lat_count = int((lat_max - lat_min) / pixel_delta_y + 1) # y lon_count = int((lon_max - lon_min) / pixel_delta_x + 1) # x # 获取地理坐标系统信息，用于选取需要的地理坐标系统 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84" proj = srs.ExportToWkt() projection = srs.ExportToPROJJSON() # （lower_left_x、lower_left_y、upper_right_x、upper_right_y） target_def = AreaDefinition("id1", "WGS84", "proj_id", projection, lon_count, lat_count, [lon_min, lat_min, lon_max, lat_max]) lon_data = lon_data.reshape(-1, 1) lat_data = lat_data.reshape(-1, 1) l1a_produc = l1a_produc.reshape(-1, 1) source_def = geometry.SwathDefinition(lons=lon_data, lats=lat_data) lalo_step = [pixel_delta_x, -pixel_delta_y] radius = TransImgL1A.get_radius_of_influence(lalo_step, src_meta='radar2geo') geo_produc = pr.bilinear.resample_bilinear(l1a_produc, source_def, target_def, radius=radius, neighbours=32, nprocs=8, fill_value=np.nan, epsilon=0) ImageHandle.ImageHandler.write_img(geo_produc_path, proj, gt, geo_produc, np.nan) def l1a_2_geo_int(self, ori_geo_path, l1a_produc_path, geo_produc_path, method='nearest'): ori_geo_tif = ImageHandle.ImageHandler.get_data(ori_geo_path) # l1a_produc = ImageHandle.ImageHandler.get_data(l1a_produc_path) l1a_produc = ImageHandle.ImageHandler.get_band_array(l1a_produc_path, 1) pixel_delta_y = (self._max_lat - self._min_lat) / (self._end_r - self._begin_r) # 0.001 pixel_delta_x = (self._max_lon - self._min_lon) / (self._end_c - self._begin_c) lon_data = ori_geo_tif[0, :, :].reshape(-1) lat_data = ori_geo_tif[1, :, :].reshape(-1) l1a_produc = l1a_produc.reshape(-1) idx = np.logical_not(np.isnan(lon_data)) lat_data = lat_data[idx] lon_data = lon_data[idx] l1a_produc = l1a_produc[idx] idx = np.logical_not(np.isnan(lat_data)) lat_data = lat_data[idx] lon_data = lon_data[idx] l1a_produc = l1a_produc[idx] gt = [self._min_lon, pixel_delta_x, 0.0, self._max_lat, 0.0, -pixel_delta_y] [lat_min, lat_max, lon_min, lon_max] = [self._min_lat, self._max_lat, self._min_lon, self._max_lon] lat_count = int((lat_max - lat_min) / pixel_delta_y + 1) # y lon_count = int((lon_max - lon_min) / pixel_delta_x + 1) # x # 获取地理坐标系统信息，用于选取需要的地理坐标系统 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84" proj = srs.ExportToWkt() projection = srs.ExportToPROJJSON() # （lower_left_x、lower_left_y、upper_right_x、upper_right_y） target_def = AreaDefinition("id1", "WGS84", "proj_id", projection, lon_count, lat_count, [lon_min, lat_min, lon_max, lat_max]) lon_data = lon_data.reshape(-1, 1) lat_data = lat_data.reshape(-1, 1) l1a_produc = l1a_produc.reshape(-1, 1) source_def = geometry.SwathDefinition(lons=lon_data, lats=lat_data) lalo_step = [pixel_delta_x, -pixel_delta_y] radius = TransImgL1A.get_radius_of_influence(lalo_step, src_meta='radar2geo') if method == 'linear': geo_produc = pr.bilinear.resample_bilinear(l1a_produc, source_def, target_def, radius=radius, neighbours=32, nprocs=8, fill_value=0, epsilon=0) elif method == 'nearest': geo_produc = pr.kd_tree.resample_nearest(source_def, l1a_produc, target_def, epsilon=0, radius_of_influence=50000, fill_value=0, nprocs=8 ) geo_produc = geo_produc[:,:,0] ImageHandle.ImageHandler.write_img(geo_produc_path, proj, gt, geo_produc) @property def mask(self): return self._mask if __name__ == '__main__': # ori_sim_path = r"I:\坐标转换\坐标转换接口\L1A数据（l1a_img_path数据）\RPC_ori_sim.tif" # roi_Extend = [[120.53, 31.5], [120.61, 31.5], [120.61, 31.45], [120.53, 31.45]] # conver_path = r"I:\坐标转换\坐标转换接口\裁剪后辅助数据（geo_img_path数据）\Covering_cut.tif" # ndvi_path = r"I:\坐标转换\坐标转换接口\裁剪后辅助数据（geo_img_path数据）\NDVI_cut.tif" # out_path = r"I:\坐标转换\SAR2GEO\test" # # tr = TransImgL1A(ori_sim_path,roi_Extend) # tr.l1a_2_geo("I:/cut.tif", "I:/salinity.tif", "I:/salinity_geo2.tif") ori_sim = r'F:\20230605\ortho\RD_ori_sim.tif' product_tif = r'D:\micro\WorkSpace\soil_test\proces\oh2004tif\oh2004_mv.tif' result = r'D:\micro\WorkSpace\soil_test\proces\oh2004tif\soil_geo.tif' method = 'linear' """ 31.14;31.50;120.34;120.75 """ # roi_Extend = [[102.12, 33.879], [102.327, 33.879], [102.327, 33.66], [102.12, 31.45]] ori_sim_data = ImageHandle.ImageHandler.get_data(ori_sim) lon = ori_sim_data[0,:,:] lat = ori_sim_data[1,:,:] print(np.nanmin(lon)) print(np.nanmax(lon)) print(np.nanmin(lat)) print(np.nanmax(lat)) roi_Extend = [[121.4627, 41.255276], [121.91531, 41.255276], [121.91531, 40.86639], [121.4627, 40.86639]] # roi_Extend = [[108.51960117899473, 38.192443138079895], [109.62308480328566, 38.192443138079895], [109.62308480328566, 37.69300142375064], [108.51960117899473, 37.69300142375064]] tr = TransImgL1A(ori_sim, roi_Extend) tr.l1a_2_geo_int(ori_sim, product_tif, result, method) pass """ import numpy as np from pyresample import kd_tree, geometry area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', {'a': '6378144.0', 'b': '6356759.0', 'lat_0': '50.00', 'lat_ts': '50.00', 'lon_0': '8.00', 'proj': 'stere'}, 800, 800, [-1370912.72, -909968.64, 1029087.28, 1490031.36]) data = np.fromfunction(lambda y, x: y*x, (50, 10)) lons = np.fromfunction(lambda y, x: 3 + x, (50, 10)) lats = np.fromfunction(lambda y, x: 75 - y, (50, 10)) swath_def = geometry.SwathDefinition(lons=lons, lats=lats) result = kd_tree.resample_nearest(swath_def, data,area_def, radius_of_influence=50000, epsilon=0.5) """