问题:
按照区域范围(min_lat, max_lat)~(min_lon, max_lon)读取一个NC文件,经、纬数据切片如何获取维数是一个大问题,切片数量不当,容易导致微小的误差,在后续的矩阵运算带来大的灾难。
例如:依据小间距网格计算结果,需与大间距网格数据进行误差比较时,一般会对小间距风格数据进行插值,此时因为插值后的网格精度是确定的,所以需要根据区域范围,计算出插值后的经纬向维数值1,利用维数值进行遍历插值。
而根据区域范围直接读取大间距网格数据,需要确定切片取值的经纬向维数值2,1与2的值需要严格一致,否则会造成矩阵运算错误,提示计算矩阵维数不一致错误信息。
过程:
1、经纬向维数值1计算公式
纬向维数:lat_steps = int(round((max_lat - min_lat)/grid_size,1)) +1
经向维数:lon_steps = int(round((max_lon - min_lon)/grid_size,1)) +1
此处假设经纬向网格间距一致,均为grid_size值
采向round取小数1位,防止浮点计算误差,采用int方式取整,亦可用math.ceil,math.floor取整,方式不影响结果
经向计算维数存在一些复杂性,因为涉及到跨越日经线或零度线的问题。
如经向值域为(-180,180],则需考虑跨越日经线,如经向范围为:170~-160,30度的范围,
此时:max_lon = -160,min_lon = 170,显然:min_lon > max_lon。
则经向维数:lon_steps = int(round((max_lon + 360 - min_lon)/grid_size,1)) + 1
经向值域为(0, 360],则需考虑跨越零度经线(格林威治线),计算方式一致。
2、准备总维数:
lats = len(data.variables['lat']
lons = len(data.variables['lon']
2、切片取值方法
切片取值基本方法:[起始序号值:起始序号值+切片维数值]
需分别计算起始序号值与切片维数值。
纬向计算:
st_lat_index = max(math.floor((min_lat - nc_min_lat)/grid_size),0) # st_lat为区域起始点,nc_min_lat为文件的经向起始点
lat_steps = min(int(round((max_lat - min_lat)/grid_size,1)) +1,lats) # 与之前的计算保持一致
rdata = data.variables[feature][..., start_lat_index:start_lat_index+lat_steps, ...]
经向计算:
如果不考虑跨越日经线,计算方式完全一致,但跨越日经线带来诸多麻烦事。
st_lon_index = max(math.floor((st_lon - min_lon)/grid_size),0) #起始序号,与纬向一致
if min_lon > max_lon: #跨越日经线,策略:分段取值,再拼接,故需分别求出 #两切片维数值,同时确保两段维数总和与之前计算的经向维数值一致。
lon_steps = min(int(round((max_lon + 360 - min_lon)/grid_size,1)) + 1, lons)
lon_steps1 = lons - start_lon_index #用经向总维数 - 起始序号
lon_steps2 = lon_steps - lon_steps1
rdata1 = data.variables[feature][..., ..., start_lon_index:]
rdata2 = data.variables[feature][..., ..., : lon_steps2]
rdata=np.concatenate((rdata1,rdata2), axis=2) # axis = 1/2根据实际情况定
else: #简单模式
lon_steps = int(round((max_lon-min_lon)/grid_size,1)) + 1
rdata = data.variables[feature][..., ..., start_lon_index:start_lon_index+lon_steps]