汇总格网数据的时候报错坐标超出界限

import numpy as np
from netCDF4 import Dataset
import pickle
import glob
import pandas as pd
import geopandas as gpd

def bring(file_name,var_name,yr0):
    x=Dataset(file_name)
    lon=x.variables['lon'][:]
    lat=x.variables['lat'][:]
    lon_min, lon_max, lat_min, lat_max = 96, 122, 31, 42
    lon_idx = np.where((lon >= lon_min) & (lon <= lon_max))[0]
    lat_idx = np.where((lat >= lat_min) & (lat <= lat_max))[0]
    lon = lon[lon_idx]
    lat = lat[lat_idx]
    var = x.variables[var_name][:,lat_idx,lon_idx]
    #convert time to helpful datetime index
    dates=np.array(str(yr0) + '-01-01',dtype=np.datetime64)
    dts=np.linspace(0,len(x.variables['time'][:])-1,len(x.variables['time'][:]))
    dates=dates+dts.astype(int)
    dates=pd.DatetimeIndex(dates)
    return(lon,lat,var,dates)

folder='D:\\shuju\\'
stem='China_1km_prep_*'
files=glob.glob(folder+stem)

#number of days exceeding these thresholds per year (wet days at lowest)
for f in range(len(files)):
    file_name=files[f]
    var_name='prep'
    yr0 = int(file_name.split('prep_')[1][0:4])
    [lon,lat,pr,dates]=bring(file_name,var_name,yr0)
    print(pr.ndim)
    print(pr.shape)

    #calculate number of wet days, and amount of rainon wet days, and amount of rain in excess of threshold, on a sequence of absolute thresholds

    intermed=np.sum(pr,axis=0)
    intermed = np.expand_dims(intermed, axis=0)
    #append these wet days measures to the previous years values            
    if f==0:
        Pt=intermed
    elif f==1:
        Pt=np.stack((Pt,intermed),axis=0)
    else:
        Pt=np.concatenate((Pt,np.expand_dims(intermed,axis=0)),axis=0)
    del intermed
    del pr
    print('done year ' + str(int(file_name.split('prep_')[1][0:4])))

mask = gpd.read_file('D:\\shuju\\shp\\liujingdexian\\LYDX.gpkg',layer=0)

isos=pd.unique(mask.OBJECTID_1)

for i in range(len(isos)):
    measure = np.zeros((1, 1))
    # load region specific coordinates and weights
    crds_file = str(18) + '_crds_wghts.pkl'
    if len(glob.glob(crds_file)) != 0:
        with open(crds_file, 'rb') as f:
            [loni_c, lati_c, weights_c] = pickle.load(f)
            Z = Pt[:, lati_c, loni_c]
        # apply coordinates and weightings to obtain region specific measures
        if len(weights_c) != 0:
            measure[0, :] = np.dot(Z, np.array(weights_c)) / sum(weights_c)
        else:
            measure[0, :] = np.nan
        #save measures
        np.save('secondary_data/measures/precip/Pt_' + str(i) + '.npy',measure)

该回答通过自己思路及引用到GPTᴼᴾᴱᴺᴬᴵ搜索,得到内容具体如下。
在进行格网数据汇总时，报错坐标超出界限，这通常是因为在选择范围时，选择了超出数据范围的坐标。需要检查选择范围是否正确，并根据数据的实际范围进行调整。

具体来说，可以在以下代码中检查选择范围：

lon_min, lon_max, lat_min, lat_max = 96, 122, 31, 42
lon_idx = np.where((lon >= lon_min) & (lon <= lon_max))[0]
lat_idx = np.where((lat >= lat_min) & (lat <= lat_max))[0]

在这里，lon_min, lon_max, lat_min, lat_max 是选择的经度和纬度范围。如果超出了数据范围，则需要进行相应的调整。

另外，需要注意的是，在进行坐标选择时，需要保证选择的坐标是单调递增或递减的。如果选择的坐标不是单调的，可能会导致数据读取错误。

