import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import rasterio
import rasterstats
from rasterio.plot import show
# show()方法用來展示柵格圖形
from rasterio.plot import show_hist
# 用來展示直方圖
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

# 使用geopandas讀取矢量數(shù)據(jù)
districts = gpd.read_file('/content/drive/MyDrive/Datashpraster/Data/Districts/districts.shp')

# 使用rasterio讀取柵格數(shù)據(jù)，柵格數(shù)據(jù)和矢量數(shù)據(jù)的坐標投影需要一致
raster = rasterio.open('/content/drive/MyDrive/Datashpraster/Data/Rainfall Data Rasters/2020-4-1.tif')

# 把矢量數(shù)據(jù)和柵格數(shù)據(jù)繪制到一個axis上，這個axis不是坐標軸，而是圖形
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 20

fig, (ax1,ax2) = plt.subplots(1,2,figsize=(15,6))

show(raster, ax=ax1,title='Rainfall')
# 讀取進來的矢量數(shù)據(jù)可以直接調(diào)用gpd的plot()方法繪制
districts.plot(ax=ax1, facecolor='None', edgecolor='red')
show_hist(raster,ax=ax2,title='hist')

plt.show()

# 提取雨量柵格值到numpy數(shù)組
# 遵循GDAL規(guī)則從第一波段讀取
rainfall_data = raster.read(1)
rainfall_data

# 設(shè)置坐標變換信息
affine = raster.transform

# 準備開始進行空間分區(qū)計算
# 第一個參數(shù)是矢量分區(qū)，第二個是柵格，第三個是坐標變換信息，第四個是統(tǒng)計均值
avg_rallrain = rasterstats.zonal_stats(districts,rainfall_data,affine=affine,stats=['mean'],geojson_out=True)
# avg_rallrain

# 除了統(tǒng)計平均值之外，還有最大最小值那些

# 將Date列轉(zhuǎn)為時間型
data['Date'] = pd.to_datetime(data['Date'], infer_datetime_format=True)

# print(data)

data['Date'] = data['Date'].dt.date
print(data)

# 準備繪制圖形
fig,(ax1,ax2)= plt.subplots(2,1,figsize=(18,6))
plt.rcParams['font.size'] = 15

data.plot(x='Date', y='Average_RF_Porto', ax=ax1, kind='bar', title='Avg_Rail_Porto')
data.plot(x='Date', y='Average_RF_Faro', ax=ax2, kind='bar', title='Avg_Rail_Faro',color='red')

#自動調(diào)整圖形的分布
plt.tight_layout()
plt.show()