为什么输出结果全是空白？
修正他
import os
import glob
import rasterio
import numpy as np
import pandas as pd
import netCDF4 as nc
from rasterio.warp import calculate_default_transform, reproject, Resampling
from tqdm import tqdm

# 定义采样网格
lon_min, lon_max = 73, 135   # 中国大致经度范围
lat_min, lat_max = 18, 54    # 中国大致纬度范围
grid_resolution = 0.1  # 采样间隔为0.1度

lon_grid = np.arange(lon_min, lon_max + grid_resolution, grid_resolution)
lat_grid = np.arange(lat_min, lat_max + grid_resolution, grid_resolution)
lon_mesh, lat_mesh = np.meshgrid(lon_grid, lat_grid)

lon_flat = np.round(lon_mesh.flatten(), 4).astype(np.float32)
lat_flat = np.round(lat_mesh.flatten(), 4).astype(np.float32)
grid_df = pd.DataFrame({'Longitude': lon_flat, 'Latitude': lat_flat})

print(f"总的采样点数：{len(grid_df)}")

# 数据文件路径（请根据实际路径修改）
dem_file = r"F:\HuaweiCup\DATA\data1\Albers_105\TIFF\chdem_105.tif"
slope_file = r"F:\HuaweiCup\DATA\data1\Albers_105\TIFF\chdem_Slope.tif"
aspect_file = r"F:\HuaweiCup\DATA\data1\Albers_105\TIFF\chdem_Aspect.tif"
temperature_dir = r"F:\HuaweiCup\DATA\data2\day_avg"
precipitation_file = r"F:\HuaweiCup\DATA\data3\CHM_PRE_0.25dg_19612022.nc\CHM_PRE_0.25dg_19612022.nc"
land_cover_dir = r"F:\HuaweiCup\DATA\data4\data"
pop_dir = r"F:\HuaweiCup\DATA\data5\pop_output"
gdp_dir = r"F:\HuaweiCup\DATA\data6\GDP_output"
output_csv = r"F:\HuaweiCup\DATA\combined_data.csv"

# 重投影 DEM 数据
def reproject_raster(src_file):
    with rasterio.open(src_file) as src:
        transform, width, height = calculate_default_transform(
            src.crs, 'EPSG:4326', src.width, src.height, *src.bounds)
        data = src.read(
            1,
            out_shape=(height, width),
            resampling=Resampling.nearest
        )
        nodata = src.nodata
        src_transform = src.transform
        src_crs = src.crs

    destination = np.empty((height, width), dtype=data.dtype)
    reproject(
        source=data,
        destination=destination,
        src_transform=src_transform,
        src_crs=src_crs,
        dst_transform=transform,
        dst_crs='EPSG:4326',
        resampling=Resampling.nearest)

    return destination, transform, nodata

dem_data_wgs84, dem_transform_wgs84, dem_nodata = reproject_raster(dem_file)
slope_data_wgs84, _, slope_nodata = reproject_raster(slope_file)
aspect_data_wgs84, _, aspect_nodata = reproject_raster(aspect_file)

dem_width = dem_data_wgs84.shape[1]
dem_height = dem_data_wgs84.shape[0]
dem_lon = np.linspace(dem_transform_wgs84[2], dem_transform_wgs84[2] + dem_transform_wgs84[0] * dem_width, dem_width, endpoint=False).astype(np.float32)
dem_lat = np.linspace(dem_transform_wgs84[5], dem_transform_wgs84[5] + dem_transform_wgs84[4] * dem_height, dem_height, endpoint=False).astype(np.float32)
dem_lat = dem_lat[::-1]

dem_lon_grid, dem_lat_grid = np.meshgrid(dem_lon, dem_lat)
dem_df = pd.DataFrame({
    'Longitude': dem_lon_grid.flatten(),
    'Latitude': dem_lat_grid.flatten(),
    'Elevation': dem_data_wgs84.flatten(),
    'Slope': slope_data_wgs84.flatten(),
    'Aspect': aspect_data_wgs84.flatten()
})

dem_df = dem_df[(dem_df['Elevation'] != dem_nodata) & (dem_df['Slope'] != slope_nodata) & (dem_df['Aspect'] != aspect_nodata)]
dem_df['Longitude'] = np.round(dem_df['Longitude'], 4)
dem_df['Latitude'] = np.round(dem_df['Latitude'], 4)

