import rasterioimport numpy as npdefparse_mtl(mtl_path):"""解析Landsat8 MTL元数据，提取定标系数、太阳高度角"""    meta = {}withopen(mtl_path, 'r') as f:for line in f:            line = line.strip()if'='in line andnot line.startswith('#'):                k, v = line.split('=', 1)                k = k.strip()                v = v.strip().strip('"')try:                    v = float(v) if'.'in v or'e'in v elseint(v)except: pass                meta[k] = vreturn metadefdn_to_toa_reflectance(dn, band_num, meta):"""DN值转换为TOA表观反射率"""    mult = meta[f'REFLECTANCE_MULT_BAND_{band_num}']    add = meta[f'REFLECTANCE_ADD_BAND_{band_num}']    sun_elev_rad = np.deg2rad(meta['SUN_ELEVATION'])    toa = (dn.astype(np.float64) * mult + add) / np.sin(sun_elev_rad)return np.clip(toa, -0.2, 1.2)defdos_atmospheric_correction(toa):"""DOS暗目标大气校正"""    valid = toa[toa > 0.01]    dark_val = np.percentile(valid, 1.0)    surf_ref = toa - dark_valreturn np.clip(surf_ref, 0, None)defbuild_cloud_mask(bqa_path):"""BQA波段提取云、云阴影、冰雪掩膜"""with rasterio.open(bqa_path) as src:        bqa = src.read(1)    fill = ((bqa >> 0) & 1).astype(bool)    cloud = ((bqa >> 3) & 1).astype(bool)    shadow = ((bqa >> 4) & 1).astype(bool)    ice = ((bqa >> 6) & 1).astype(bool)    invalid = fill | cloud | shadow | icereturn invalid

