增加天仪卫星
parent
d4b5539ce6
commit
7dadcfb53b
|
|
@ -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)
|
||||
|
|
@ -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")
|
||||
|
|
@ -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")
|
||||
|
|
@ -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，请提供更多细节。
|
||||
Loading…
Reference in New Issue