NC-子数组的最大累加和

最优子数组和:动规算法详解与实现
本文解析了如何使用动态规划解决求解数组中最大连续子数组和的问题,通过代码实例展示了如何更新数组元素并找到最大和。适合理解数组问题和动规应用的开发者。

题目:
在这里插入图片描述

解决:动规

    public int maxsumofSubarray (int[] arr) {
        // write code here
        if(arr.length<1)
            return -1;
        int max = arr[0];
        for(int i =1;i<arr.length;i++){
            arr[i] = arr[i]>arr[i]+arr[i-1]?arr[i]:arr[i]+arr[i-1];
            max = Math.max(max,arr[i]);
        }
        return max;
    }
我通过以下python代码,分析NC文件,并绘制了停留时间的后向轨迹图。我还需要根据山东省16地市的SHP文件(16个单独的shp文件)的面积区域,针对每个释放点,分别计算停留时间在山东省16地市区域内的在整个NC文件区域的比重,将这个比重数据计算出来,并存储在csv文件中,请帮我编写这段代码 原代码如下: def main_analysis(nc_file, output_dir): # ===================== 参数配置 ===================== # nc_file = "grid_time_20250921000000.nc" # 替换为您的NC文件路径 # output_dir = "./lpdm_plots2/" # 输出图片目录 height_layer = 0 # 选择高度层(0为地面层,共4层) dpi = 300 # 图片分辨率 font() # ===================== 数据读取 ===================== ds = nc.Dataset(nc_file) # 提取基础变量 lons = ds.variables['longitude'][:] # (200,) lats = ds.variables['latitude'][:] # (200,) lon_grid, lat_grid = np.meshgrid(lons, lats) # 经纬度网格 (200x200) lons_min = lon_grid.min() lons_max = lon_grid.max() lats_min = lat_grid.min() lats_max = lat_grid.max() # extent = [95, 130, 22, 46] # 经度范围95-130E,纬度范围22-46N extent = [lons_min, lons_max, lats_min, lats_max] # 提取释放点名称 (RELCOM) relcom = ds.variables['RELCOM'][:] # shape: (14, 45) rel_names = [relcom[i].tobytes().decode('utf-8').strip().replace('city', '') for i in range(14)] # ------------------- 关键修改:从time变量获取起始时间 ------------------- # 读取time变量(假设为一维时间数组,单位为"hours since ..."或类似) time_var = ds.variables['time'] if time_var is None: raise ValueError("NC文件中未找到'time'变量,请检查文件结构") # 获取第一个时间点的数值(通常是自基准时间的偏移量) start_time_num = time_var[0] end_time_num = time_var[-1] # 将数值时间转换为datetime对象(关键步骤) # 注意:需使用time变量的units属性(如"hours since 2025-09-21 00:00:00") start_time = nc.num2date( start_time_num, units=time_var.units, calendar=getattr(time_var, 'calendar', 'gregorian') # 默认使用公历 ) end_time = nc.num2date( end_time_num, units=time_var.units, calendar=getattr(time_var, 'calendar', 'gregorian') # 默认使用公历 ) total_days = (start_time - end_time).days + 1 # ----------------------------------------------------------------------- # 提取停留时间数据 (1 ageclass, 14 pointspec, 216 time, 4 height, 200 lat, 200 lon) spec_data = ds.variables['spec001_mr'][:] # 提取释放点信息 rel_lons = (ds.variables['RELLNG1'][:] + ds.variables['RELLNG2'][:]) / 2 rel_lats = (ds.variables['RELLAT1'][:] + ds.variables['RELLAT2'][:]) / 2 # ===================== 绘图设置 ===================== # 地图投影设置(参考图片范围) proj = ccrs.PlateCarree() bounds = [ 0, 0.00001,0.0005, 0.001, 0.005, 0.01, 0.1, 0.3, 0.5, 0.8, 1, 2, 3 # 10 是上界,处理 >3 的数据 ] # 2. 定义对应颜色(手动调色,严格匹配图片中的颜色分布) colors = [ '#FFFFFF', '#8F8FFF', # 0.0005-0.001 浅紫 '#3D3DFF', # 0.001-0.005 稍深紫 '#010FEB', # 0.005-0.01 蓝 '#03499D', # 0.01-0.1 浅蓝(接近图片深蓝) '#06844F', # 0.1-0.3 绿 '#66D009', # 0.5-0.8 黄 '#B5E805', # 0.8-1 金黄 '#FFFF00', # 1-2 橙 '#FFA500', # 2-3 深橙 '#FF0000', # 3-10 红(覆盖 >3 区间) '#9C0202' ] # 3. 创建自定义配色方案和分级规范 cmap = ListedColormap(colors) norm = BoundaryNorm(bounds, cmap.N) # ===================== 绘制每个释放点和每个天 ===================== num_pointspec = spec_data.shape[1] # 14个释放点 for point_idx in range(num_pointspec): logger.info(f"正在绘制{pinyin_to_hanzi[rel_names[point_idx]]}的图片") # 此种方式为按天累加绘制 start_idx = 0 for day in range(total_days): # 9天(24h/day) logger.info(f"{pinyin_to_hanzi[rel_names[point_idx]]}第{day+1}张图片绘制中...") # 计算时间切片 # start_idx = day * 24 # 此种方式为按天分别绘制 end_idx = (day + 1) * 24 # 提取当前点位和当前天的停留时间数据 # shape: (24 time, 4 height, 200 lat, 200 lon) daily_data = spec_data[0, point_idx, start_idx:end_idx, :, :, :] # 固定高度层并积分时间维度 daily_cumulative = daily_data[:, height_layer, :, :].sum(axis=0) / 3600.0 # 转换为小时 # ------------------- 新增:生成透明度掩码(alpha数组) ------------------- # 设定透明度阈值(小于0.0005时完全透明) alpha_threshold = 0.00001 # 生成与daily_cumulative形状相同的alpha数组(0=透明,1=不透明) alpha = np.where(daily_cumulative < alpha_threshold, 0, 1) # ----------------------------------------------------------------------- # 创建地图画布 fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=proj) ax.set_extent(extent, crs=proj) # 添加地理要素 # ------------------ 原地图地理要素 ------------------ ax.add_feature(cfeature.COASTLINE, linewidth=0.8) # 海岸线 ax.add_feature(cfeature.BORDERS, linestyle=":", linewidth=0.5) # 国界 ax.add_feature(cfeature.OCEAN, facecolor="lightblue") # 海洋 zorderidx = 1 plot_china_province(ds, ax, 'lon', 'lat', 'fig', zorderidx) zorderidx += 1 plot_china_citys_withoutname('', ax, '', '', fig, zorderidx) zorderidx += 1 plot_gridline('', ax, 'lon', 'lat', fig, zorderidx) zorderidx += 1 # 绘制停留时间分布 im = ax.pcolormesh(lon_grid, lat_grid, daily_cumulative, cmap=cmap, norm=norm, transform=proj, alpha=alpha) # 添加颜色条 ax_height = ax.get_position().height cbar = plt.colorbar(im, ax=ax, aspect=40, # 控制长宽比 shrink=ax_height * 0.98, # 长度 fraction=0.03, pad=0.04, ticks=bounds) cbar.set_label('停留时间(小时)', size=10, labelpad=1) cbar.ax.set_yticklabels([f'{x:.5g}' for x in bounds]) # 标记释放点位置 # ax.plot(rel_lons[point_idx], rel_lats[point_idx], 'ro', markersize=8, # transform=proj, label='Release Point') # 创建自定义五角星路径 star_vertices = [ (0, 1), (-0.2245, 0.3090), (-0.9511, 0.3090), (-0.3633, -0.1180), (-0.5878, -0.8090), (0, -0.3819), (0.5878, -0.8090), (0.3633, -0.1180), (0.9511, 0.3090), (0.2245, 0.3090), (0, 1) ] star_codes = [PlotPath.MOVETO] + [PlotPath.LINETO] * 9 + [PlotPath.CLOSEPOLY] star_path = PlotPath(star_vertices, star_codes) # 使用自定义路径 plt.scatter( rel_lons[point_idx], rel_lats[point_idx], color='none', s=100, # 需要更大尺寸 marker=star_path, edgecolor='black', zorder=5, transform=ccrs.PlateCarree() ) # 调整地图边界宽度 for spine in ax.spines.values(): spine.set_linewidth(1.0) # 确保x轴和y轴单位长度一致 ax.set_aspect('equal') # -------------------- 新增:顶端样式文本 -------------------- # 左上角日期(示例:2025092016) target_time = start_time - timedelta(days=day) target_time_str = target_time.strftime("%Y%m%d%H") path_target_time_str = target_time.strftime("%Y-%m-%d") ax.text(0.02, 1.04, f'{target_time_str}', transform=ax.transAxes, fontsize=14, color='black', va='top', ha='left') # 中间地点(示例:hanzhong) ax.text(0.5, 1.04, f'{pinyin_to_hanzi[rel_names[point_idx]]}', transform=ax.transAxes, fontsize=14, color='black', va='top', ha='center') # 右上角模型信息(示例:4day-LPDM) ax.text(0.98, 1.04, f'{day + 1:02d}DAY-LPDM', transform=ax.transAxes, fontsize=14, color='black', va='top', ha='right') savepath = output_dir / f'{path_target_time_str}' / f'{pinyin_to_code[rel_names[point_idx]]}' if not os.path.exists(savepath): os.makedirs(savepath) # 保存图片 plt.savefig(savepath / f"lpdm_{pinyin_to_code[rel_names[point_idx]]}_{path_target_time_str}.png", dpi=dpi, bbox_inches='tight') plt.close(fig) logger.info(f"完成!共生成 {num_pointspec * 9} 张图片,保存至 {output_dir} 目录") ds.close()
最新发布
10-21
<think>我们需要计算每个释放点(14个)的停留时间在山东省16个地市的区域内的比重(即每个地市区域内停留时间总和占整个nc文件区域停留时间总和的比重),并将结果存入CSV文件。 步骤: 1. 读取山东省16地市的SHP文件(每个地市一个单独的shp文件)。假设这些shp文件存放在同一个目录下,我们可以使用geopandas读取。 2. 对于每个释放点,我们都需要处理其停留时间数据(按天累加后的数据,即每个点每天一张图对应的数据,但注意原代码是按天分别绘制,这里我们需要按天分别计算比重吗?还是按整个时间段?根据问题,应该是整个NC文件区域,但问题描述中并没有明确是按天还是累计。根据需求“停留时间在山东省16地市区域内的在整个NC文件区域的比重”,这里应该是整个时间段(即9天)的累计停留时间?还是按天分别计算?由于问题没有明确,这里我们假设需要按天分别计算(因为原代码是按天分别绘图),并且每个释放点每天计算16个地市的比重。如果不需要按天,而是整个时间段,那么只需要在时间维度上累加即可。但问题要求的是“整个NC文件区域”,注意整个NC文件区域是指我们的网格区域,而地市区域是其中的一部分。 根据问题描述,我们需要计算的是:针对每个释放点,分别计算停留时间在山东省16地市区域内的在整个NC文件区域的比重。这里没有提到天,所以可能是整个时间段(9天)的累计。但原代码中每天都会生成一张图,并且每天的数据是独立的(按24小时累加)。因此,我们需要明确: - 方案1:按天计算,即每天每个释放点在每个地市的比重,那么结果将是14个释放点×9天×16地市的数据- 方案2:整个时间段(9天)的累积,那么每个释放点只计算16个地市的数据(14×16)。 由于问题描述中并没有明确,但是要求存储在csv文件中,我们可以考虑两种方案: - 如果要求每天的数据,那么csv文件可能包含:释放点名称、日期、地市名称、停留时间在该地市的比重(相对于整个区域)。 - 如果要求整个时间段,那么csv文件包含:释放点名称、地市名称、比重。 考虑到原代码是按天绘制,且停留时间是按天累加的,我推测可能需要按天计算。但问题描述中说“整个NC文件区域”,这里的整个区域可能是指整个网格区域(即200x200网格),而且没有提到时间维度。为了明确,我们再看问题:“针对每个释放点,分别计算停留时间在山东省16地市区域内的在整个NC文件区域的比重”。这里没有提时间,所以可能是整个时间段(即9天)的累计。 因此,我决定按照整个时间段(9天)的累计来计算。即对于每个释放点,我们将其9天(216小时)的停留时间数据在高度层0上累加(注意原代码中只用了高度层0),然后计算这个累计数据在每个地市区域内的总和,再除以整个网格区域的总和(即整个NC文件区域的总和),得到比重。 但是注意:原代码中每天的停留时间数据是独立的,我们也可以按天计算再累加,效果一样(因为累加是线性的)。所以我们可以先计算每个释放点的总停留时间(9天×高度0)。 具体步骤: 步骤1:读取所有地市的shp文件,合并为一个GeoDataFrame(每个地市一个多边形,并带有地市名称)。注意:每个shp文件代表一个地市,我们假设文件名就是地市名称(不带扩展名),或者shp文件中有一个字段代表地名(如'name')。由于问题没有说明,我们假设每个shp文件都有一个字段可以表示地名(例如,字段名为'name')。如果文件名就是地名,那么我们可以通过文件名获取地名。 步骤2:对于每个释放点(point_idx从0到13): a. 提取该释放点在整个时间段(0到215小时)的停留时间数据,并选择高度层0,然后沿时间轴累加(求和),得到网格数据(200x200)。 b. 计算整个网格区域的总停留时间:total = np.sum(grid_data) 注意:这里没有考虑网格面积权重?问题要求的是比重,即分子是地市内停留时间,分母是整个网格区域停留时间。所以不需要面积权重(因为都是数值累加)。 c. 对于每个地市: - 使用地市的多边形创建一个掩膜(mask),该掩膜是一个200x200的布尔数组,表示网格点是否在该地市多边形内。 - 注意:我们的网格是经纬度网格(lat_grid, lon_grid),每个网格点是一个坐标。我们需要判断每个网格点是否在地市多边形内。 - 计算该地市区域内停留时间总和:masked_data = grid_data * mask 然后求和:sum_in_city = np.sum(masked_data) - 比重 = sum_in_city / total (注意total可能为0,需要处理) 步骤3:将结果组织成DataFrame(列:释放点名称,地市名称,比重),然后保存为CSV。 但是,这里有一个问题:整个网格区域(200x200)超出了山东省的范围,但是我们的数据是整个网格区域,分母是整个网格区域的总和(包括山东省以外)。这样计算出来的比重可能很小(因为山东省只占网格的一部分)。但是问题要求的是“在整个NC文件区域的比重”,所以分母是整个网格。 另外,关于掩膜创建:我们可以使用geopandas和rasterio的geometry_mask?或者直接用geopandas的contains方法循环每个点?但这样效率低。我们可以使用rasterio的rasterize,或者用matplotlib.path.Path。 由于我们已经有了网格点坐标(lon_grid, lat_grid),我们可以使用geopandas的GeoSeries.contains方法,但是要注意,我们的网格点坐标是二维的(200x200),我们可以将其展平,然后判断每个点是否在某个地市的多边形内,然后重塑形状。 具体方法如下: city_gdf = geopandas.read_file(city_shp_path) # 读取一个地市的shp文件(注意:一个文件代表一个地市,所以每个文件只有一个要素?) # 我们假设每个文件只有一个要素(即一个地市),如果有多个,则取第一个。 polygon = city_gdf.geometry.iloc[0] # 获取该地市的多边形(注意:可能是MultiPolygon,这里假设是Polygon) # 将网格点坐标展平 points = list(zip(lon_grid.ravel(), lat_grid.ravel())) # 使用polygon.contains判断每个点(注意:GeoSeries.contains只接受GeoSeries,所以我们可以创建一个点的GeoSeries) # 但是这样循环效率低(200*200=40000个点),我们可以使用向量化方法:用shapely.prepared.prep来加速。 更高效的方法:使用rasterio的rasterize函数,但需要转换到相同投影?这里我们都是经纬度,所以可以使用。 但是,由于我们只有16个地市,每个地市只需要做一次掩膜,而且网格点只有40000个,循环16次*40000次点(64万次)在现代化计算机上是可以接受的。 另一种方法是:事先将16个地市合并为一个GeoDataFrame,然后一次性为每个地市创建掩膜(避免重复读取文件)。然后,我们循环16个地市,分别创建掩膜。 步骤: 1. 读取所有地市的shp文件,合并为一个GeoDataFrame,并记录地市名称(从文件名或者属性字段获取)。 2. 创建一组掩膜(16个,每个掩膜是200x200的布尔数组),表示每个网格点是否在该地市多边形内。 假设:16个shp文件存放在一个目录(shp_dir)下,文件名如:济南市.shp、青岛市.shp等。 代码步骤: 1. 读取所有shp文件,创建地市列表和对应的几何图形列表。 cities = [] # 存储地市名称(不带扩展名) geometries = [] # 存储每个地市的几何图形(第一个要素) for shp_file in os.listdir(shp_dir): if shp_file.endswith('.shp'): city_name = shp_file[:-4] # 去掉扩展名 gdf = gpd.read_file(os.path.join(shp_dir, shp_file)) # 假设每个文件只有一个要素 if len(gdf) > 0: geom = gdf.geometry.iloc[0] cities.append(city_name) geometries.append(geom) 然后创建一个总的GeoDataFrame: city_gdf = gpd.GeoDataFrame({'city_name': cities, 'geometry': geometries}) 2. 准备网格点坐标(200x200),并转换为点几何对象的列表(注意:这里我们不需要创建整个网格的GeoSeries,因为我们要循环每个地市,对每个地市判断网格点是否在其内部) 3. 为每个地市创建掩膜(布尔数组): 掩膜的形状是(200,200),初始化为False。 遍历city_gdf的每一行(即每个地市): 获取该地市的多边形:polygon = row['geometry'] 然后,遍历每个网格点(双重循环),判断点是否在多边形内(使用polygon.contains(Point(lon, lat))),将结果存入掩膜数组。 但是这样效率低(200*200*16=64万次),所以我们可以用向量化方法?或者使用预先准备的shapely.prepared.polygon。 优化:使用shapely.prepared创建预处理的几何体,然后使用列表推导式。 prep_polygon = prep(polygon) points = [Point(lon, lat) for (lon, lat) in zip(lon_grid.ravel(), lat_grid.ravel())] # 然后使用map或者列表推导式判断 mask = np.array([prep_polygon.contains(p) for p in points]).reshape(lon_grid.shape) 但是这样仍然要循环40000次,但是预处理会加速。 或者,我们可以使用一个函数,输入一个坐标点列表,返回一个布尔数组。 4. 对于每个释放点,计算其总停留时间网格(9天累加)。 5. 对于每个地市,用掩膜提取该地市内的网格值,然后求和,再除以整个网格的总和(该释放点),得到比重。 6. 将结果存入DataFrame:释放点名称(rel_names[point_idx])、地市名称、比重。 7. 保存为CSV。 注意:整个网格区域的总停留时间(分母)是np.sum(grid_data),这个grid_data是二维的(200x200),对所有网格点求和。 由于网格点数量多,且计算掩膜需要时间,我们可以预先计算16个地市的掩膜(只需要计算一次,然后用于所有释放点)。 另外,如果网格点在地市边界外,那么就不计入。如果网格点正好落在多边形内,才计入。 代码实现: 我们将添加一个新的函数来计算掩膜数组(针对一个多边形): from shapely.prepared import prep from shapely.geometry import Point def create_mask_from_polygon(polygon, lon_grid, lat_grid): # 准备多边形 prep_polygon = prep(polygon) # 获取网格形状 ny, nx = lon_grid.shape mask = np.zeros((ny, nx), dtype=bool) # 遍历每个网格点 for i in range(ny): for j in range(nx): point = Point(lon_grid[i, j], lat_grid[i, j]) if prep_polygon.contains(point): mask[i, j] = True return mask 但是,双重循环效率低(200*200=40000次),我们可以尝试向量化?或者使用其他库(如rasterio.mask)?但是这里我们使用简单方法,因为16*40000=64万次,在普通电脑上可能需要几秒到几十秒。 另一种高效方法:使用matplotlib.path.Path(适用于多边形): from matplotlib.path import Path def create_mask_from_polygon(polygon, lon_grid, lat_grid): # 获取多边形的边界点(注意:可能包含多个部分,这里简化处理,只取第一个外边界?) # 对于多边形,我们可以获取其exterior坐标 if polygon.geom_type == 'Polygon': coords = polygon.exterior.coords[:] elif polygon.geom_type == 'MultiPolygon': # 取第一个多边形(通常一个地市是一个整体) coords = polygon.geoms[0].exterior.coords[:] else: raise ValueError("Unsupported geometry type") path = Path(coords) # 将网格展平 points = np.vstack((lon_grid.ravel(), lat_grid.ravel())).T mask = path.contains_points(points).reshape(lon_grid.shape) return mask 这种方法可能更快,因为内部是向量化操作。 我们选择使用matplotlib.path.Path。 步骤: 1. 读取所有地市的shp,并创建掩膜数组(16个掩膜,每个200x200)。 2. 计算每个释放点的总停留时间网格(9天累加,高度层0)。 3. 遍历每个释放点,然后遍历16个地市,计算该地市的停留时间总和以及整个网格总和,然后计算比重。 4. 收集每个释放点每个地市的比重数据。 注意:如果整个网格的总和(分母)为0,则比重设为0(或NaN)以避免除以0错误。 代码实现细节: 我们将主要步骤封装在main_analysis函数后面,或者单独写一个函数,然后在main_analysis中调用。 由于原代码已经关闭了nc文件,我们需要在关闭之前提取数据。所以我们把这段代码放在原代码中绘制循环之后,但在ds.close()之前。 但是注意原代码已经对每个释放点每天绘制了图片,我们不需要改动绘图部分,只需要在绘图结束后,进行这一部分计算。 因此,我们在原代码的绘图循环之后,添加计算比重的代码,然后保存csv。 具体位置:在绘制完所有图片后(即最外层循环point_idx结束之后),在ds.close()之前。 步骤代码: 1. 设置shp文件目录(需要作为参数传入?或者硬编码?)这里我们假设作为参数传入,修改函数签名: def main_analysis(nc_file, output_dir, shp_dir): 2. 读取shp文件,创建地市掩膜列表(含地市名称和掩膜数组)。 3. 计算每个释放点的总停留时间网格(9天累加)。 4. 遍历释放点和地市,计算比重,并保存结果到csv。 代码示例: 注意:以下代码需要插入原代码的适当位置(在绘制循环之后,ds.close()之前) 我们将代码封装在一个函数中,但为了清晰,直接写在main_analysis中。 考虑到代码较长,我们只写出关键部分。 注意:由于停留时间数据是spec001_mr[0, point_idx, :, height_layer, :, :],然后沿时间轴(第0轴)求和。 开始写代码: </think>以下是一个完整的解决方案,用于计算每个释放点在山东省16个地市的停留时间比重并导出为CSV文件: ```python import geopandas as gpd import numpy as np import pandas as pd import netCDF4 as nc import os from shapely.geometry import Point from tqdm import tqdm def calculate_residence_ratios(nc_file, shp_dir, output_csv): """ 计算每个释放点在山东省16地市的停留时间比重 :param nc_file: NC文件路径 :param shp_dir: 包含16地市SHP文件的目录 :param output_csv: 输出CSV文件路径 """ # 读取NC文件 ds = nc.Dataset(nc_file) # 读取经纬度网格 lons = ds.variables['longitude'][:] # (200,) lats = ds.variables['latitude'][:] # (200,) lon_grid, lat_grid = np.meshgrid(lons, lats) # (200, 200) # 获取释放点信息 relcom = ds.variables['RELCOM'][:] rel_names = [relcom[i].tobytes().decode('utf-8').strip().replace('city', '') for i in range(14)] heights = [0] # 只使用地面层高度层 # 加载停留时间数据 (1 ageclass, 14 pointspec, 216 time, 4 height, 200 lat, 200 lon) spec_data = ds.variables['spec001_mr'][:] # 准备结果数据结构 results = [] city_names = [] # 加载所有地市SHP文件 print("加载地市边界数据...") city_shapes = {} for shp_file in os.listdir(shp_dir): if shp_file.endswith('.shp'): city_name = os.path.splitext(shp_file)[0] city_names.append(city_name) gdf = gpd.read_file(os.path.join(shp_dir, shp_file)) # 确保所有几何图形在同一个CRS中 gdf = gdf.to_crs("EPSG:4326") # 转换为WGS84坐标系 city_shapes[city_name] = gdf.geometry.iloc[0] # 为每个释放点创建网格掩膜 print("创建网格掩膜...") masks = {} for city_name, geom in city_shapes.items(): mask = np.zeros_like(lon_grid, dtype=bool) for i in range(len(lats)): for j in range(len(lons)): point = Point(lon_grid[i, j], lat_grid[i, j]) if geom.contains(point): mask[i, j] = True masks[city_name] = mask # 计算每个释放点的停留时间比重 print("计算停留时间比重...") for point_idx in tqdm(range(len(rel_names))): # 累加所有时间和高度层的停留时间 total_residence = np.zeros_like(lon_grid) for height in heights: # 选择当前高度层和时间序列的所有数据 height_data = spec_data[0, point_idx, :, height, :, :] # 沿时间轴求和 time_sum = np.sum(height_data, axis=0) total_residence += time_sum # 计算整个网格的总停留时间 total_grid_sum = np.sum(total_residence) # 计算每个地市的比重 for city_name in city_names: city_mask = masks[city_name] city_sum = np.sum(total_residence[city_mask]) if total_grid_sum > 0: ratio = city_sum / total_grid_sum else: ratio = 0.0 results.append({ 'Release Point': rel_names[point_idx], 'City': city_name, 'Residence Ratio': ratio }) # 转换为DataFrame并保存CSV df = pd.DataFrame(results) df.to_csv(output_csv, index=False) print(f"结果已保存至: {output_csv}") # 关闭NC文件 ds.close() return df # 使用示例 if __name__ == "__main__": nc_file = "grid_time_20250921000000.nc" # 替换为实际NC文件路径 shp_dir = "shandong_cities_shp" # 山东16地市SHP文件目录 output_csv = "residence_ratios.csv" # 输出CSV文件路径 calculate_residence_ratios(nc_file, shp_dir, output_csv) ``` ### 代码详解: 1. **数据准备阶段**: - 读取NC文件中的经纬度网格和停留时间数据 - 解析14个释放点名称 - 加载16个地市的SHP边界文件并转换为WGS84坐标系 2. **网格掩膜创建**: - 为每个地市创建布尔掩膜矩阵(200×200) - 遍历每个网格点,判断是否位于地市边界内 - 存储每个地市的掩膜矩阵供后续计算使用 3. **停留时间计算**: - 对每个释放点: - 累加所有时间步和地面层高度的停留时间 - 计算整个网格区域的停留时间总和 - 使用各地市掩膜提取对应区域内的停留时间 - 计算比重:城市停留时间 / 网格总停留时间 4. **结果输出**: - 将结果组织为DataFrame(释放点、城市、比重三列) - 保存为CSV文件,格式示例: ``` Release Point,City,Residence Ratio jinan,Jinan,0.125 jinan,Qingdao,0.087 ... ``` ### 注意事项: 1. SHP文件命名应直接使用城市名(如"济南市.shp"),代码会使用文件名作为城市标识 2. 坐标系自动转换为WGS84(EPSG:4326)确保与NC文件坐标系一致 3. 仅使用地面层(height_layer=0)数据,如需其他高度层可修改heights参数 4. 使用tqdm库显示进度条,处理大量数据时更直观 5. 处理大型网格时可能耗时较长(200×200×16≈64万次空间判断) ### 优化建议: 1. 对于大型数据集可使用Dask进行并行计算 2. 使用更高效的空间索引技术(如R-tree) 3. 缓存掩膜矩阵避免重复计算
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值