diff --git a/.gitignore b/.gitignore index c699223..0bc6b75 100644 --- a/.gitignore +++ b/.gitignore @@ -669,3 +669,5 @@ cython_debug/ +/WindSpeedModel/TestData +/WindFileData/output diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..30312d0 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "Manual-Labeling-Client/BaseCommonLibrary"] + path = Manual-Labeling-Client/BaseCommonLibrary + url = http://192.168.3.188:3000/EngineeringDepartment/BaseCommonLibrary.git diff --git a/Manual-Labeling-Client/BaseCommonLibrary b/Manual-Labeling-Client/BaseCommonLibrary new file mode 160000 index 0000000..1c1cbc3 --- /dev/null +++ b/Manual-Labeling-Client/BaseCommonLibrary @@ -0,0 +1 @@ +Subproject commit 1c1cbc3d7f3418f9e26d5a2892eace6a87e0c5e1 diff --git a/WindFileData/ERA5ToWindDataConverter.py b/WindFileData/ERA5ToWindDataConverter.py new file mode 100644 index 0000000..e69de29 diff --git a/WindFileData/WindDataHandler.py b/WindFileData/WindDataHandler.py new file mode 100644 index 0000000..5d45d5f --- /dev/null +++ b/WindFileData/WindDataHandler.py @@ -0,0 +1,68 @@ +class WindDataHandler: + def __init__(self, dll_path='./LAMPWindData.dll'): + self.dll_path = dll_path + self.wind_lib = ctypes.CDLL(dll_path) + self._setup_prototypes() + self.current_file_info = None + + def _setup_prototypes(self): + """设置 DLL 函数的参数和返回类型""" + # ... 函数原型设置代码见上一节，此处省略 ... + + def get_data_file_info(self, filepath): + """获取风场数据文件信息""" + # 将Python字符串转换为字节字符串 + filepath_bytes = filepath.encode('utf-8') + info = self.wind_lib.getDataFileInfo(filepath_bytes) + self.current_file_info = info + return info + + def get_wind_data_file_time_arr(self, filepath, info): + """获取时间数组""" + filepath_bytes = filepath.encode('utf-8') + if info.Num <= 0: + return [] + + # 创建足够大的数组来存储时间戳 + time_arr = (ctypes.c_int64 * info.Num)() + result = self.wind_lib.get_WindDataFileTimeArr(filepath_bytes, info, time_arr) + + if result == 0: # 假设返回0表示成功 + return list(time_arr) + else: + print(f"Error getting time array: {result}") + return [] + + def read_wind_data_file(self, filepath, info, time_id): + """读取指定时间ID的风场数据 (U/V分量)""" + filepath_bytes = filepath.encode('utf-8') + layer_size = info.Height * info.Width + + # 创建数组来存储U和V分量数据 + u_arr = (ctypes.c_double * layer_size)() + v_arr = (ctypes.c_double * layer_size)() + + result = self.wind_lib.Read_WindDataFile(filepath_bytes, info, time_id, u_arr, v_arr) + + if result == 0: + # 将ctypes数组转换为Python列表或numpy数组（如果可用） + return list(u_arr), list(v_arr) + else: + print(f"Error reading wind data: {result}") + return None, None + + def write_wind_data_file(self, filepath, info, time_id, u_arr, v_arr): + """将风场数据写入文件""" + filepath_bytes = filepath.encode('utf-8') + layer_size = info.Height * info.Width + + # 确保输入数据长度正确 + if len(u_arr) != layer_size or len(v_arr) != layer_size: + raise ValueError(f"UArr and VArr must have length {layer_size}") + + # 将Python列表转换为ctypes数组 + u_arr_ctypes = (ctypes.c_double * layer_size)(*u_arr) + v_arr_ctypes = (ctypes.c_double * layer_size)(*v_arr) + + result = self.wind_lib.Write_WindDataFile(filepath_bytes, info, time_id, u_arr_ctypes, v_arr_ctypes) + return result \ No newline at end of file diff --git a/WindSpeedModel/ImageHandle.py b/WindSpeedModel/ImageHandle.py new file mode 100644 index 0000000..27949e9 --- /dev/null +++ b/WindSpeedModel/ImageHandle.py @@ -0,0 +1,956 @@ +""" +@Project ：microproduct +@File ：ImageHandle.py +@Function ：实现对待处理SAR数据的读取、格式标准化和处理完后保存文件功能 +@Author ：LMM +@Date ：2021/10/19 14:39 +@Version ：1.0.0 +""" +import os +from PIL import Image +import time +from osgeo import gdal +from osgeo import osr +import numpy as np +from PIL import Image +import cv2 +import logging +from xml.etree.ElementTree import ElementTree, Element +import math +import shutil +import tarfile +logger = logging.getLogger("mylog") + + +class ImageHandler: + + @staticmethod + + def extract_tarfile(self,tar_file, target_dir): + with tarfile.open(tar_file, 'r') as tar_ref: + for member in tar_ref.getmembers(): + if member.isdir(): + continue + filename = os.path.basename(member.name) + source = tar_ref.extractfile(member) + target = open(os.path.join(target_dir, filename), "wb") + with source, target: + shutil.copyfileobj(source, target) + + """ + 影像读取、编辑、保存 + """ + def __init__(self): + pass + @staticmethod + def get_dataset(filename): + """ + :param filename: tif路径 + :return: 图像句柄 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + return dataset + + def get_scope(self, filename): + """ + :param filename: tif路径 + :return: 图像范围 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + im_scope = self.cal_img_scope(dataset) + del dataset + return im_scope + + @staticmethod + def get_projection(filename): + """ + :param filename: tif路径 + :return: 地图投影信息 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + im_proj = dataset.GetProjection() + del dataset + return im_proj + + @staticmethod + def get_geotransform(filename): + """ + :param filename: tif路径 + :return: 从图像坐标空间（行、列），也称为（像素、线）到地理参考坐标空间（投影或地理坐标）的仿射变换 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + geotransform = dataset.GetGeoTransform() + del dataset + return geotransform + + def get_invgeotransform(filename): + """ + :param filename: tif路径 + :return: 从地理参考坐标空间（投影或地理坐标）的到图像坐标空间（行、列 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + geotransform = dataset.GetGeoTransform() + geotransform=gdal.InvGeoTransform(geotransform) + del dataset + return geotransform + + @staticmethod + def get_bands(filename): + """ + :param filename: tif路径 + :return: 影像的波段数 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + bands = dataset.RasterCount + del dataset + return bands + + @staticmethod + def geo2lonlat(dataset, x, y): + """ + 将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定） + :param dataset: GDAL地理数据 + :param x: 投影坐标x + :param y: 投影坐标y + :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat) + """ + prosrs = osr.SpatialReference() + prosrs.ImportFromWkt(dataset.GetProjection()) + geosrs = prosrs.CloneGeogCS() + ct = osr.CoordinateTransformation(prosrs, geosrs) + coords = ct.TransformPoint(x, y) + return coords[:2] + + @staticmethod + def get_band_array(filename, num=1): + """ + :param filename: tif路径 + :param num: 波段序号 + :return: 对应波段的矩阵数据 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + bands = dataset.GetRasterBand(num) + array = bands.ReadAsArray(0, 0, bands.XSize, bands.YSize) + + # if 'int' in str(array.dtype): + # array[np.where(array == -9999)] = np.inf + # else: + # array[np.where(array < -9000.0)] = np.nan + + del dataset + return array + + @staticmethod + def write_imgArray(filename, im_data, no_data='null'): + """ + 影像保存 + :param filename: 保存的路径 + :param im_proj: + :param im_geotrans: + :param im_data: + :param no_data: 把无效值设置为 nodata + :return: + """ + + gdal_dtypes = { + 'int8': gdal.GDT_Byte, + 'unit16': gdal.GDT_UInt16, + 'int16': gdal.GDT_Int16, + 'unit32': gdal.GDT_UInt32, + 'int32': gdal.GDT_Int32, + 'float32': gdal.GDT_Float32, + 'float64': gdal.GDT_Float64, + } + if not gdal_dtypes.get(im_data.dtype.name, None) is None: + datatype = gdal_dtypes[im_data.dtype.name] + else: + datatype = gdal.GDT_Float32 + + # 判读数组维数 + if len(im_data.shape) == 3: + im_bands, im_height, im_width = im_data.shape + else: + im_bands, (im_height, im_width) = 1, im_data.shape + + # 创建文件 + if os.path.exists(os.path.split(filename)[0]) is False: + os.makedirs(os.path.split(filename)[0]) + + driver = gdal.GetDriverByName("GTiff") # 数据类型必须有，因为要计算需要多大内存空间 + dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) + + if im_bands == 1: + # outRaster.GetRasterBand(1).WriteArray(array) # 写入数组数据 + outband = dataset.GetRasterBand(1) + outband.WriteArray(im_data) + if no_data != 'null': + outband.SetNoDataValue(no_data) + outband.FlushCache() + else: + for i in range(im_bands): + outband = dataset.GetRasterBand(1 + i) + outband.WriteArray(im_data[i]) + outband.FlushCache() + # outRaster.GetRasterBand(i + 1).WriteArray(array[i]) + del dataset + + @staticmethod + def get_data(filename): + """ + :param filename: tif路径 + :return: 获取所有波段的数据 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + im_data = dataset.ReadAsArray(0, 0, im_width, im_height) + del dataset + return im_data + + @staticmethod + def get_scopes(ori_sim): + ori_sim_data = 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) + + scopes = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]] + return scopes + + @staticmethod + def get_dataset(filename): + """ + :param filename: tif路径 + :return: 获取所有波段的数据 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + return dataset + + @staticmethod + def get_all_band_array(filename): + """ + （大气延迟算法） + 将ERA-5影像所有波段存为一个数组, 波段数在第三维度 get_data（）->（37，8，8） + :param filename： 影像路径 get_all_band_array ->（8，8，37） + :return: 影像数组 + """ + dataset = gdal.Open(filename) + x_size = dataset.RasterXSize + y_size = dataset.RasterYSize + nums = dataset.RasterCount + array = np.zeros((y_size, x_size, nums), dtype=float) + if nums == 1: + bands_0 = dataset.GetRasterBand(1) + array = bands_0.ReadAsArray(0, 0, x_size, y_size) + else: + for i in range(0, nums): + bands = dataset.GetRasterBand(i+1) + arr = bands.ReadAsArray(0, 0, x_size, y_size) + array[:, :, i] = arr + return array + + @staticmethod + def get_img_width(filename): + """ + :param filename: tif路径 + :return: 影像宽度 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + width = dataset.RasterXSize + + del dataset + return width + + @staticmethod + def get_img_height(filename): + """ + :param filename: tif路径 + :return: 影像高度 + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + height = dataset.RasterYSize + del dataset + return height + + @staticmethod + def read_img(filename): + """ + 影像读取 + :param filename: + :return: + """ + gdal.AllRegister() + img_dataset = gdal.Open(filename) # 打开文件 + + if img_dataset is None: + msg = 'Could not open ' + filename + logger.error(msg) + return None, None, None + + im_proj = img_dataset.GetProjection() # 地图投影信息 + if im_proj is None: + return None, None, None + im_geotrans = img_dataset.GetGeoTransform() # 仿射矩阵 + + im_width = img_dataset.RasterXSize # 栅格矩阵的行数 + im_height = img_dataset.RasterYSize # 栅格矩阵的行数 + im_arr = img_dataset.ReadAsArray(0, 0, im_width, im_height) + del img_dataset + return im_proj, im_geotrans, im_arr + + def cal_img_scope(self, dataset): + """ + 计算影像的地理坐标范围 + 根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换） + :param dataset :GDAL地理数据 + :return: list[point_upleft, point_upright, point_downleft, point_downright] + """ + if dataset is None: + return None + + img_geotrans = dataset.GetGeoTransform() + if img_geotrans is None: + return None + + width = dataset.RasterXSize # 栅格矩阵的列数 + height = dataset.RasterYSize # 栅格矩阵的行数 + + point_upleft = self.trans_rowcol2geo(img_geotrans, 0, 0) + point_upright = self.trans_rowcol2geo(img_geotrans, width, 0) + point_downleft = self.trans_rowcol2geo(img_geotrans, 0, height) + point_downright = self.trans_rowcol2geo(img_geotrans, width, height) + + return [point_upleft, point_upright, point_downleft, point_downright] + + @staticmethod + def get_scope_ori_sim(filename): + """ + 计算影像的地理坐标范围 + 根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换） + :param dataset :GDAL地理数据 + :return: list[point_upleft, point_upright, point_downleft, point_downright] + """ + gdal.AllRegister() + dataset = gdal.Open(filename) + if dataset is None: + return None + + width = dataset.RasterXSize # 栅格矩阵的列数 + height = dataset.RasterYSize # 栅格矩阵的行数 + + band1 = dataset.GetRasterBand(1) + array1 = band1.ReadAsArray(0, 0, band1.XSize, band1.YSize) + + band2 = dataset.GetRasterBand(2) + array2 = band2.ReadAsArray(0, 0, band2.XSize, band2.YSize) + + if array1[0, 0] < array1[0, width-1]: + point_upleft = [array1[0, 0], array2[0, 0]] + point_upright = [array1[0, width-1], array2[0, width-1]] + else: + point_upright = [array1[0, 0], array2[0, 0]] + point_upleft = [array1[0, width-1], array2[0, width-1]] + + + if array1[height-1, 0] < array1[height-1, width-1]: + point_downleft = [array1[height - 1, 0], array2[height - 1, 0]] + point_downright = [array1[height - 1, width - 1], array2[height - 1, width - 1]] + else: + point_downright = [array1[height - 1, 0], array2[height - 1, 0]] + point_downleft = [array1[height - 1, width - 1], array2[height - 1, width - 1]] + + + if(array2[0, 0] < array2[height - 1, 0]): + #上下调换顺序 + tmp1 = point_upleft + point_upleft = point_downleft + point_downleft = tmp1 + + tmp2 = point_upright + point_upright = point_downright + point_downright = tmp2 + + return [point_upleft, point_upright, point_downleft, point_downright] + + + @staticmethod + def trans_rowcol2geo(img_geotrans,img_col, img_row): + """ + 据GDAL的六参数模型仿射矩阵将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换） + :param img_geotrans: 仿射矩阵 + :param img_col:图像纵坐标 + :param img_row:图像横坐标 + :return: [geo_x,geo_y] + """ + geo_x = img_geotrans[0] + img_geotrans[1] * img_col + img_geotrans[2] * img_row + geo_y = img_geotrans[3] + img_geotrans[4] * img_col + img_geotrans[5] * img_row + return [geo_x, geo_y] + + @staticmethod + def write_era_into_img(filename, im_proj, im_geotrans, im_data): + """ + 影像保存 + :param filename: + :param im_proj: + :param im_geotrans: + :param im_data: + :return: + """ + gdal_dtypes = { + 'int8': gdal.GDT_Byte, + 'unit16': gdal.GDT_UInt16, + 'int16': gdal.GDT_Int16, + 'unit32': gdal.GDT_UInt32, + 'int32': gdal.GDT_Int32, + 'float32': gdal.GDT_Float32, + 'float64': gdal.GDT_Float64, + } + if not gdal_dtypes.get(im_data.dtype.name, None) is None: + datatype = gdal_dtypes[im_data.dtype.name] + else: + datatype = gdal.GDT_Float32 + + # 判读数组维数 + if len(im_data.shape) == 3: + im_height, im_width, im_bands = im_data.shape # shape[0] 行数 + else: + im_bands, (im_height, im_width) = 1, im_data.shape + + # 创建文件 + if os.path.exists(os.path.split(filename)[0]) is False: + os.makedirs(os.path.split(filename)[0]) + + driver = gdal.GetDriverByName("GTiff") # 数据类型必须有，因为要计算需要多大内存空间 + dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) + dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数 + dataset.SetProjection(im_proj) # 写入投影 + + if im_bands == 1: + dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据 + else: + for i in range(im_bands): + dataset.GetRasterBand(i + 1).WriteArray(im_data[:, :, i]) + # dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) + del dataset + + # 写GeoTiff文件 + + @staticmethod + def lat_lon_to_pixel(raster_dataset_path, location): + """From zacharybears.com/using-python-to-translate-latlon-locations-to-pixels-on-a-geotiff/.""" + gdal.AllRegister() + raster_dataset = gdal.Open(raster_dataset_path) + if raster_dataset is None: + return None + ds = raster_dataset + gt = ds.GetGeoTransform() + srs = osr.SpatialReference() + srs.ImportFromWkt(ds.GetProjection()) + srs_lat_lon = srs.CloneGeogCS() + ct = osr.CoordinateTransformation(srs_lat_lon, srs) + new_location = [None, None] + # Change the point locations into the GeoTransform space + (new_location[1], new_location[0], holder) = ct.TransformPoint(location[1], location[0]) + # Translate the x and y coordinates into pixel values + Xp = new_location[0] + Yp = new_location[1] + dGeoTrans = gt + dTemp = dGeoTrans[1] * dGeoTrans[5] - dGeoTrans[2] * dGeoTrans[4] + Xpixel = (dGeoTrans[5] * (Xp - dGeoTrans[0]) - dGeoTrans[2] * (Yp - dGeoTrans[3])) / dTemp + Yline = (dGeoTrans[1] * (Yp - dGeoTrans[3]) - dGeoTrans[4] * (Xp - dGeoTrans[0])) / dTemp + del raster_dataset + return (Xpixel, Yline) + + @staticmethod + def write_img(filename, im_proj, im_geotrans, im_data, no_data='0'): + """ + 影像保存 + :param filename: 保存的路径 + :param im_proj: + :param im_geotrans: + :param im_data: + :param no_data: 把无效值设置为 nodata + :return: + """ + + gdal_dtypes = { + 'int8': gdal.GDT_Byte, + 'unit16': gdal.GDT_UInt16, + 'int16': gdal.GDT_Int16, + 'unit32': gdal.GDT_UInt32, + 'int32': gdal.GDT_Int32, + 'float32': gdal.GDT_Float32, + 'float64': gdal.GDT_Float64, + } + if not gdal_dtypes.get(im_data.dtype.name, None) is None: + datatype = gdal_dtypes[im_data.dtype.name] + else: + datatype = gdal.GDT_Float32 + flag = False + # 判读数组维数 + if len(im_data.shape) == 3: + im_bands, im_height, im_width = im_data.shape + flag = True + else: + im_bands, (im_height, im_width) = 1, im_data.shape + + # 创建文件 + if os.path.exists(os.path.split(filename)[0]) is False: + os.makedirs(os.path.split(filename)[0]) + + driver = gdal.GetDriverByName("GTiff") # 数据类型必须有，因为要计算需要多大内存空间 + dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) + + dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数 + + dataset.SetProjection(im_proj) # 写入投影 + + if im_bands == 1: + # outRaster.GetRasterBand(1).WriteArray(array) # 写入数组数据 + if flag: + outband = dataset.GetRasterBand(1) + outband.WriteArray(im_data[0]) + if no_data != 'null': + outband.SetNoDataValue(np.double(no_data)) + outband.FlushCache() + else: + outband = dataset.GetRasterBand(1) + outband.WriteArray(im_data) + if no_data != 'null': + outband.SetNoDataValue(np.double(no_data)) + outband.FlushCache() + else: + for i in range(im_bands): + outband = dataset.GetRasterBand(1 + i) + outband.WriteArray(im_data[i]) + if no_data != 'null': + outband.SetNoDataValue(np.double(no_data)) + outband.FlushCache() + # outRaster.GetRasterBand(i + 1).WriteArray(array[i]) + del dataset + + # 写GeoTiff文件 + + @staticmethod + def write_img_envi(filename, im_proj, im_geotrans, im_data, no_data='null'): + """ + 影像保存 + :param filename: 保存的路径 + :param im_proj: + :param im_geotrans: + :param im_data: + :param no_data: 把无效值设置为 nodata + :return: + """ + + gdal_dtypes = { + 'int8': gdal.GDT_Byte, + 'unit16': gdal.GDT_UInt16, + 'int16': gdal.GDT_Int16, + 'unit32': gdal.GDT_UInt32, + 'int32': gdal.GDT_Int32, + 'float32': gdal.GDT_Float32, + 'float64': gdal.GDT_Float64, + } + if not gdal_dtypes.get(im_data.dtype.name, None) is None: + datatype = gdal_dtypes[im_data.dtype.name] + else: + datatype = gdal.GDT_Float32 + + # 判读数组维数 + if len(im_data.shape) == 3: + im_bands, im_height, im_width = im_data.shape + else: + im_bands, (im_height, im_width) = 1, im_data.shape + + # 创建文件 + if os.path.exists(os.path.split(filename)[0]) is False: + os.makedirs(os.path.split(filename)[0]) + + driver = gdal.GetDriverByName("ENVI") # 数据类型必须有，因为要计算需要多大内存空间 + dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) + + dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数 + + dataset.SetProjection(im_proj) # 写入投影 + + if im_bands == 1: + # outRaster.GetRasterBand(1).WriteArray(array) # 写入数组数据 + outband = dataset.GetRasterBand(1) + outband.WriteArray(im_data) + if no_data != 'null': + outband.SetNoDataValue(no_data) + outband.FlushCache() + else: + for i in range(im_bands): + outband = dataset.GetRasterBand(1 + i) + outband.WriteArray(im_data[i]) + outband.FlushCache() + # outRaster.GetRasterBand(i + 1).WriteArray(array[i]) + del dataset + + @staticmethod + def write_img_rpc(filename, im_proj, im_geotrans, im_data, rpc_dict): + """ + 图像中写入rpc信息 + """ + # 判断栅格数据的数据类型 + if 'int8' in im_data.dtype.name: + datatype = gdal.GDT_Byte + elif 'int16' in im_data.dtype.name: + datatype = gdal.GDT_Int16 + else: + datatype = gdal.GDT_Float32 + + # 判读数组维数 + if len(im_data.shape) == 3: + im_bands, im_height, im_width = im_data.shape + else: + im_bands, (im_height, im_width) = 1, im_data.shape + + # 创建文件 + driver = gdal.GetDriverByName("GTiff") + dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) + + dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数 + dataset.SetProjection(im_proj) # 写入投影 + + # 写入RPC参数 + for k in rpc_dict.keys(): + dataset.SetMetadataItem(k, rpc_dict[k], 'RPC') + + if im_bands == 1: + dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据 + else: + for i in range(im_bands): + dataset.GetRasterBand(i + 1).WriteArray(im_data[i]) + + del dataset + + + def transtif2mask(self,out_tif_path, in_tif_path, threshold): + """ + :param out_tif_path:输出路径 + :param in_tif_path:输入的路径 + :param threshold:阈值 + """ + im_proj, im_geotrans, im_arr, im_scope = self.read_img(in_tif_path) + im_arr_mask = (im_arr < threshold).astype(int) + self.write_img(out_tif_path, im_proj, im_geotrans, im_arr_mask) + + def write_quick_view(self, tif_path, color_img=False, quick_view_path=None): + """ + 生成快视图,默认快视图和影像同路径且同名 + :param tif_path:影像路径 + :param color_img:是否生成随机伪彩色图 + :param quick_view_path:快视图路径 + """ + if quick_view_path is None: + quick_view_path = os.path.splitext(tif_path)[0]+'.jpg' + + n = self.get_bands(tif_path) + if n == 1: # 单波段 + t_data = self.get_data(tif_path) + else: # 多波段，转为强度数据 + t_data = self.get_data(tif_path) + t_data = t_data.astype(float) + t_data = np.sqrt(t_data[0] ** 2 + t_data[1] ** 2) + t_data[np.isnan(t_data)] = 0 + t_data[np.where(t_data == -9999)] = 0 + t_r = self.get_img_height(tif_path) + t_c = self.get_img_width(tif_path) + if t_r > 10000 or t_c > 10000: + q_r = int(t_r / 10) + q_c = int(t_c / 10) + elif 1024 < t_r < 10000 or 1024 < t_c < 10000: + if t_r > t_c: + q_r = 1024 + q_c = int(t_c/t_r * 1024) + else: + q_c = 1024 + q_r = int(t_r/t_c * 1024) + else: + q_r = t_r + q_c = t_c + + if color_img is True: + # 生成伪彩色图 + img = np.zeros((t_r, t_c, 3), dtype=np.uint8) # (高，宽，维度) + u = np.unique(t_data) + for i in u: + if i != 0: + w = np.where(t_data == i) + img[w[0], w[1], 0] = np.random.randint(0, 255) # 随机生成一个0到255之间的整数 可以通过挑参数设定不同的颜色范围 + img[w[0], w[1], 1] = np.random.randint(0, 255) + img[w[0], w[1], 2] = np.random.randint(0, 255) + + img = cv2.resize(img, (q_c, q_r)) # (宽，高) + cv2.imwrite(quick_view_path, img) + # cv2.imshow("result4", img) + # cv2.waitKey(0) + else: + # 灰度图 + min = np.percentile(t_data, 2) # np.nanmin(t_data) + max = np.percentile(t_data, 98) # np.nanmax(t_data) + # if (max - min) < 256: + t_data = (t_data - min) / (max - min) * 255 + out_img = Image.fromarray(t_data) + out_img = out_img.resize((q_c, q_r)) # 重采样 + out_img = out_img.convert("L") # 转换成灰度图 + out_img.save(quick_view_path) + + def get_center_scopes(self, dataset): + if dataset is None: + return None + + img_geotrans = dataset.GetGeoTransform() + if img_geotrans is None: + return None + + width = dataset.RasterXSize # 栅格矩阵的列数 + height = dataset.RasterYSize # 栅格矩阵的行数 + + x_split = int(width/5) + y_split = int(height/5) + img_col_start = x_split * 1 + img_col_end = x_split * 3 + img_row_start = y_split * 1 + img_row_end = y_split *3 + cols = img_col_end - img_col_start + rows = img_row_end - img_row_start + if cols > 10000 or rows > 10000: + img_col_end = img_col_start + 10000 + img_row_end = img_row_start + 10000 + + point_upleft = self.trans_rowcol2geo(img_geotrans, img_col_start, img_row_start) + point_upright = self.trans_rowcol2geo(img_geotrans, img_col_end, img_row_start) + point_downleft = self.trans_rowcol2geo(img_geotrans, img_col_start, img_row_end) + point_downright = self.trans_rowcol2geo(img_geotrans, img_col_end, img_row_end) + + return [point_upleft, point_upright, point_downleft, point_downright] + def write_view(self, tif_path, color_img=False, quick_view_path=None): + """ + 生成快视图,默认快视图和影像同路径且同名 + :param tif_path:影像路径 + :param color_img:是否生成随机伪彩色图 + :param quick_view_path:快视图路径 + """ + if quick_view_path is None: + quick_view_path = os.path.splitext(tif_path)[0]+'.jpg' + + n = self.get_bands(tif_path) + if n == 1: # 单波段 + t_data = self.get_data(tif_path) + else: # 多波段，转为强度数据 + t_data = self.get_data(tif_path) + t_data = t_data.astype(float) + t_data = np.sqrt(t_data[0] ** 2 + t_data[1] ** 2) + t_data[np.isnan(t_data)] = 0 + t_data[np.where(t_data == -9999)] = 0 + t_r = self.get_img_height(tif_path) + t_c = self.get_img_width(tif_path) + q_r = t_r + q_c = t_c + + if color_img is True: + # 生成伪彩色图 + img = np.zeros((t_r, t_c, 3), dtype=np.uint8) # (高，宽，维度) + u = np.unique(t_data) + for i in u: + if i != 0: + w = np.where(t_data == i) + img[w[0], w[1], 0] = np.random.randint(0, 255) # 随机生成一个0到255之间的整数 可以通过挑参数设定不同的颜色范围 + img[w[0], w[1], 1] = np.random.randint(0, 255) + img[w[0], w[1], 2] = np.random.randint(0, 255) + + img = cv2.resize(img, (q_c, q_r)) # (宽，高) + cv2.imwrite(quick_view_path, img) + # cv2.imshow("result4", img) + # cv2.waitKey(0) + else: + # 灰度图 + min = np.percentile(t_data, 2) # np.nanmin(t_data) + max = np.percentile(t_data, 98) # np.nanmax(t_data) + # if (max - min) < 256: + t_data = (t_data - min) / (max - min) * 255 + out_img = Image.fromarray(t_data) + out_img = out_img.resize((q_c, q_r)) # 重采样 + out_img = out_img.convert("L") # 转换成灰度图 + out_img.save(quick_view_path) + + return quick_view_path + + def limit_field(self, out_path, in_path, min_value, max_value): + """ + :param out_path:输出路径 + :param in_path:主mask路径，输出影像采用主mask的地理信息 + :param min_value + :param max_value + """ + proj = self.get_projection(in_path) + geotrans = self.get_geotransform(in_path) + array = self.get_band_array(in_path, 1) + array[array < min_value] = min_value + array[array > max_value] = max_value + self.write_img(out_path, proj, geotrans, array) + return True + + def band_merge(self, lon, lat, ori_sim): + lon_arr = self.get_data(lon) + lat_arr = self.get_data(lat) + temp = np.zeros((2, lon_arr.shape[0], lon_arr.shape[1]), dtype=float) + temp[0, :, :] = lon_arr[:, :] + temp[1, :, :] = lat_arr[:, :] + self.write_img(ori_sim, '', [0.0, 1.0, 0.0, 0.0, 0.0, 1.0], temp, '0') + + + def get_scopes(self, ori_sim): + ori_sim_data = self.get_data(ori_sim) + lon = ori_sim_data[0, :, :] + lat = ori_sim_data[1, :, :] + + min_lon = np.nanmin(np.where((lon != 0) & ~np.isnan(lon), lon, np.inf)) + max_lon = np.nanmax(np.where((lon != 0) & ~np.isnan(lon), lon, -np.inf)) + min_lat = np.nanmin(np.where((lat != 0) & ~np.isnan(lat), lat, np.inf)) + max_lat = np.nanmax(np.where((lat != 0) & ~np.isnan(lat), lat, -np.inf)) + + scopes = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]] + return scopes + + @staticmethod + def dem_merged(in_dem_path, out_dem_path): + ''' + DEM重采样函数，默认坐标系为WGS84 + agrs: + in_dem_path: 输入的DEM文件夹路径 + meta_file_path: 输入的xml元文件路径 + out_dem_path: 输出的DEM文件夹路径 + ''' + # 读取文件夹中所有的DEM + dem_file_paths = [os.path.join(in_dem_path, dem_name) for dem_name in os.listdir(in_dem_path) if + dem_name.find(".tif") >= 0 and dem_name.find(".tif.") == -1] + spatialreference = osr.SpatialReference() + spatialreference.SetWellKnownGeogCS("WGS84") # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree + spatialproj = spatialreference.ExportToWkt() # 导出投影结果 + # 将DEM拼接成一张大图 + mergeFile = gdal.BuildVRT(os.path.join(out_dem_path, "mergedDEM_VRT.tif"), dem_file_paths) + out_DEM = os.path.join(out_dem_path, "mergedDEM.tif") + gdal.Warp(out_DEM, + mergeFile, + format="GTiff", + dstSRS=spatialproj, + dstNodata=-9999, + outputType=gdal.GDT_Float32) + time.sleep(3) + # gdal.CloseDir(out_DEM) + return out_DEM + + @staticmethod + def landcover_merged(in_dem_path, out_dem_path): + ''' + DEM重采样函数，默认坐标系为WGS84 + agrs: + in_dem_path: 输入的DEM文件夹路径 + meta_file_path: 输入的xml元文件路径 + out_dem_path: 输出的DEM文件夹路径 + ''' + # 读取文件夹中所有的DEM + dem_file_paths = [os.path.join(in_dem_path, dem_name) for dem_name in os.listdir(in_dem_path) if + dem_name.find(".tif") >= 0 and dem_name.find(".tif.") == -1] + spatialreference = osr.SpatialReference() + spatialreference.SetWellKnownGeogCS("WGS84") # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree + spatialproj = spatialreference.ExportToWkt() # 导出投影结果 + # 将DEM拼接成一张大图 + mergeFile = gdal.BuildVRT(os.path.join(out_dem_path, "mergedDEM_VRT.tif"), dem_file_paths) + out_DEM = os.path.join(out_dem_path, "mergedDEM.tif") + gdal.Warp(out_DEM, + mergeFile, + format="GTiff", + dstSRS=spatialproj, + dstNodata=-9999, + outputType=gdal.GDT_Byte) + time.sleep(3) + # gdal.CloseDir(out_DEM) + return out_DEM + + @staticmethod + def get_inc_angle(inc_xml, rows, cols, out_path): + tree = ElementTree() + tree.parse(inc_xml) # 影像头文件 + root = tree.getroot() + values = root.findall('incidenceValue') + angle_value = [value.text for value in values] + angle_value = np.array(angle_value) + inc_angle = np.tile(angle_value, (rows, 1)) + ImageHandler.write_img(out_path, '', [0.0, 1.0, 0.0, 0.0, 0.0, 1.0], inc_angle) + + pass + + +if __name__ == '__main__': + cols = 7086 + rows = 8064 + inc_xml = r"D:\micro\WorkSpace\Dem\Temporary\processing\product\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686-DEM.tiff" + # ImageHandler().write_quick_view(inc_xml) + # fn = r'E:\202306hb\result\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892-cal-SMC.tif' + # out = r'E:\202306hb\result\soil.tif' + # # + # im_proj, im_geotrans, im_arr = ImageHandler.read_img(fn) + # im_arr[np.where(np.isnan(im_arr))] = 0 + # # h,w = im_arr.shape + # # arr = np.random.rand(h,w)*0.4 + # ImageHandler.write_img(out, im_proj, im_geotrans, im_arr, '-1') + +# path = r'D:\BaiduNetdiskDownload\GZ\lon.rdr' +# path2 = r'D:\BaiduNetdiskDownload\GZ\lat.rdr' +# path3 = r'D:\BaiduNetdiskDownload\GZ\lon_lat.tif' +# s = ImageHandler().band_merge(path, path2, path3) +# print(s) +# pass + fileDir = r'D:\micro\WorkSpace\ortho\Temporary\VH_db' + outDir = r'D:\micro\WorkSpace\ortho\Temporary\VH_db\test' + files = os.listdir(fileDir) + for file in files: + tifFile = os.path.join(fileDir, file) + outFile = os.path.join(outDir, file) + im_proj, im_geotrans, im_arr = ImageHandler.read_img(tifFile) + im_arr[np.isnan(im_arr)] = 0 + im_arr[np.isinf(im_arr)] = 0 + ImageHandler.write_img(outFile, im_proj, im_geotrans, im_arr, '0') diff --git a/WindSpeedModel/WindSpeedInversion.py b/WindSpeedModel/WindSpeedInversion.py new file mode 100644 index 0000000..6479902 --- /dev/null +++ b/WindSpeedModel/WindSpeedInversion.py @@ -0,0 +1,9 @@ +""" +@Project ：windProject +@File ：WindSpeedInversion.py +@Function ： +@Author ：CZH +@Date ：2025/10/19 14:39 +@Version ：1.0.0 +""" + diff --git a/WindSpeedModel/WindSpeedInversionBatch.py b/WindSpeedModel/WindSpeedInversionBatch.py new file mode 100644 index 0000000..e69de29