卫星图像另类数据：用卫星图像预测公司业绩与股价

# 示例：使用Sentinel API下载免费卫星图像
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt

# 连接Copernicus Open Access Hub
api = SentinelAPI('username', 'password', 'https://apihub.copernicus.eu/apihub')

# 定义感兴趣区域（ROI）：沃尔玛总部周围10km
roi = geojson_to_wkt(read_geojson('walmart_headquarters.geojson'))

# 查询2026年5月的Sentinel-2图像
products = api.query(
    roi,
    date=('20260501', '20260531'),
    platformname='Sentinel-2',
    producttype='S2MSI1C',
    cloudcoverpercentage=(0, 10)  # 云量<10%
)

# 下载图像
api.download_all(products)
print(f"下载了 {len(products)} 景 Sentinel-2 图像")

import rasterio
import numpy as np
from rasterio.plot import show

def load_satellite_image(image_path):
    """
    加载卫星图像并提取RGB波段
    """
    with rasterio.open(image_path) as src:
        # 读取RGB波段（通常是第3、2、1波段）
        red = src.read(3)
        green = src.read(2)
        blue = src.read(1)
        
        # 堆叠为RGB图像
        rgb = np.dstack((red, green, blue))
        
        # 归一化到0-255
        rgb_normalized = (rgb / rgb.max() * 255).astype(np.uint8)
        
        # 获取地理元数据
        meta = src.meta
        bounds = src.bounds
        
    return rgb_normalized, meta, bounds

# 示例：加载沃尔玛停车场卫星图像
image, meta, bounds = load_satellite_image('walmart_parking_lot_20260515.tif')
show(image)
print(f"图像分辨率: {meta['width']}x{meta['height']}")
print(f"地理范围: {bounds}")

import torch
from torchvision.models.detection import fasterrcnn_resnet50_fpn
from torchvision.transforms import functional as F

def detect_parking_cars(image_path, confidence_threshold=0.7):
    """
    使用Faster R-CNN检测停车场车辆
    """
    # 加载预训练模型
    model = fasterrcnn_resnet50_fpn(pretrained=True)
    model.eval()
    
    # 读取图像
    image, _, _ = load_satellite_image(image_path)
    image_tensor = F.to_tensor(image)
    
    # 推理
    with torch.no_grad():
        predictions = model([image_tensor])
    
    # 解析预测结果
    cars = []
    for box, label, score in zip(predictions[0]['boxes'], 
                                  predictions[0]['labels'], 
                                  predictions[0]['scores']):
        if label == 3 and score > confidence_threshold:  # COCO数据集中label=3是汽车
            cars.append({
                'bbox': box.numpy(),
                'confidence': score.numpy()
            })
    
    return cars

# 示例：分析沃尔玛停车场
image_path = 'walmart_parking_lot_20260515.tif'
detected_cars = detect_parking_cars(image_path)

# 计算停车场占用率
total_parking_spaces = 500  # 假设停车场有500个车位
occupancy_rate = len(detected_cars) / total_parking_spaces
print(f"检测到 {len(detected_cars)} 辆车")
print(f"停车场占用率: {occupancy_rate:.2%}")

# 可视化检测结果
import cv2
image, _, _ = load_satellite_image(image_path)
for car in detected_cars:
    box = car['bbox'].astype(int)
    cv2.rectangle(image, (box[0], box[1]), (box[2], box[3]), (0, 255, 0), 2)

cv2.imwrite('detection_result.jpg', image)

def build_parking_flow_index(image_paths, store_metadata):
    """
    构建停车场车流指数
    """
    flow_data = []
    
    for img_path, date in image_paths:
        # 检测车辆
        cars = detect_parking_cars(img_path)
        
        # 计算占用率
        store_id = extract_store_id(img_path)
        total_spaces = store_metadata[store_id]['parking_spaces']
        occupancy = len(cars) / total_spaces
        
        flow_data.append({
            'date': date,
            'store_id': store_id,
            'car_count': len(cars),
            'occupancy_rate': occupancy
        })
    
    # 构建面板数据
    df = pd.DataFrame(flow_data)
    df['week'] = pd.to_datetime(df['date']).dt.isocalendar().week
    
    # 计算同店同比增长（YoY）
    df = df.sort_values(['store_id', 'date'])
    df['occupancy_yoy'] = df.groupby('store_id')['occupancy_rate'].pct_change(periods=52)
    
    return df

# 示例：构建沃尔玛车流指数
walmart_images = [
    ('walmart_2024W20.tif', '2024-05-15'),
    ('walmart_2025W20.tif', '2025-05-14'),
    ('walmart_2026W20.tif', '2026-05-13')
]
parking_index = build_parking_flow_index(walmart_images, store_metadata)

