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 import math # 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, l1a_height, l1a_width): self._begin_r, self._begin_c, self._end_r, self._end_c = 0, 0, 0, 0 self.ori2geo_img = None self._mask = None self.l1a_height = l1a_height self.l1a_width = l1a_width 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 get_lonlat_points(self): lon = self.ori2geo_img[0, :, :][np.where(self._mask == 1)] lat = self.ori2geo_img[1, :, :][np.where(self._mask == 1)] data = [(row, col) for (row, col) in zip(lon, lat)] return data ###################### # 插值方法 ###################### def init_trans_para(self, sim_ori_path, roi): """裁剪L1a_img --裁剪L1A影像 --- 修改 ori_sim 变换为 sim_ori Args: src_img_path (_type_): 原始L1A影像 cuted_img_path (_type_): 待裁剪对象 roi (_type_): 裁剪roi """ ori2geo_img_height = ImageHandle.ImageHandler.get_img_height(sim_ori_path) ori2geo_img_width = ImageHandle.ImageHandler.get_img_width(sim_ori_path) ori2geo_img = ImageHandle.ImageHandler.get_data(sim_ori_path) ori2geo_gt = ImageHandle.ImageHandler.get_geotransform(sim_ori_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 # 根据 min_lon max_lon # 根据 min_lat max_lat (x_min, y_min) = ImageHandle.ImageHandler.lat_lon_to_pixel(sim_ori_path, (min_lon, min_lat)) (x_max, y_max) = ImageHandle.ImageHandler.lat_lon_to_pixel(sim_ori_path, (max_lon, max_lat)) xmin = x_min if x_min < x_max else x_max xmax = x_min if x_min > x_max else x_max ymin = y_min if y_min < y_max else y_max ymax = y_min if y_min > y_max else y_max xmin = int(math.floor(xmin)) # 列号 xmax = int(math.ceil(xmax)) # 因为python 的索引机制 # xmax = int(math.ceil(xmax)) + 1 # 因为python 的索引机制 ymin = int(math.floor(ymin)) # 行号 ymax = int(math.ceil(ymax)) # 因为pytohn的索引机制 # ymax = int(math.ceil(ymax)) + 1 # 因为pytohn的索引机制 # 处理最大最小范围 xmin = 0 if 0 > xmin else xmin ymin = 0 if 0 > ymin else ymin xmax = xmax if ori2geo_img_width > xmax else ori2geo_img_width ymax = ymax if ori2geo_img_height > ymax else ori2geo_img_height # 判断条件 xmax = xmax + 1 if xmax == xmin else xmax ymax = ymax + 1 if ymax == ymin else ymax if ymax <= ymin or xmax <= xmin or ymax > ori2geo_img_height or xmax > ori2geo_img_width or xmin < 0 or ymin < 0 or xmin > ori2geo_img_width or ymin > ori2geo_img_height or ymax < 0 or xmax < 0: msg = 'csv_roi:' + str(roi) + 'not in box,please revise csv data！' print(msg) else: r_arr = ori2geo_img[0, ymin:ymax, xmin:xmax] c_arr = ori2geo_img[1, ymin:ymax, xmin:xmax] # 构建坐标矩阵 ori2geo_mask_r_count = ymax - ymin ori2geo_mask_c_count = xmax - xmin lon_lat_arr = np.ones((2, ori2geo_mask_r_count, ori2geo_mask_c_count)) col_arr = np.arange(xmin, xmax) * np.ones((ori2geo_mask_r_count, ori2geo_mask_c_count)) row_arr = ((np.arange(ymin, ymax)) * np.ones((ori2geo_mask_c_count, ori2geo_mask_r_count))).T img_geotrans = ori2geo_gt lon_arr = img_geotrans[0] + img_geotrans[1] * col_arr + img_geotrans[2] * row_arr lat_arr = img_geotrans[3] + img_geotrans[4] * col_arr + img_geotrans[5] * row_arr lon_lat_arr[0, :, :] = lon_arr lon_lat_arr[1, :, :] = lat_arr # print("csv_roi:") # print(roi) r_min = np.floor(np.nanmin(r_arr)) # 获取 L1A 的行列号范围 r_max = np.ceil(np.nanmax(r_arr)) + 1 c_min = np.floor(np.nanmin(c_arr)) c_max = np.ceil(np.nanmax(c_arr)) + 1 # 判断是否越界 r_min = 0 if r_min < 0 else r_min r_max = self.l1a_height if r_max > self.l1a_height else r_max c_min = 0 if c_min < 0 else c_min c_max = self.l1a_width if c_max > self.l1a_width else c_max # 判断条件 r_max = r_max + 1 if r_min == r_max else r_max c_max = c_max + 1 if c_min == c_max else c_max if r_max <= r_min or c_max <= c_min or r_max > self.l1a_height or c_max > self.l1a_width or r_min < 0 or c_min < 0 or c_min > self.l1a_width or r_min > self.l1a_height or r_max < 0 or c_max < 0: msg = 'csv_roi:' + str(roi) + 'not in box,please revise csv data！' else: pass mask_geo = SAR_GEO.cut_L1A_img(lon_lat_arr, point_list) # 在地理坐标系下裁剪对应影像 mask_geo = mask_geo.reshape(-1) r_arr = r_arr.reshape(-1) c_arr = c_arr.reshape(-1) mask_geo_idx = np.where(mask_geo == 1)[0] if mask_geo_idx.shape[0] == 0: msg = 'csv_roi:' + str(roi) + 'not in box,please revise csv data！' print(msg) else: r_idx = r_arr[mask_geo_idx] c_idx = c_arr[mask_geo_idx] r_idx = r_idx - r_min # offset row c_idx = c_idx - c_min # offset col r_count = r_max - r_min # 行数 c_count = c_max - c_min # 列数 # mask_l1a = np.zeros((int(r_count), int(c_count))) * np.nan # 创建目标大小的行列号 mask = SAR_GEO.gereratorMask(r_idx.astype(np.float64), c_idx.astype(np.float64).astype(np.float64), mask_l1a) # 这个函数修改了 self._begin_r = int(r_min) self._end_r = int(r_max) self._begin_c = int(c_min) self._end_c = int(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, self._begin_c:self._end_c] 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, row, col): lonlats = [] data_roi = [] rowcols = [] measdata_list = [] data_sim = ImageHandle.ImageHandler.get_all_band_array(ori_sim_path) 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) for lonlat in lonlats: (x, y) = ImageHandle.ImageHandler.lat_lon_to_pixel(ori_sim_path, lonlat) rowcols.append([x, y]) for data, rowcol in zip(data_roi, rowcols): img_x = round(data_sim[round(rowcol[1]), round(rowcol[0]), 0]) img_y = round(data_sim[round(rowcol[1]), round(rowcol[0]), 1]) if (img_x > 0 and img_x < row and img_y > 0 and img_y < col): measdata_list.append([img_x, img_y, float(data[3])]) # 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 class TransImgL1A_ori: 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.NumpyBilinearResampler(source_def, target_def, neighbours=32, radius_of_influence=radius, epsilon=0).resample(l1a_produc) # 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'D:\micro\WorkSpace\SurfaceRoughness\Temporary\preprocessed\ori_sim_preprocessed.tif' product_tif = r'D:\micro\WorkSpace\SurfaceRoughness\Temporary\SurfaceRoughnessProduct_temp.tif' result = r'D:\micro\WorkSpace\SurfaceRoughness\Temporary\SurfaceRoughnessProduct.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, :, :] min_lon = np.nanmin(lon) max_lon = np.nanmax(lon) min_lat = np.nanmin(lat) max_lat = np.nanmax(lat) print(np.nanmin(lon)) print(np.nanmax(lon)) print(np.nanmin(lat)) print(np.nanmax(lat)) # roi_Extend = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]] roi_Extend = [[116.17328, 43.727577], [116.652504, 43.727577], [116.652504, 44.119164], [116.17328, 44.119164]] # 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) """