以下是修改后的完整代码：

import numpy as np
from netCDF4 import Dataset
import pickle
import glob
import pandas as pd
import geopandas as gpd

def bring(file_name,var_name,yr0):
    x=Dataset(file_name)
    lon=x.variables['lon'][:]
    lat=x.variables['lat'][:]
    lon_min, lon_max, lat_min, lat_max = 96, 122, 31, 42
    lon_idx = np.where((lon >= lon_min) & (lon <= lon_max))[0]
    lat_idx = np.where((lat >= lat_min) & (lat <= lat_max))[0]
    lon = lon[lon_idx]
    lat = lat[lat_idx]
    var = x.variables[var_name][:,lat_idx,lon_idx]
    #convert time to helpful datetime index
    dates=np.array(str(yr0) + '-01-01',dtype=np.datetime64)
    dts=np.linspace(0,len(x.variables['time'][:])-1,len(x.variables['time'][:]))
    dates=dates+dts.astype(int)
    return(lon,lat,var,dates)

folder='D:\\shuju\\'
stem='China_1km_prep_*'
files=glob.glob(folder+stem)

#number of days exceeding these thresholds per year (wet days at lowest)
Pt = np.zeros((len(files), 1, 1))
for f in range(len(files)):
    file_name=files[f]
    var_name='prep'
    yr0 = int(file_name.split('prep_')[1][0:4])
    [lon,lat,pr,dates]=bring(file_name,var_name,yr0)
    print(pr.ndim)
    print(pr.shape)

    #calculate number of wet days, and amount of rainon wet days, and amount of rain in excess of threshold, on a sequence of absolute thresholds

    intermed=np.sum(pr,axis=0)
    intermed = np.expand_dims(intermed, axis=0)
    #append these wet days measures to the previous years values            
    Pt[f, :, :] = intermed
    del intermed
    del pr
    print('done year ' + str(int(file_name.split('prep_')[1][0:4])))

mask = gpd.read_file('D:\\shuju\\shp\\liujingdexian\\LYDX.gpkg',layer=0)

isos=pd.unique(mask.OBJECTID_1)

for i in range(len(isos)):
    measure = np.zeros((1, 1))
    # load region specific coordinates and weights
    crds_file = str(18) + '_crds_wghts.pkl'
    if len(glob.glob(crds_file)) != 0:
        with open(crds_file, 'rb') as f:
            [loni_c, lati_c, weights_c] = pickle.load(f)
            Z = Pt[:, lati_c, loni_c]
        # apply coordinates and weightings to obtain region specific measures
        if len(weights_c) != 0:
            measure[0, :] = np.dot(Z, np.array(weights_c)) / sum(weights_c)
        else:
            measure[0, :] = np.nan
        #save measures
        np.save('secondary_data/measures/precip/Pt_' + str(i) + '.npy',measure)

这里对代码进行了一些修改，主要是将 Pt 数组的大小调整为 (len(files), 1, 1)，以便在后续的计算中使用。另外，也删除了一些不必要的变量和语句。

如果以上回答对您有所帮助，点击一下采纳该答案～谢谢

引用chatGPT作答，这个错误通常是由于在调用bring函数时指定的lon_min、lon_max、lat_min和lat_max超出了数据的范围。请确保您传递给bring函数的经纬度边界值在数据的有效范围内，并且与数据本身的经纬度值匹配。

如果您已经确认经纬度边界值正确，那么您可以尝试打印出经纬度值以及经纬度边界值，以便更好地理解数据的范围和问题的根源。

例如，您可以在bring函数中添加以下代码，以在调用where函数之前打印经纬度值和范围：

print(lon)
print(lat)
print(lon_min, lon_max, lat_min, lat_max)

此外，您还可以将数据的范围可视化，以更好地了解数据的分布和问题的根源。可以使用matplotlib或其他可视化库将数据可视化，或将数据导入到GIS软件中进行可视化。