print("沃尔玛停车场车流指数（2024-2026）：")
print(parking_index[['date', 'occupancy_rate', 'occupancy_yoy']])

import numpy as np
from scipy import ndimage

def estimate_tank_level(image, tank_coordinates):
    """
    通过阴影分析估算储油罐液位
    """
    # 提取储油罐ROI
    x, y, r = tank_coordinates  # 圆心(x,y)，半径r
    tank_roi = image[y-r:y+r, x-r:x+r]
    
    # 边缘检测（Canny）
    edges = cv2.Canny(tank_roi, 100, 200)
    
    # 检测浮顶阴影（假设阴影在罐体左侧）
    shadow_region = tank_roi[:, :r//2]
    shadow_length = np.sum(shadow_region < 50)  # 阴影区域像素值较低
    
    # 根据阴影长度推算液位（需要标定）
    # 假设线性关系：液位 = a * 阴影长度 + b
    a, b = 0.8, 10  # 标定系数
    liquid_level = a * shadow_length + b
    
    # 计算库存（假设圆柱形储罐）
    tank_height = 20  # 米
    tank_radius = 15   # 米
    volume = np.pi * tank_radius**2 * (liquid_level / tank_height * tank_height)
    
    return liquid_level, volume

# 示例：监测库欣（Cushing）储油基地
cushing_image = load_sar_image('cushing_oklahoma_20260610.tif')  # SAR图像

# 库欣基地有约100个大型储油罐
tank_coordinates = [
    (1250, 800, 50),   # 罐1
    (1350, 820, 52),   # 罐2
    # ... 其他储罐坐标
]

total_inventory = 0
for coords in tank_coordinates:
    level, volume = estimate_tank_level(cushing_image, coords)
    total_inventory += volume
    print(f"储罐液位: {level:.1f}m, 库存: {volume:.0f}桶")

print(f"库欣总库存: {total_inventory/1e6:.1f}百万桶")

# 与EIA数据对比
eia_report = pd.read_csv('eia_cushing_inventory.csv')
latest_eia = eia_report.iloc[-1]['inventory_mbbl']
print(f"EIA最新报告: {latest_eia:.1f}百万桶")
print(f"卫星估算偏差: {(total_inventory/1e6 - latest_eia)/latest_eia:.2%}")

def calculate_ndvi(nir_band, red_band):
    """
    计算归一化植被指数（NDVI）
    NDVI = (NIR - Red) / (NIR + Red)
    """
    ndvi = (nir_band.astype(float) - red_band.astype(float)) / \
           (nir_band + red_band)
    
    # NDVI范围：-1 to 1，植被通常0.2-0.8
    ndvi = np.clip(ndvi, -1, 1)
    
    return ndvi

def monitor_crop_health(image_path, field_boundary):
    """
    监测农田作物健康状况
    """
    # 加载多光谱图像
    with rasterio.open(image_path) as src:
        # Sentinel-2波段：B4(Red), B8(NIR)
        red = src.read(4).astype(float)
        nir = src.read(8).astype(float)
        
        # 创建掩膜（仅保留农田区域）
        mask = rasterio.features.geometry_mask(
            [field_boundary],
            out_shape=red.shape,
            transform=src.transform,
            invert=True
        )
        
        # 计算NDVI
        ndvi = calculate_ndvi(nir, red)
        ndvi_masked = np.ma.array(ndvi, mask=~mask)
        
    # 统计指标
    mean_ndvi = ndvi_masked.mean()
    health_status = classify_crop_health(mean_ndvi)
    
    return ndvi, mean_ndvi, health_status

def classify_crop_health(ndvi_value):
    """
    根据NDVI判断作物健康状态
    """
    if ndvi_value < 0.2:
        return "裸土/植被稀疏"
    elif ndvi_value < 0.4:
        return "植被生长初期"
    elif ndvi_value < 0.7:
        return "健康植被"
    else:
        return "非常茂密/可能病害"

# 示例：监测爱荷华州玉米田
iowa_corn_field = read_geojson('iowa_corn_field.geojson')
image_path = 'sentinel2_iowa_20260605.tif'

ndvi, mean_ndvi, status = monitor_crop_health(image_path, iowa_corn_field)
print(f"平均NDVI: {mean_ndvi:.3f}")
print(f"作物健康状态: {status}")

# 可视化NDVI分布
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 8))
plt.imshow(ndvi, cmap='RdYlGn', vmin=0, vmax=1)
plt.colorbar(label='NDVI')
plt.title('爱荷华州玉米田NDVI分布（2026-06-05）')
plt.savefig('corn_ndvi_map.png', dpi=300, bbox_inches='tight')

