NC裁切
步骤如下：
- 读取nc文件：获取其中的经纬度数据；
- 判断裁切范围：判断各个经纬度坐标点是否在用户输入的wkt范围内，并生成掩膜矩阵；
- 写入变量：根据掩膜矩阵写入变量，不同维度的变量写入方式有所差异。
读取nc文件
def get_nc_content(self, input_nc, lon_parm, lat_parm, wkt_range):
def fill_nan_value(param, data):
attrs = nc_data.variables[param].ncattrs()
if "_FillValue" in attrs:
fill_value = nc_data.variables[param]._FillValue
data[data == fill_value] = np.nan
if "scale_factor" in attrs:
scale_factor = nc_data.variables[param].scale_factor
data = data * scale_factor
return data
logger.info(f"开始读取nc文件: {input_nc}")
nc_data = Dataset(input_nc, 'r')
keys = nc_data.variables.keys()
logger.debug(keys)
# 获取到NC文件中的经度
longitude = nc_data.variables[lon_parm][:].data
# 处理填充值和尺度因子
longitude = fill_nan_value(lon_parm, longitude)
# 获取到NC文件中的纬度
latitude = nc_data.variables[lat_parm][:].data
# 处理填充值和尺度因子
latitude = fill_nan_value(lat_parm, latitude)
判断裁切范围
续上述代码。读取nc文件和判断裁切范围应该拆开来写，这里合并到一起了。
# 加载wkt为多边形
polygon = wkt.loads(wkt_range)
# 初始化掩膜矩阵
matrix = None
if len(longitude.shape) == 1: # 经纬度为一维数据
length = longitude.size
# 创建掩膜矩阵
matrix = np.full((length, 1), np.nan)
logger.info(f"开始判断经纬度坐标是否在h5文件中: {length}×1")
for i in range(length):
# 根据客户输入的坐标来判断坐标是否在NC文件中
flag = NCFile.whether_contain(longitude[i], latitude[i], polygon)
if flag:
matrix[i, :] = 1
elif len(longitude.shape) == 2:
m, n = longitude.shape
# 创建掩膜矩阵
matrix = np.full((m, n), np.nan)
logger.info(f"开始判断经纬度坐标是否在h5文件中: {m}×{n}")
for i in range(0, m):
for j in range(0, n):
# 根据客户输入的坐标来判断坐标是否在h5文件中
flag = NCFile.whether_contain(longitude[i, j], latitude[i, j], polygon)
if flag:
matrix[i, j] = 1
logger.info("读取完毕")
return nc_data, matrix
写入变量
注意变量中的尺度因子、填充值和偏移量。先创建维度，再进行写入，我这里写入数据时变量会溢出，原因不明，仅作参考。
def write_nc(self, output_nc, nc_data, matrix):
logger.info(f"开始写入: {output_nc}")
new_nc_file = Dataset(output_nc, 'w', format='NETCDF4')
# 写入描述信息 global attributes
inattrs = nc_data.ncattrs()
for attr in inattrs:
new_nc_file.setncattr(attr, nc_data.getncattr(attr))
keys = nc_data.variables.keys()
# 记录有数据的行列号
t = np.where(matrix == 1)
min_hang, max_hang = min(t[0]), max(t[0])
# 创建维度
new_nc_file.createDimension("time", max_hang - min_hang + 1)
new_nc_file.createDimension("meas_ind")
logger.info("created Dimension:time，meas_ind")
if "wvf_ind" in keys: # SDR多两个维度 waveform wvf_ind
new_nc_file.createDimension('waveform', 2560)
new_nc_file.createDimension('wvf_ind')
logger.info("created Dimension:waveform，wvf_ind")
for key in keys:
new_variable = new_nc_file.createVariable(key,
nc_data.variables[key].dtype,
nc_data.variables[key].dimensions)
# 写入变量描述信息
inattrs1 = nc_data.variables[key].ncattrs()
for attr1 in inattrs1:
# if attr1 != "_FillValue":
# 写入尺度因子，显示为Nan，其他参数显示乱码，不写入
if attr1 == "scale_factor" or attr1 == "valid_min" or attr1 == "valid_max":
new_nc_file.variables[key].setncattr(attr1, nc_data.variables[key].getncattr(attr1))
# if "add_offset" in inattrs1: # 偏移量
# logger.debug("有add_offset", key)
# 注意！有部分变量写入时会出现变量溢出，导致生成的文件中变量为 -21478.3648。
# 读取源文件数据，不做任何处理，原封不动得写入也会导致变量溢出。只能手动调整变量类型
if "_FillValue" in inattrs1 and "scale_factor" in inattrs1:
# 写入一维变量和二维变量
fillvalue = nc_data.variables[key]._FillValue
scalefactor = nc_data.variables[key].scale_factor
new_variable.missing_value = fillvalue
a = nc_data.variables[key][:].data
a1 = np.where(a == fillvalue, np.nan, a) # 处理其中的nan值
new_data_clip = a1[min_hang:max_hang + 1]
# 比例偏移 例如scalefactor为0.1，实际值存入自动扩大10倍（除以0.1）,取出自动缩小10倍（乘以0.1）
write_data = np.where(np.isnan(new_data_clip), fillvalue * scalefactor, new_data_clip)
new_nc_file.variables[key][:] = write_data
elif "_FillValue" in inattrs1:
# print("只有_FillValue", key, NC_data.variables[key].dimensions, NC_data.variables[key].dtype)
fillvalue = nc_data.variables[key]._FillValue
new_variable.missing_value = fillvalue
a = nc_data.variables[key][:].data
a1 = np.where(a == fillvalue, np.nan, a) # 处理其中的nan值
new_data_clip = a1[min_hang:max_hang + 1]
write_data = np.where(np.isnan(new_data_clip), fillvalue, new_data_clip) # 置空值为_FillValue
new_nc_file.variables[key][:] = write_data
elif "scale_factor" in inattrs1: # lon lat int32类型
# print("只有scale_factor", key, NC_data.variables[key].dimensions, NC_data.variables[key].dtype)
scalefactor = nc_data.variables[key].scale_factor
new_variable.missing_value = -2147483648
a = nc_data.variables[key][:].data
new_data_clip = a[min_hang:max_hang + 1]
# int32类型，置空值为-2147483648
write_data = np.where(np.isnan(new_data_clip), -2147483648 * scalefactor, new_data_clip)
new_nc_file.variables[key][:] = write_data
else:
if key == 'time':
data_1 = nc_data.variables[key][:].data
data_2 = data_1[min_hang:max_hang + 1]
new_nc_file.variables[key][:] = data_2
elif key == "meas_ind" or key == "wvf_ind": # 坐标轴变量直接写入
new_nc_file.variables[key][:] = nc_data.variables[key][:].data
else:
# 二维变量 waveforms_20hz_ku waveforms_20hz_c
# 一维变量 time_day time_sec time_microsec int32类型 可归为一类
new_variable.missing_value = -2147483648
a = nc_data.variables[key][:].data
a1 = np.where(a == 2147483647, np.nan, a)
new_data_clip = a1[min_hang:max_hang + 1]
# int32类型，置空值为-2147483648
write_data = np.where(np.isnan(new_data_clip), -2147483648, new_data_clip)
new_nc_file.variables[key][:] = write_data
logger.info("写入完毕")
# 关闭NC文件
nc_data.close()
new_nc_file.close()
参考文章
浙公网安备 33010602011771号