From 7dadcfb53b46c777cdac04e5692b3fc691191f2e Mon Sep 17 00:00:00 2001 From: chenzenghui <3045316072@qq.com> Date: Mon, 25 Aug 2025 10:29:15 +0800 Subject: [PATCH] =?UTF-8?q?=E5=A2=9E=E5=8A=A0=E5=A4=A9=E4=BB=AA=E5=8D=AB?= =?UTF-8?q?=E6=98=9F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- tools/DataFormatConversion/nc2bin.py | 117 ++++++++++++++++++ tools/DataFormatConversion/tiff2bin.py | 41 +++++++ tools/ImageDataOperator/resample_to_match.py | 123 +++++++++++++++++++ tools/SpacetySliceDataTools/test.py | 97 +++++++++++++++ 4 files changed, 378 insertions(+) create mode 100644 tools/DataFormatConversion/nc2bin.py create mode 100644 tools/DataFormatConversion/tiff2bin.py create mode 100644 tools/ImageDataOperator/resample_to_match.py create mode 100644 tools/SpacetySliceDataTools/test.py diff --git a/tools/DataFormatConversion/nc2bin.py b/tools/DataFormatConversion/nc2bin.py new file mode 100644 index 0000000..70352c0 --- /dev/null +++ b/tools/DataFormatConversion/nc2bin.py @@ -0,0 +1,117 @@ +# -*- coding: utf-8 -*- +import argparse +import os +from osgeo import gdal, osr +import numpy as np +import netCDF4 as nc +from datetime import datetime, timezone +import time + +def nc2ENVIBin(nc_file_path,output_envi_path,utctime): + # Input parameters (replace these with your actual values) + time_index = -1 # Example: 0 for the first time step; adjust if needed + + # Step 1: Read the NetCDF file using netCDF4 + ds=None + filename = "数据/测试文件.nc" + with open(nc_file_path, 'rb') as f: + ds = nc.Dataset('in_memory', mode='r', memory=f.read()) + print(f"Available variables: {list(ds.variables.keys())}") # Optional: Print variables for verification + valid_time_list=ds.variables["valid_time"][:] + for timeID in range(0,len(valid_time_list)): + # print(valid_time_list[timeID],utctime) + if valid_time_list[timeID] == utctime: + time_index = timeID + break + + if time_index == -1: + return -1 + else: + print(f"time_index={time_index}") + + # Extract the variable data (assuming 3D: time, lat, lon) + data_U10 = ds.variables["u10"][time_index, :, :].filled(np.nan) # Handle masked values as NaN + data_V10 = ds.variables["v10"][time_index, :, :].filled(np.nan) # Handle masked values as NaN + + # Extract lat and lon (adjust variable names if your file uses different ones, e.g., 'latitude', 'longitude') + lat = ds.variables['latitude'][:] + lon = ds.variables['longitude'][:] + + # Calculate geotransform parameters + x_res = (lon[-1] - lon[0]) / (len(lon) - 1) + y_res = (lat[0] - lat[-1]) / (len(lat) - 1) # Note: lat often decreases north to south + top_left_x = lon[0] - x_res / 2 + top_left_y = lat[0] + y_res / 2 + geotransform = (top_left_x, x_res, 0, top_left_y, 0, -y_res) # Negative y_res for north-up + + # Step 2: Create ENVI file using GDAL + driver = gdal.GetDriverByName('ENVI') + envi_ds = driver.Create( + output_envi_path, + data_U10.shape[1], # Width (lon) + data_U10.shape[0], # Height (lat) + 2, # Single band + gdal.GDT_Float32 # Assuming float data; adjust to gdal.GDT_Float64 if needed + ) + + # Set geotransform and projection (EPSG:4326 for WGS84 lat/lon) + envi_ds.SetGeoTransform(geotransform) + srs = osr.SpatialReference() + srs.ImportFromEPSG(4326) + envi_ds.SetProjection(srs.ExportToWkt()) + + # Write the data to the band + envi_band = envi_ds.GetRasterBand(1) + envi_band.SetNoDataValue(np.nan) # Set NaN as no-data value (optional) + envi_band.WriteArray(data_U10) + + envi_band = envi_ds.GetRasterBand(2) + envi_band.SetNoDataValue(np.nan) # Set NaN as no-data value (optional) + envi_band.WriteArray(data_V10) + + # Add ENVI-specific metadata to the header (optional, but recommended) + envi_ds.SetMetadataItem('wavelength units', 'None', 'ENVI') # Example; customize as needed + envi_ds.SetMetadataItem('data type', '4', 'ENVI') # 4 = float32; matches GDT_Float32 + + # Clean up + envi_ds.FlushCache() + envi_ds = None + ds.close() + + print(f"Conversion complete: {output_envi_path} (with .hdr)") + + +def getParams(): + parser = argparse.ArgumentParser() + parser.add_argument('-i','--infile',default=r'D:\风场\20250807_SAR_EC_2024年\ERA5风场\total_precipitation202410.nc', help='输入风场文件路径') + parser.add_argument('-o', '--outfile',default=r'D:\风场\20250807_SAR_EC_2024年\ERA5风场\total_precipitation202410.bin', help='输出Bin文件路径') + parser.add_argument('-t', '--time',type=int, default=2024100108,help='yyyymmddhh') + args = parser.parse_args() + return args + +def timestr2utc(timeint): + timestr = str(timeint) + + # 解析字符串为datetime对象（默认无时区信息） + dt = datetime.strptime(timestr, "%Y%m%d%H") + + # 显式标记为UTC时间（两种等效方式） + dt_utc = dt.replace(tzinfo=timezone.utc) # 方式1：使用标准库 + # 或使用pytz库：dt_utc = pytz.utc.localize(dt) + timestamp = time.mktime(dt_utc.timetuple()) + utc_timestamp = int(timestamp) # 当前UTC时间戳（秒级） + print("UTC时间:", dt_utc) + print("ISO格式:", dt_utc.isoformat()) # 输出：2024-08-01T02:00:00+00:00 + print("UTC时间戳:", utc_timestamp) # 输出：2024-08-01T02:00:00+00:00 + return utc_timestamp + + +if __name__ == '__main__': + parser = getParams() + inNcFilePath=parser.infile + outbinpath=parser.outfile + searchTime=parser.time + print('infile=',inNcFilePath) + print('outfile=',outbinpath) + time_utc=timestr2utc(searchTime) + nc2ENVIBin(inNcFilePath,outbinpath,time_utc) \ No newline at end of file diff --git a/tools/DataFormatConversion/tiff2bin.py b/tools/DataFormatConversion/tiff2bin.py new file mode 100644 index 0000000..dc5fd72 --- /dev/null +++ b/tools/DataFormatConversion/tiff2bin.py @@ -0,0 +1,41 @@ +import argparse +import os +from osgeo import gdal, osr +import numpy as np + +def tiff_to_bin_with_metadata(tiff_path, bin_path): + # 转换格式为ENVI格式的BIN（自动生成.hdr元数据文件） + options = gdal.TranslateOptions( + format='ENVI', # ENVI格式生成.hdr文件存储地理信息 + ) + # 临时输出ENVI格式文件（含.hdr） + # temp_path = bin_path.replace('.bin', '_temp.bin') + gdal.Translate(bin_path, tiff_path, options=options) + + # # 重命名文件并清理临时文件 + # os.rename(temp_path, bin_path) # 主数据文件 + # os.rename(temp_path + '.hdr', bin_path + '.hdr') # 元数据文件 + + # 可选：额外保存投影文件（.prj） + dataset = gdal.Open(tiff_path) + if dataset.GetProjection(): + with open(bin_path.replace('.bin', '.prj'), 'w') as f_prj: + f_prj.write(dataset.GetProjection()) + dataset = None + + +def getParams(): + parser = argparse.ArgumentParser() + parser.add_argument('-i','--infile',default=r'K:\风场数据\样例风场数据\4254_20240801095415\s1a-iw-grd-vv-20240801t095415-20240801t095444-055018-06b3e9-001-DB-GEO.tif', help='输入原始tiff') + parser.add_argument('-o', '--outfile',default=r'K:\风场数据\样例风场数据\4254_20240801095415\s1a-iw-grd-vv-20240801t095415-20240801t095444-055018-06b3e9-001-DB-GEO.bin', help='输出Bin文件路径') + args = parser.parse_args() + return args + +if __name__ == '__main__': + parser = getParams() + intiffPath=parser.infile + outbinPath=parser.outfile + print('infile=',intiffPath) + print('outfile=',outbinPath) + tiff_to_bin_with_metadata(intiffPath, outbinPath) + print("convert tiff to bin file successfully") \ No newline at end of file diff --git a/tools/ImageDataOperator/resample_to_match.py b/tools/ImageDataOperator/resample_to_match.py new file mode 100644 index 0000000..8a5ad99 --- /dev/null +++ b/tools/ImageDataOperator/resample_to_match.py @@ -0,0 +1,123 @@ +import argparse +import os +from osgeo import gdal, osr +import numpy as np +os.environ['GDAL_CACHEMAX'] = '2048' +def resample_to_match(image_A_path, image_B_path, output_path, resample_method=gdal.GRA_Bilinear): + """ + 将图像A重采样至与图像B相同的坐标网格和分辨率 + :param image_A_path: 待重采样的图像A路径 + :param image_B_path: 参考图像B路径 + :param output_path: 输出文件路径 + :param resample_method: 重采样方法（默认双线性插值） + """ + gdal.AllRegister() + # 1. 打开参考图像B，获取其空间参数 + ds_B = gdal.Open(image_B_path, gdal.GA_ReadOnly) + if ds_B is None: + raise RuntimeError("无法打开参考图像B") + + # 提取B的投影、仿射变换、尺寸和波段数 + proj_B = ds_B.GetProjection() + geotrans_B = ds_B.GetGeoTransform() + width_B = ds_B.RasterXSize + height_B = ds_B.RasterYSize + bands_B = ds_B.RasterCount + data_type_B = ds_B.GetRasterBand(1).DataType + + + ds_A = gdal.Open(image_A_path, gdal.GA_ReadOnly) + if ds_A is None: + raise RuntimeError("无法打开输入图像A") + + bands_A=ds_A.RasterCount + # 2. 创建输出文件，完全复制B的空间参数 + driver = gdal.GetDriverByName('ENVI') + ds_out = driver.Create( + output_path, + width_B, height_B, + bands_A, + data_type_B + ) + ds_out.SetGeoTransform(geotrans_B) + ds_out.SetProjection(proj_B) + + # 3. 执行重采样 + + # 配置WarpOptions参数（关键修改点） + warp_opts = gdal.WarpOptions( + srcSRS=ds_A.GetProjection(), # 输入投影 + dstSRS=proj_B, # 目标投影 + resampleAlg=resample_method, # 重采样方法（如gdal.GRA_Bilinear） + warpMemoryLimit=2 * 1024**3, # 单块内存限制为2GB（单位：字节） + # 分块处理相关参数 + multithread=True, # 启用多线程加速 + ) + + # 执行Warp操作（自动分块处理） + gdal.Warp( + ds_out, # 输出数据集 + ds_A, # 输入数据集 + options=warp_opts + ) + + + # 4. 释放资源 + ds_A = ds_B = ds_out = None + print(f"重采样完成，结果已保存至: {output_path}") + + + # 4. 释放资源 + ds_A = ds_B = ds_out = None + print(f"重采样完成，结果已保存至: {output_path}") + +def getParams(): + parser = argparse.ArgumentParser() + parser.add_argument('-i','--infile',default=r'F:\Data\SAR\IDL\SAR_EC_20240801\total_precipitation2024080110.bin', help='输入文件') + parser.add_argument('-o', '--outfile',default=r'F:\Data\SAR\IDL\SAR_EC_20240801\total_precipitation2024080110_resample.bin', help='输出Bin文件路径') + parser.add_argument('-r', '--reffile',default=r'F:\Data\SAR\IDL\SAR_EC_20240801\4254-DB-GEO.tif', help='参考影像') + group = parser.add_mutually_exclusive_group() + group.add_argument( + '--nearest', + action='store_const', + const='nearest', + dest='method', + help='启用最邻近插值（默认）' + ) + group.add_argument( + '--bilinear', + action='store_const', + const='bilinear', + dest='method', + help='启用双线性插值' + ) + parser.set_defaults(method='nearest') + args = parser.parse_args() + return args + +if __name__ == '__main__': + parser = getParams() + intiffPath=parser.infile + outbinPath=parser.outfile + refbinPath=parser.reffile + methodstr=parser.method + print('infile=',intiffPath) + print('outfile=',outbinPath) + print('reffile=',refbinPath) + print('method=',methodstr) + + resamplemethodgdal=None + if methodstr == 'nearest': + resamplemethodgdal=gdal.GRA_NearestNeighbour + elif methodstr == 'bilinear': + resamplemethodgdal=gdal.GRA_Bilinear + else: + resamplemethodgdal = gdal.GRA_NearestNeighbour + + resample_to_match( + image_A_path=intiffPath, + image_B_path=refbinPath, + output_path=outbinPath, + resample_method=resamplemethodgdal # 可选：gdal.GRA_NearestNeighbour（最邻近） + ) + print("resample to match successfully") \ No newline at end of file diff --git a/tools/SpacetySliceDataTools/test.py b/tools/SpacetySliceDataTools/test.py new file mode 100644 index 0000000..4db8093 --- /dev/null +++ b/tools/SpacetySliceDataTools/test.py @@ -0,0 +1,97 @@ +import rasterio +from rasterio.windows import Window +import geopandas as gpd +from shapely.geometry import box +import numpy as np +import os + +# 参数设置 +image_path = 'path/to/your/image.tif' # 替换为您的影像路径（20000x20000） +shapefile_path = 'path/to/your/shapefile.shp' # 替换为您的图斑Shapefile路径 +output_dir = 'path/to/output/tiles' # 输出切片目录 +tile_size = 1024 # 固定切片大小 +max_overlap_ratio = 0.2 # 重叠小于20% +min_step = int(tile_size * (1 - max_overlap_ratio)) # 最小步长 ≈820 + +# 创建输出目录 +os.makedirs(output_dir, exist_ok=True) + +# 读取影像 +with rasterio.open(image_path) as src: + image_width, image_height = src.width, src.height + assert image_width == 20000 and image_height == 20000, "影像大小不匹配" + transform = src.transform # 地理变换矩阵 + + # 读取图斑数据（假设CRS一致） + gdf = gpd.read_file(shapefile_path) + # 确保图斑CRS与影像一致，如果不一致需转换：gdf = gdf.to_crs(src.crs) + + # 初始化切片列表，避免重叠过多 + generated_windows = [] + + for idx, row in gdf.iterrows(): + geom = row.geometry # 获取单个图斑多边形 + + # 计算图斑的边界框 (minx, miny, maxx, maxy) - 像素坐标 + minx, miny, maxx, maxy = geom.bounds + bbox_width = maxx - minx + bbox_height = maxy - miny + + # 如果图斑太大，无法放入单个1024x1024切片，需处理（这里简单跳过，您可修改） + if bbox_width > tile_size or bbox_height > tile_size: + print(f"图斑 {idx} 太大 ({bbox_width}x{bbox_height})，跳过或需多切片处理") + continue + + # 计算覆盖完整图斑的最小窗口：居中扩展到1024x1024 + center_x = (minx + maxx) / 2 + center_y = (miny + maxy) / 2 + half_tile = tile_size / 2 + win_minx = max(0, center_x - half_tile) + win_miny = max(0, center_y - half_tile) + win_maxx = min(image_width, center_x + half_tile) + win_maxy = min(image_height, center_y + half_tile) + + # 调整以确保正好1024x1024（如果边界不足，padding或调整） + if win_maxx - win_minx < tile_size: + win_minx = max(0, win_maxx - tile_size) + if win_maxy - win_miny < tile_size: + win_miny = max(0, win_maxy - tile_size) + + proposed_window = Window(win_minx, win_miny, tile_size, tile_size) + + # 检查与现有窗口的重叠，确保<20% + too_much_overlap = False + for existing_win in generated_windows: + # 计算重叠区域 + overlap_x = max(0, min(win_minx + tile_size, existing_win.col_off + tile_size) - max(win_minx, existing_win.col_off)) + overlap_y = max(0, min(win_miny + tile_size, existing_win.row_off + tile_size) - max(win_miny, existing_win.row_off)) + overlap_area = overlap_x * overlap_y + if overlap_area > (tile_size * tile_size * max_overlap_ratio): + too_much_overlap = True + break + + if too_much_overlap: + # 如果重叠太大，调整位置（这里简单跳过，您可实现更复杂的调整，如偏移步长） + print(f"图斑 {idx} 的窗口重叠过多，跳过或调整") + continue + + # 添加到列表 + generated_windows.append(proposed_window) + + # 裁剪并保存切片 + tile_data = src.read(window=proposed_window) + tile_meta = src.meta.copy() + tile_meta.update({ + 'width': proposed_window.width, + 'height': proposed_window.height, + 'transform': rasterio.windows.transform(proposed_window, src.transform) + }) + output_path = os.path.join(output_dir, f'tile_{idx}.tif') + with rasterio.open(output_path, 'w', **tile_meta) as dst: + dst.write(tile_data) + + print(f"生成切片 {output_path}，覆盖图斑 {idx}") + + print(f"总共生成 {len(generated_windows)} 个切片") + +# 注意：这个代码确保每个切片包含至少一个完整图斑，但可能有图斑被多个切片覆盖（如果重叠允许）。如果需要覆盖所有图斑的最优 tiling，请提供更多细节。