dem_df['Longitude'] = dem_df['Longitude'].astype(np.float32)
dem_df['Latitude'] = dem_df['Latitude'].astype(np.float32)
dem_df['Elevation'] = dem_df['Elevation'].astype(np.float32)
dem_df['Slope'] = dem_df['Slope'].astype(np.float32)
dem_df['Aspect'] = dem_df['Aspect'].astype(np.float32)

grid_df = pd.merge(grid_df, dem_df, on=['Longitude', 'Latitude'], how='left')

print("完成高程数据的合并。")

# 准备写入 CSV 文件
if os.path.exists(output_csv):
    os.remove(output_csv)

headers = ['Longitude', 'Latitude', 'Elevation', 'Slope', 'Aspect', 'Year',
           'Max_Temperature', 'Max_Precipitation', 'Land_Cover', 'Population', 'GDP']
with open(output_csv, 'w', encoding='utf-8') as f:
    f.write(','.join(headers) + '\n')

# 逐年处理数据
start_year = 1979  # 根据数据实际年份调整
end_year = 2018    # 根据数据实际年份调整

# 打开降水数据集
precip_dataset = nc.Dataset(precipitation_file)
precipitation = precip_dataset.variables['pre'][:]  # time x lat x lon
time_var = precip_dataset.variables['time']
time_units = time_var.units
dates = nc.num2date(time_var[:], units=time_units)
years_precip = np.array([date.year for date in dates])
unique_years_precip = np.unique(years_precip)
year_indices_precip = {year: np.where(years_precip == year)[0] for year in unique_years_precip}

lats_precip = precip_dataset.variables['latitude'][:]
lons_precip = precip_dataset.variables['longitude'][:]
lon_mesh_precip, lat_mesh_precip = np.meshgrid(lons_precip, lats_precip)
lons_precip_rounded = np.round(lon_mesh_precip.flatten(), 4).astype(np.float32)
lats_precip_rounded = np.round(lat_mesh_precip.flatten(), 4).astype(np.float32)