from sklearn.ensemble import RandomForestRegressor

def predict_crop_yield(ndvi_time_series, weather_data, historical_yield):
    """
    基于NDVI时间序列和气象数据预测作物产量
    """
    # 构建特征工程
    features = []
    
    for field_id in ndvi_time_series['field_id'].unique():
        field_ndvi = ndvi_time_series[ndvi_time_series['field_id'] == field_id]
        
        # NDVI特征
        ndvi_max = field_ndvi['ndvi'].max()  # 生长季峰值
        ndvi_mean = field_ndvi['ndvi'].mean()  # 平均生物量
        ndvi_duration = (field_ndvi['ndvi'] > 0.6).sum()  # 健康植被持续时间
        
        # 气象特征
        field_weather = weather_data[weather_data['field_id'] == field_id]
        temp_sum = field_weather['temperature'].sum()  # 积温
        precipitation = field_weather['precipitation'].sum()  # 总降水
        
        features.append([
            ndvi_max, ndvi_mean, ndvi_duration,
            temp_sum, precipitation
        ])
    
    # 训练随机森林模型
    X = np.array(features)
    y = historical_yield['yield_bushels_per_acre'].values
    
    model = RandomForestRegressor(n_estimators=100, random_state=42)
    model.fit(X, y)
    
    # 预测2026年产量
    predicted_yield = model.predict(X)
    
    return model, predicted_yield

# 示例：预测2026年美国玉米产量
us_corn_ndvi = pd.read_csv('us_corn_ndvi_2026.csv')
us_corn_weather = pd.read_csv('us_corn_weather_2026.csv')
us_corn_historical = pd.read_csv('us_corn_yield_2016_2025.csv')

model, predicted_2026 = predict_crop_yield(
    us_corn_ndvi, us_corn_weather, us_corn_historical
)

print("2026年美国玉米产量预测：")
print(f"预测平均值: {predicted_2026.mean():.0f} 蒲式耳/英亩")
print(f"USDA 2025实际值: {us_corn_historical['yield_bushels_per_acre'].iloc[-1]:.0f}")

def satellite_event_driven_strategy(satellite_signals, financial_reports, threshold=2.0):
    """
    卫星图像事件驱动策略
    """
    trades = []
    
    for signal in satellite_signals:
        # 计算信号强度（Z-Score）
        z_score = (signal['value'] - signal['historical_mean']) / signal['historical_std']
        
        if abs(z_score) > threshold:
            # 生成交易信号
            if z_score > 0:
                direction = 'long'  # 车流/库存超预期
            else:
                direction = 'short'  # 车流/库存低于预期
            
            # 等待财报验证
            earnings_date = financial_reports[
                (financial_reports['ticker'] == signal['ticker']) &
                (financial_reports['date'] > signal['date'])
            ]['date'].min()
            
            trades.append({
                'ticker': signal['ticker'],
                'entry_date': signal['date'],
                'direction': direction,
                'z_score': z_score,
                'earnings_date': earnings_date,
                'holding_days': (earnings_date - signal['date']).days
            })
    
    return pd.DataFrame(trades)

# 示例：基于停车场车流的交易信号
parking_signals = [
    {'ticker': 'WMT', 'date': '2026-05-15', 'value': 0.85, 'historical_mean': 0.65, 'historical_std': 0.08},
    {'ticker': 'TGT', 'date': '2026-05-15', 'value': 0.55, 'historical_mean': 0.68, 'historical_std': 0.07}
]

trades = satellite_event_driven_strategy(parking_signals, earnings_calendar, threshold=2.0)
print("卫星图像事件驱动策略交易信号：")
print(trades[['ticker', 'direction', 'z_score', 'holding_days']])

def build_satellite_factor(satellite_data, returns, lookback=252):
    """
    构建卫星图像另类数据因子
    """
    # 标准化卫星数据
    satellite_data['satellite_zscore'] = satellite_data.groupby('date')['value'].transform(
        lambda x: (x - x.mean()) / x.std()
    )
    
    # 计算因子收益率
    factor_returns = []
    
    for date in satellite_data['date'].unique():
        if date < satellite_data['date'].min() + pd.Timedelta(days=lookback):
            continue
        
        # 根据卫星数据排序股票
        daily_data = satellite_data[satellite_data['date'] == date]
        daily_data = daily_data.sort_values('satellite_zscore', ascending=False)
        
        # 多空组合：做多卫星信号最强的10%，做空最弱的10%
        long_portfolio = daily_data.head(int(len(daily_data) * 0.1))['ticker'].tolist()
        short_portfolio = daily_data.tail(int(len(daily_data) * 0.1))['ticker'].tolist()
        
        # 计算组合收益
        long_return = returns[long_portfolio].loc[date:date+pd.Timedelta(days=20)].mean().mean()
        short_return = returns[short_portfolio].loc[date:date+pd.Timedelta(days=20)].mean().mean()
        
        factor_return = long_return - short_return
        factor_returns.append({
            'date': date,
            'factor_return': factor_return,
            'long_return': long_return,
            'short_return': short_return
        })
    
    return pd.DataFrame(factor_returns)