import torchimport torch.nn as nnimport torch.nn.functional as FclassChannelAttention(nn.Module):def__init__(self, in_ch, reduction=16):super().__init__()self.avg_pool = nn.AdaptiveAvgPool2d(1)self.max_pool = nn.AdaptiveMaxPool2d(1)self.fc = nn.Sequential(            nn.Linear(in_ch, in_ch//16),            nn.ReLU(),            nn.Linear(in_ch//16, in_ch)        )self.sigmoid = nn.Sigmoid()defforward(self, x):        b,c,_,_ = x.shape        avg = self.fc(self.avg_pool(x).view(b,c))max = self.fc(self.max_pool(x).view(b,c))        att = self.sigmoid((avg+max).view(b,c,1,1))return x * attclassSpatialAttention(nn.Module):def__init__(self, k=7):super().__init__()self.conv = nn.Conv2d(2,1,k,padding=k//2,bias=False)self.sigmoid = nn.Sigmoid()defforward(self, x):        avg = torch.mean(x, dim=1, keepdim=True)max,_ = torch.max(x, dim=1, keepdim=True)        att = self.sigmoid(self.conv(torch.cat([avg,max],dim=1)))return x * attclassCBAM(nn.Module):def__init__(self, c):super().__init__()self.ca = ChannelAttention(c)self.sa = SpatialAttention()defforward(self,x):returnself.sa(self.ca(x))

class ASPP(nn.Module):    def __init__(self, in_ch, out_ch, dilations=[1,6,12,18]):super().__init__()        self.branches = nn.ModuleList()        for d in dilations:            if d == 1:                self.branches.append(nn.Conv2d(in_ch, out_ch, 1, bias=False))            else:                self.branches.append(nn.Conv2d(in_ch, out_ch,3,padding=d,dilation=d,bias=False))        self.global_pool = nn.Sequential(nn.AdaptiveAvgPool2d(1),nn.Conv2d(in_ch,out_ch,1))        self.fuse = nn.Conv2d(out_ch*5, out_ch, 1)    def forward(self, x):        h,w = x.shape[2:]        feats = []        for b in self.branches:            feats.append(b(x))        gp = self.global_pool(x)        gp = F.interpolate(gp, (h,w), mode='bilinear',align_corners=False)        feats.append(gp)        return self.fuse(torch.cat(feats,dim=1))

classCombinedLoss(nn.Module):def__init__(self, bce_w=0.5, dice_w=0.5, smooth=1e-6):super().__init__()self.bce = nn.BCELoss()self.bce_w = bce_wself.dice_w = dice_wself.smooth = smoothdefdice_loss(self, pred, target):        pred = pred.view(-1)        target = target.view(-1)        inter = (pred*target).sum()        dice = (2*inter + self.smooth)/(pred.sum()+target.sum()+self.smooth)return1 - dicedefforward(self,pred,target):        loss_bce = self.bce(pred, target)        loss_dice = self.dice_loss(pred, target)returnself.bce_w * loss_bce + self.dice_w * loss_dice

classMetricsCalculator:def__init__(self):self.reset()defreset(self):self.tp=self.fp=self.tn=self.fn=0defupdate(self, pred, target, th=0.5):        pred_bin = (pred>th).float()        tar_bin = (target>0.5).float()self.tp += (pred_bin*tar_bin).sum().item()self.fp += (pred_bin*(1-tar_bin)).sum().item()self.tn += ((1-pred_bin)*(1-tar_bin)).sum().item()self.fn += ((1-pred_bin)*tar_bin).sum().item()defget_metrics(self, smooth=1e-6):        iou = (self.tp+smooth)/(self.tp+self.fp+self.fn+smooth)        prec = (self.tp+smooth)/(self.tp+self.fp+smooth)        recall = (self.tp+smooth)/(self.tp+self.fn+smooth)        f1 = 2*prec*recall/(prec+recall+smooth)        oa = (self.tp+self.tn)/(self.tp+self.fp+self.tn+self.fn)return {"IoU":iou,"F1":f1,"Precision":prec,"Recall":recall,"OA":oa}

deffull_image_predict(model, ms_path, out_tif, tile_size=256, device="cuda"):"""整景Landsat大图分块滑窗预测，避免显存溢出"""import rasteriowith rasterio.open(ms_path) as src:        data = src.read()        h,w = src.height, src.width        profile = src.profile标准化参数（训练集统计值）    mean = np.array([0.0481,0.0716,0.1185,0.1826,0.2406,0.2090]).reshape(6,1,1)    std = np.array([0.0172,0.0250,0.0410,0.0460,0.0639,0.0677]).reshape(6,1,1)    pred_full = np.zeros((h,w), dtype=np.float32)滑窗遍历for y inrange(0, h, tile_size):for x inrange(0, w, tile_size):            tw = min(tile_size, w-x)            th = min(tile_size, h-y)            tile = data[:, y:y+th, x:x+tw].astype(np.float32)            tile = (tile - mean) / std            tensor = torch.from_numpy(tile).unsqueeze(0).to(device)with torch.no_grad():                res = model(tensor).squeeze().cpu().numpy()            pred_full[y:y+th, x:x+tw] = res# 输出概率栅格    profile.update(count=1, dtype="float32")with rasterio.open(out_tif, "w", **profile) as dst:        dst.write(pred_full, 1)# 输出二值建筑图    binary = (pred_full>0.5).astype(np.uint8)    bin_path = out_tif.replace(".tif","_binary.tif")    profile.update(dtype="uint8")with rasterio.open(bin_path, "w", **profile) as dst:        dst.write(binary,1)return pred_full, binary

import geopandas as gpdimport rasteriofrom rasterio.mask import maskdefclip_raster_by_shp(raster_path, shp_path, out_path):"""矢量行政边界裁剪栅格，自动投影匹配"""    gdf = gpd.read_file(shp_path)with rasterio.open(raster_path) as src:if gdf.crs != src.crs:            gdf = gdf.to_crs(src.crs)        geom = [i.__geo_interface__ for i in gdf.geometry]        img, trans = mask(src, geom, crop=True, nodata=0)        meta = src.meta.copy()        meta.update(height=img.shape[1], width=img.shape[2], transform=trans)with rasterio.open(out_path, "w",**meta) as dst:        dst.write(img)统计建筑面积    valid = img[0] != 0    building = (img[0]>0.5) & valid    pixel_area = 30*30 / 1e6    total_km2 = valid.sum() * pixel_area    build_km2 = building.sum() * pixel_area    coverage = build_km2 / total_km2print(f"总面积:{total_km2:.2f}km²,建筑面积:{build_km2:.2f}km²,覆盖率:{coverage:.2%}")return img, total_km2, build_km2, coverage