for year in tqdm(range(start_year, end_year + 1), desc='逐年处理数据'):
    yearly_df = grid_df.copy()
    yearly_df['Year'] = year

    # --- 处理气温数据 ---
    year_dir = os.path.join(temperature_dir, f"{year}_avg")
    if os.path.exists(year_dir):
        tiff_files = sorted(glob.glob(os.path.join(year_dir, '*_avg.tif')))
        max_temp_year = None
        for tif_file in tiff_files:
            with rasterio.open(tif_file) as src:
                data = src.read(1)
                nodata = src.nodata
                transform = src.transform
                if nodata is not None:
                    data = np.where(data == nodata, np.nan, data)
                else:
                    data = data.astype(float)
                if max_temp_year is None:
                    max_temp_year = data
                else:
                    max_temp_year = np.fmax(max_temp_year, data)
        if max_temp_year is not None:
            rows, cols = np.indices(max_temp_year.shape)
            xs, ys = rasterio.transform.xy(transform, rows, cols)
            lons = np.array(xs).flatten()
            lats = np.array(ys).flatten()
            lons = np.round(lons, 4).astype(np.float32)
            lats = np.round(lats, 4).astype(np.float32)
            temp_df = pd.DataFrame({
                'Longitude': lons,
                'Latitude': lats,
                'Max_Temperature': max_temp_year.flatten().astype(np.float32)
            })
            temp_df = temp_df[~np.isnan(temp_df['Max_Temperature'])]
            yearly_df = pd.merge(yearly_df, temp_df, on=['Longitude', 'Latitude'], how='left')
        else:
            yearly_df['Max_Temperature'] = np.nan
    else:
        yearly_df['Max_Temperature'] = np.nan

    # --- 处理降水数据 ---
    if year in year_indices_precip:
        indices = year_indices_precip[year]
        precip_year = precipitation[indices, :, :]
        missing_value = -99.9
        precip_year = np.where(precip_year == missing_value, np.nan, precip_year)
        max_precip_year = np.nanmax(precip_year, axis=0)
        precip_df = pd.DataFrame({
            'Longitude': lons_precip_rounded,
            'Latitude': lats_precip_rounded,
            'Max_Precipitation': max_precip_year.flatten().astype(np.float32)
        })
        precip_df = precip_df[~np.isnan(precip_df['Max_Precipitation'])]
        yearly_df = pd.merge(yearly_df, precip_df, on=['Longitude', 'Latitude'], how='left')
    else:
        yearly_df['Max_Precipitation'] = np.nan

    # --- 处理土地利用数据 ---
    land_cover_arrays = {}
    cover_types = {
        'cropland': 'cropland',
        'forest': 'forest',
        'grass': 'grass',
        'shrub': 'shrub',
        'wetland': 'wetland'
    }
    for cover_type_en in cover_types.keys():
        tif_file = os.path.join(land_cover_dir, f"{cover_type_en}-{year}.tif")
        if not os.path.exists(tif_file):
            continue
        with rasterio.open(tif_file) as src:
            data = src.read(1)
            nodata = src.nodata
            transform = src.transform
            if nodata is not None:
                data = np.where(data == nodata, np.nan, data)
            else:
                data = data.astype(float)
            land_cover_arrays[cover_type_en] = data
    if land_cover_arrays:
        stacked_data = np.stack(list(land_cover_arrays.values()))
        max_indices = np.nanargmax(stacked_data, axis=0)
        cover_types_list = list(land_cover_arrays.keys())
        dominant_cover = np.array(cover_types_list)[max_indices]
        rows, cols = np.indices(dominant_cover.shape)
        xs, ys = rasterio.transform.xy(transform, rows, cols)
        lons = np.array(xs).flatten()
        lats = np.array(ys).flatten()
        lons = np.round(lons, 4).astype(np.float32)
        lats = np.round(lats, 4).astype(np.float32)
        land_cover_df = pd.DataFrame({
            'Longitude': lons,
            'Latitude': lats,
            'Land_Cover': dominant_cover.flatten()
        })
        yearly_df = pd.merge(yearly_df, land_cover_df, on=['Longitude', 'Latitude'], how='left')
    else:
        yearly_df['Land_Cover'] = np.nan

    # --- 处理人口数据 ---
    if 1990 <= year <= 2015:
        tif_file = os.path.join(pop_dir, f"pop{year}", f"pop{year}.tif")
        if os.path.exists(tif_file):
            with rasterio.open(tif_file) as src:
                data = src.read(1)
                nodata = src.nodata
                transform = src.transform
                if nodata is not None:
                    data = np.where(data == nodata, np.nan, data)
                else:
                    data = data.astype(float)
                rows, cols = np.indices(data.shape)
                xs, ys = rasterio.transform.xy(transform, rows, cols)
                lons = np.array(xs).flatten()
                lats = np.array(ys).flatten()
                lons = np.round(lons, 4).astype(np.float32)
                lats = np.round(lats, 4).astype(np.float32)
                pop_df = pd.DataFrame({
                    'Longitude': lons,
                    'Latitude': lats,
                    'Population': data.flatten().astype(np.float32)
                })
                pop_df = pop_df[~np.isnan(pop_df['Population'])]
                yearly_df = pd.merge(yearly_df, pop_df, on=['Longitude', 'Latitude'], how='left')
        else:
            yearly_df['Population'] = np.nan
    else:
        yearly_df['Population'] = np.nan

    # --- 处理GDP数据 ---
    if 1990 <= year <= 2015:
        tif_file = os.path.join(gdp_dir, f"gdp{year}.tif")
        if os.path.exists(tif_file):
            with rasterio.open(tif_file) as src:
                data = src.read(1)
                nodata = src.nodata
                transform = src.transform
                if nodata is not None:
                    data = np.where(data == nodata, np.nan, data)
                else:
                    data = data.astype(float)
                rows, cols = np.indices(data.shape)
                xs, ys = rasterio.transform.xy(transform, rows, cols)
                lons = np.array(xs).flatten()
                lats = np.array(ys).flatten()
                lons = np.round(lons, 4).astype(np.float32)
                lats = np.round(lats, 4).astype(np.float32)
                gdp_df = pd.DataFrame({
                    'Longitude': lons,
                    'Latitude': lats,
                    'GDP': data.flatten().astype(np.float32)
                })
                gdp_df = gdp_df[~np.isnan(gdp_df['GDP'])]
                yearly_df = pd.merge(yearly_df, gdp_df, on=['Longitude', 'Latitude'], how='left')
        else:
            yearly_df['GDP'] = np.nan
    else:
        yearly_df['GDP'] = np.nan

    # 重新排列列的顺序
    yearly_df = yearly_df[headers]

    # 追加到 CSV 文件
    yearly_df.to_csv(output_csv, mode='a', header=False, index=False, encoding='utf-8')

# 关闭降水数据集
precip_dataset.close()

print(f"所有数据处理完毕，合并后的数据已保存到 {output_csv}")