# 示例：构建停车场车流因子
parking_data = pd.read_csv('parking_flow_weekly.csv')
stock_returns = pd.read_csv('stock_returns_weekly.csv', index_col='date', parse_dates=True)

satellite_factor = build_satellite_factor(parking_data, stock_returns, lookback=52)
print("卫星图像因子表现：")
print(f"因子年化收益率: {satellite_factor['factor_return'].mean() * 52:.2%}")
print(f"因子夏普比率: {satellite_factor['factor_return'].mean() / satellite_factor['factor_return'].std() * np.sqrt(52):.2f}")

import statsmodels.api as sm

def predict_walmart_revenue(parking_data, financial_data, macro_data):
    """
    预测沃尔玛季度营收
    """
    # 合并数据
    merged = []
    
    for quarter in financial_data['quarter'].unique():
        # 卫星特征：季度平均车流、同比增长
        quarter_parking = parking_data[
            (parking_data['date'] >= quarter) &
            (parking_data['date'] < pd.to_datetime(quarter) + pd.offsets.QuarterEnd(1))
        ]
        parking_feature = quarter_parking['occupancy_rate'].mean()
        parking_yoy = quarter_parking['occupancy_yoy'].mean()
        
        # 宏观特征
        quarter_macro = macro_data[macro_data['quarter'] == quarter].iloc[0]
        consumer_confidence = quarter_macro['consumer_confidence']
        gas_price = quarter_macro['gas_price']
        
        # 目标变量
        revenue = financial_data[financial_data['quarter'] == quarter]['revenue'].values[0]
        
        merged.append([
            quarter, parking_feature, parking_yoy,
            consumer_confidence, gas_price, revenue
        ])
    
    df = pd.DataFrame(merged, columns=[
        'quarter', 'parking_flow', 'parking_yoy',
        'consumer_confidence', 'gas_price', 'revenue'
    ])
    
    # 线性回归模型
    X = df[['parking_flow', 'parking_yoy', 'consumer_confidence', 'gas_price']]
    X = sm.add_constant(X)
    y = df['revenue']
    
    model = sm.OLS(y, X).fit()
    
    print("沃尔玛营收预测模型回归结果：")
    print(model.summary())
    
    # 预测2026 Q2营收
    parking_2026q2 = parking_data[
        (parking_data['date'] >= '2026-04-01') &
        (parking_data['date'] < '2026-07-01')
    ]['occupancy_rate'].mean()
    
    prediction_features = pd.DataFrame([{
        'const': 1,
        'parking_flow': parking_2026q2,
        'parking_yoy': 0.05,  # 假设同比增长5%
        'consumer_confidence': 110,
        'gas_price': 3.2
    }])
    
    predicted_revenue = model.predict(prediction_features)[0]
    print(f"\n2026 Q2 沃尔玛营收预测: ${predicted_revenue/1e9:.2f}B")
    print(f"2025 Q2 实际营收: ${financial_data[financial_data['quarter']=='2026-04-01']['revenue'].values[0]/1e9:.2f}B")
    
    return model, predicted_revenue

# 运行预测
model, prediction = predict_walmart_revenue(
    parking_data=pd.read_csv('walmart_parking_weekly.csv'),
    financial_data=pd.read_csv('walmart_quarterly.csv'),
    macro_data=pd.read_csv('us_macro_quarterly.csv')
)

# 示例：使用Dask进行分布式图像处理
import dask
from dask import delayed
from dask.distributed import Client

def process_satellite_image_distributed(image_paths):
    """
    分布式处理卫星图像
    """
    # 启动Dask集群
    client = Client(n_workers=8, threads_per_worker=2)
    
    # 并行处理
    tasks = []
    for img_path in image_paths:
        task = delayed(detect_parking_cars)(img_path)
        tasks.append(task)
    
    results = dask.compute(*tasks)
    
    client.close()
    return results

# 处理1000景图像
image_paths = [f'walmart_parking_{i}.tif' for i in range(1000)]
results = process_satellite_image_distributed(image_paths)
print(f"处理了 {len(results)} 景图像")