diff --git a/app_run.py b/app_run.py index 292e480..0cea8a1 100644 --- a/app_run.py +++ b/app_run.py @@ -1,7 +1,4 @@ import logging -import sys -import os -from datetime import datetime import tkinter as tk from tkinter import messagebox, ttk from tkcalendar import DateEntry @@ -9,9 +6,7 @@ # 導入你的處理模組 from src.api.sentinel_api import S5PFetcher -from src.processing.data_processor import S5Processor -from src.utils.logger import setup_logging -from src.config.settings import setup_directory_structure +from src.config import setup_directory_structure class SatelliteApp: diff --git a/main.py b/main.py index f3843d8..47d6e1d 100644 --- a/main.py +++ b/main.py @@ -1,27 +1,29 @@ """主程式""" import logging +from datetime import datetime from src.api.sentinel_api import S5PFetcher from src.processing.data_processor import S5Processor -from src.utils.logger import setup_logging -from src.utils.richer import rich_print -from src.utils.catalog import ProductClassInput, ProductTypeInput, ProductType -from src.config.settings import setup_directory_structure, FILTER_BOUNDARY -from src.visualization.gif_nc import animate_data + +from src.config.richer import rich_print +from src.config.catalog import ClassInput, TypeInput, PRODUCT_CONFIGS +from src.config import setup_directory_structure, FILTER_BOUNDARY logger = logging.getLogger(__name__) -def fetch_data(file_class: ProductClassInput, - file_type: ProductTypeInput, - start_date: str, - end_date: str): +def fetch_data(file_class: ClassInput, + file_type: TypeInput, + start_date: str | datetime, + end_date: str | datetime): """下載數據的工作流程""" try: + rich_print( + f"正在獲取 sentinel-5p 衛星數據 ({PRODUCT_CONFIGS[file_type].display_name}) from {start_date} to {end_date} ...") + fetcher = S5PFetcher(max_workers=3) - rich_print(f"正在獲取 sentinel-5p 衛星數據 ({ProductType[file_type].display_name}) from {start_date} to {end_date} ...") products = fetcher.fetch_data( file_class=file_class, file_type=file_type, @@ -33,7 +35,8 @@ def fetch_data(file_class: ProductClassInput, if products: if rich_print("是否要下載數據？", confirm=True): - rich_print(f"開始下載 sentinel-5p 衛星數據 ({ProductType[file_type].display_name}) from {start_date} to {end_date} ...") + rich_print( + f"開始下載 sentinel-5p 衛星數據 ({PRODUCT_CONFIGS[file_type].display_name}) from {start_date} to {end_date} ...") fetcher.parallel_download(products) rich_print("數據下載完成！") else: @@ -47,14 +50,15 @@ def fetch_data(file_class: ProductClassInput, logger.error(error_message) -def process_data(file_class: ProductClassInput, - file_type: ProductTypeInput, - start_date: str, - end_date: str): +def process_data(file_class: ClassInput, + file_type: TypeInput, + start_date: str | datetime, + end_date: str | datetime): """處理數據的工作流程""" try: if rich_print("是否要處理數據？", confirm=True): - rich_print(f"建立 sentinel-5p 衛星數據 ({ProductType[file_type].display_name}) 處理器 ...") + rich_print( + f"正在處理 sentinel-5p 衛星數據 ({PRODUCT_CONFIGS[file_type].display_name}) from {start_date} to {end_date} ...") processor = S5Processor( interpolation_method='kdtree', @@ -62,14 +66,13 @@ def process_data(file_class: ProductClassInput, mask_qc_value=0.5 ) - rich_print(f"正在處理 sentinel-5p 衛星數據 ({ProductType[file_type].display_name}) from {start_date} to {end_date} ...") - processor.process_each_data( file_class=file_class, file_type=file_type, start_date=start_date, end_date=end_date, - use_taiwan_mask=True) + ) + rich_print("數據完成處理") else: rich_print("已取消處理操作") @@ -82,20 +85,21 @@ def process_data(file_class: ProductClassInput, def main(): # 設定參數 - start, end = '2024-01-01', '2024-11-30' + start, end = '2024-03-01', '2024-03-02' + file_class: ClassInput = 'OFFL' + file_type: TypeInput = 'NO2___' # 設定輸入輸出配置 - setup_logging() - setup_directory_structure(file_type='NO2___', start_date=start, end_date=end) + setup_directory_structure(file_type=file_type, start_date=start, end_date=end) # 下載數據 - # fetch_data(file_class='OFFL', file_type='NO2___', start_date=start, end_date=end) + # fetch_data(file_class=file_class, file_type=file_type, start_date=start, end_date=end) - # 處理數據 - process_data(file_class='OFFL', file_type='NO2___', start_date=start, end_date=end) + # 處理與繪製數據 + process_data(file_class=file_class, file_type=file_type, start_date=start, end_date=end) # 動畫 - # animate_data(file_type='NO2___', start_date=start, end_date=end) + # animate_data(file_type=file_type, start_date=start, end_date=end) if __name__ == "__main__": diff --git a/src/api/sentinel_api.py b/src/api/sentinel_api.py index 67b18a7..bb22e23 100644 --- a/src/api/sentinel_api.py +++ b/src/api/sentinel_api.py @@ -25,8 +25,8 @@ RAW_DATA_DIR, DEFAULT_TIMEOUT ) -from src.utils.richer import console, rich_print, DisplayManager -from src.utils.catalog import ProductClassInput, ProductTypeInput, ProductType +from src.config.richer import console, rich_print, DisplayManager +from src.config.catalog import ClassInput, TypeInput, PRODUCT_CONFIGS logger = logging.getLogger(__name__) @@ -64,8 +64,8 @@ def __init__(self, max_workers: int = 5): } def fetch_data(self, - file_class: ProductClassInput, - file_type: ProductTypeInput, + file_class: ClassInput, + file_type: TypeInput, start_date: str, end_date: str, boundary: tuple[float, float, float, float] | None = None, @@ -98,7 +98,7 @@ def fetch_data(self, file_class = '' if file_class == '*' else file_class file_type = '' if file_type == '*' else file_type - self.file_type = ProductType[file_type] + self.file_type = file_type # 構建基本篩選條件 base_filter = ( diff --git a/src/config/__init__.py b/src/config/__init__.py new file mode 100644 index 0000000..a4b0241 --- /dev/null +++ b/src/config/__init__.py @@ -0,0 +1,100 @@ +import logging +from pathlib import Path +from datetime import datetime + +from src.config.settings import RAW_DATA_DIR, PROCESSED_DATA_DIR, FIGURE_DIR, LOGS_DIR, FILTER_BOUNDARY +from src.config.catalog import TypeInput + + +__all__ = ['setup_directory_structure', 'FILTER_BOUNDARY'] + + +def setup_logging(): + """設置日誌配置""" + # 確保日誌目錄存在 + log_dir = Path(LOGS_DIR) + log_dir.mkdir(parents=True, exist_ok=True) + + # 創建日誌檔案路徑 + log_file = log_dir / f"Satellite_S5P_{datetime.now().strftime('%Y%m')}.log" + + # 配置基本設定 + logging.basicConfig( + level=logging.INFO, + format='%(asctime)s - %(name).10s - %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + handlers=[ + logging.FileHandler(log_file, encoding='utf-8'), + logging.StreamHandler() # 同時輸出到控制台 + ] + ) + + +def setup_directory_structure(file_type: TypeInput, + start_date: str | datetime, + end_date: str | datetime): + """確保所有必要的目錄存在""" + directories = [RAW_DATA_DIR, PROCESSED_DATA_DIR, FIGURE_DIR, LOGS_DIR] + for directory in directories: + try: + directory.mkdir(parents=True, exist_ok=True) + except Exception as e: + raise + + setup_logging() + + """依照開始和結束時間設定資料夾結構""" + start = datetime.strptime(start_date, '%Y-%m-%d') + end = datetime.strptime(end_date, '%Y-%m-%d') + + # 遍歷範圍內的所有月份 + current_date = start + while current_date <= end: + product_type = file_type + year = current_date.strftime('%Y') + month = current_date.strftime('%m') + + # 構建每個月份的 figure、processed 和 raw 路徑 + for base_dir in [FIGURE_DIR, PROCESSED_DATA_DIR, RAW_DATA_DIR]: + month_dir = base_dir / product_type / year / month + month_dir.mkdir(parents=True, exist_ok=True) + + # 移動到下個月 + if month == "12": + current_date = current_date.replace(year=current_date.year + 1, month=1, day=1) + else: + current_date = current_date.replace(month=current_date.month + 1, day=1) + + +""" I/O structure +Main Folder (Sentinel_data) +├── figure +│ ├── NO2 +│ │ ├── 2023 +│ │ │ ├── 01 +│ │ │ └── ... +│ │ └── 2024 +│ │ ├── 01 +│ │ └── ... +│ └── ... +├── logs +│ └── Satellite_S5P_202411.log +├── processed +│ ├── NO2 +│ │ ├── 2023 +│ │ │ ├── 01 +│ │ │ └── ... +│ │ └── 2024 +│ │ ├── 01 +│ │ └── ... +│ └── ... +└── raw + ├── NO2 + │ ├── 2023 + │ │ ├── 01 + │ │ └── ... + │ └── 2024 + │ ├── 01 + │ └── ... + └── ... +""" \ No newline at end of file diff --git a/src/config/catalog.py b/src/config/catalog.py new file mode 100644 index 0000000..d54ca1d --- /dev/null +++ b/src/config/catalog.py @@ -0,0 +1,149 @@ +from dataclasses import dataclass +from typing import Literal, TypeAlias +from enum import Enum + + +ClassInput = Literal['NRTI', 'OFFL', 'RPRO'] +TypeInput = Literal['O3____', 'O3_TCL', 'O3__PR', 'CH4___', 'CO____', 'NO2___', 'HCHO__', 'SO2___', 'CLOUD_', 'FRESCO', 'AER_LH', 'AER_AI'] + + +@dataclass +class ProductConfig: + """產品配置""" + display_name: str # 顯示名稱（帶下標） + dataset_name: str # 數據集名稱 + vmin: float | None # 最小值 + vmax: float | None # 最大值 + units: str # 單位 + title: str # 標題 + cmap: str = 'viridis' # 預設色階 + + +class ProductType(str, Enum): + """Available product types""" + # Ozone products + O3____ = 'O3____' + O3_TCL = 'O3_TCL' # No data ?? + O3__PR = 'O3__PR' # start from 2022 have total_vertical_columns and tropospheric_column and ozone_profile + + # Other gas products + CH4___ = 'CH4___' + CO____ = 'CO____' + NO2___ = 'NO2___' + HCHO__ = 'HCHO__' + SO2___ = 'SO2___' + + # Cloud and aerosol products + CLOUD_ = 'CLOUD_' + FRESCO = 'FRESCO' + AER_LH = 'AER_LH' + AER_AI = 'AER_AI' + + +class ProductLevel(str, Enum): + Level0 = 'L0__' + Level1B = 'L1B_' + Level2 = 'L2__' + + +class ProductClass(str, Enum): + """可用的處理類型 + + Available classes: + Main processing types: + - NRTI : Near-real time processing (最快但精確度較低) + - OFFL : Offline processing (標準處理) + - RPRO : Reprocessing (重新處理的歷史數據) + + Testing types: + - TEST : Internal testing + - OGCA : On-ground calibration + - GSOV : Ground segment validation + - OPER : Operational processing + """ + # Main processing types + NRTI = 'NRTI' + OFFL = 'OFFL' + RPRO = 'RPRO' + + # Testing types + TEST = 'TEST' + OGCA = 'OGCA' + GSOV = 'GSOV' + OPER = 'OPER' + + +PRODUCT_CONFIGS: dict[str, ProductConfig] = { + 'O3____': ProductConfig( + display_name='O\u2083', + dataset_name='ozone_total_vertical_column', + vmin=0.1, vmax=0.2, + units=f'O$_3$ Total Vertical Column (mol/m$^2$)', + title=f'O$_3$ Total Vertical Column', + cmap='RdBu_r' + ), + 'O3_TCL': ProductConfig( + display_name='O\u2083_TCL', + dataset_name='ozone_tropospheric_column', + vmin=0.1, vmax=0.2, + units=f'O$_3$ Tropospheric Column (mol/m$^2$)', + title=f'O$_3$ Tropospheric Column', + cmap='RdBu_r' + ), + 'O3__PR': ProductConfig( + display_name='O\u2083_PR', + dataset_name='ozone_tropospheric_column', + vmin=0.01, vmax=0.04, + units=f'O$_3$ Tropospheric Column (mol/m$^2$)', + title=f'O$_3$ Tropospheric Column', + cmap='jet' + ), + 'CH4___': ProductConfig( + display_name='CH\u2084', + dataset_name='methane_mixing_ratio', + vmin=None, vmax=None, + units=f'CH$_4$ Methane Mixing Ratio', + title=f'CH$_4$ Methane Mixing Ratio', + cmap='RdBu_r' + ), + 'CO____': ProductConfig( + display_name='CO', + dataset_name='carbonmonoxide_total_column', + vmin=None, vmax=None, + units=f'CO Total Column (mol/m$^2$)', + title=f'CO Total Column', + cmap='RdBu_r' + ), + 'NO2___': ProductConfig( + display_name='NO\u2082', + dataset_name='nitrogendioxide_tropospheric_column', + vmin=-2e-4, vmax=2e-4, + units=f'NO$_2$ Tropospheric Column (mol/m$^2$)', + title=f'NO$_2$ Tropospheric Column', + cmap='RdBu_r' + ), + 'HCHO__': ProductConfig( + display_name='HCHO', + dataset_name='formaldehyde_tropospheric_vertical_column', + vmin=-4e-4, vmax=4e-4, + units=f'HCHO Tropospheric Vertical Column (mol/m$^2$)', + title=f'HCHO Tropospheric Vertical Column', + cmap='RdBu_r' + ), + 'SO2___': ProductConfig( + display_name='SO\u2082', + dataset_name='sulfurdioxide_total_vertical_column', + vmin=-4e-3, vmax=4e-3, + units=f'SO$_2$ Total Vertical Column (mol/m$^2$)', + title=f'SO$_2$ Total Vertical Column', + cmap='RdBu_r' + ), + 'AER_AI': ProductConfig( + display_name='Aerosol Index', + dataset_name='aerosol_index_340_380', + vmin=None, vmax=None, + units=f'Aerosol Index', + title=f'Aerosol Index', + cmap='viridis' + ) +} diff --git a/src/utils/richer.py b/src/config/richer.py similarity index 99% rename from src/utils/richer.py rename to src/config/richer.py index 96665a4..c73258e 100644 --- a/src/utils/richer.py +++ b/src/config/richer.py @@ -7,7 +7,7 @@ # 定義常數 -PANEL_WIDTH = 90 +PANEL_WIDTH = 100 console = Console(force_terminal=True, color_system="auto", width=PANEL_WIDTH) # 使用您想要的寬度 diff --git a/src/config/settings.py b/src/config/settings.py index 96b45aa..076b8b7 100644 --- a/src/config/settings.py +++ b/src/config/settings.py @@ -1,12 +1,5 @@ """API 設定和常數""" from pathlib import Path -from datetime import datetime -import logging - -from src.utils.catalog import ProductTypeInput - - -logger = logging.getLogger(__name__) # API URLs COPERNICUS_TOKEN_URL = "https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token" @@ -33,97 +26,3 @@ FILTER_BOUNDARY = (120, 122, 22, 25.5) # (118, 124, 20, 27) FIGURE_BOUNDARY = (119, 123, 21, 26) # (100, 145, 0, 45) - - -def setup_directory_structure(file_type: ProductTypeInput, start_date: str, end_date: str): - """依照開始和結束時間設定資料夾結構""" - start = datetime.strptime(start_date, '%Y-%m-%d') - end = datetime.strptime(end_date, '%Y-%m-%d') - - # 遍歷範圍內的所有月份 - current_date = start - while current_date <= end: - product_type = file_type - year = current_date.strftime('%Y') - month = current_date.strftime('%m') - - # 構建每個月份的 figure、processed 和 raw 路徑 - for base_dir in [FIGURE_DIR, PROCESSED_DATA_DIR, RAW_DATA_DIR]: - month_dir = base_dir / product_type / year / month - month_dir.mkdir(parents=True, exist_ok=True) - - # 移動到下個月 - if month == "12": - current_date = current_date.replace(year=current_date.year + 1, month=1, day=1) - else: - current_date = current_date.replace(month=current_date.month + 1, day=1) - - -def ensure_directories(): - """確保所有必要的目錄存在""" - directories = [RAW_DATA_DIR, PROCESSED_DATA_DIR, FIGURE_DIR, LOGS_DIR] - for directory in directories: - try: - directory.mkdir(parents=True, exist_ok=True) - logger.info(f"確保目錄存在: {directory}") - except Exception as e: - logger.error(f"創建目錄失敗 {directory}: {str(e)}") - raise - - -ensure_directories() - -""" I/O structure -Main Folder (Sentinel_data) -├── figure -│ ├── NO2 -│ │ ├── 2023 -│ │ │ ├── 01 -│ │ │ └── ... -│ │ └── 2024 -│ │ ├── 01 -│ │ └── ... -│ ├── SO2 -│ │ ├── 2023 -│ │ │ ├── 01 -│ │ │ └── ... -│ │ └── 2024 -│ │ ├── 01 -│ │ └── ... -│ └── ... -├── logs -│ └── Satellite_S5P_202411.log -├── processed -│ ├── NO2 -│ │ ├── 2023 -│ │ │ ├── 01 -│ │ │ └── ... -│ │ └── 2024 -│ │ ├── 01 -│ │ └── ... -│ ├── SO2 -│ │ ├── 2023 -│ │ │ ├── 01 -│ │ │ └── ... -│ │ └── 2024 -│ │ ├── 01 -│ │ └── ... -│ └── ... -└── raw - ├── NO2 - │ ├── 2023 - │ │ ├── 01 - │ │ └── ... - │ └── 2024 - │ ├── 01 - │ └── ... - ├── SO2 - │ ├── 2023 - │ │ ├── 01 - │ │ └── ... - │ └── 2024 - │ ├── 01 - │ └── ... - └── ... - -""" diff --git a/src/processing/data_processor.py b/src/processing/data_processor.py index 4633743..a797298 100644 --- a/src/processing/data_processor.py +++ b/src/processing/data_processor.py @@ -1,17 +1,16 @@ """src/processing/no2_processor.py""" +import logging import numpy as np import xarray as xr -import logging -from pathlib import Path from datetime import datetime from dateutil.relativedelta import relativedelta from src.processing.interpolators import DataInterpolator from src.processing.taiwan_frame import TaiwanFrame from src.config.settings import RAW_DATA_DIR, PROCESSED_DATA_DIR, FIGURE_DIR, FIGURE_BOUNDARY -from src.visualization.sample_plot_nc import plot_global_var -from src.utils.catalog import ProductClassInput, ProductTypeInput, ProductType -from src.utils.richer import DisplayManager +from src.visualization.plot_nc import plot_global_var +from src.config.catalog import ClassInput, TypeInput, PRODUCT_CONFIGS +from src.config.richer import DisplayManager logger = logging.getLogger(__name__) @@ -48,82 +47,58 @@ def create_grid(self, lon: np.ndarray, lat: np.ndarray): # 創建網格矩陣 return np.meshgrid(grid_lon, grid_lat) - def extract_data(self, dataset: xr.Dataset, use_taiwan_mask: bool = False): - """從數據集中提取數據 + def extract_data(self, dataset: xr.Dataset, extract_range: tuple[float, float, float, float] = None): + """提取數據，可選擇是否限定範圍 - Parameters: - ----------- - dataset : xr.Dataset - 輸入的數據集 - use_taiwan_mask : bool - 是否只提取台灣區域的數據 + Args: + dataset: xarray Dataset + extract_range: 可選的tuple (min_lon, max_lon, min_lat, max_lat)，如果提供則提取指定範圍 """ - if use_taiwan_mask: - return self._extract_data_taiwan(dataset) - return self._extract_data_global(dataset) - - def _extract_data_global(self, dataset: xr.Dataset): - """提取全球範圍的數據""" - time = np.datetime64(dataset.time.values[0], 'D') - lat = dataset.latitude.values[0] - lon = dataset.longitude.values[0] - shape = dataset.longitude.values[0].shape - var = dataset.nitrogendioxide_tropospheric_column.values[0] - qa = dataset.qa_value.values[0] - - # 根據 QA 值過濾數據 - mask = qa < self.mask_qc_value - var[mask] = np.nan - - info_dict = { - 'time': f"{time}", - 'shape': f"{shape}", - 'latitude': f'{lat.min():.2f} to {lat.max():.2f}', - 'longitude': f'{lon.min():.2f} to {lon.max():.2f}', - } - - return (lon, lat, var), info_dict - - def _extract_data_taiwan(self, dataset: xr.Dataset, extract_range: tuple[float, float, float, float] = None): - """提取台灣區域的數據""" - min_lon, max_lon, min_lat, max_lat = FIGURE_BOUNDARY - # 設定條件 + # 初始處理 time = np.datetime64(dataset.time.values[0], 'D') - mask_lon = ((dataset.longitude >= min_lon) & (dataset.longitude <= max_lon)) - mask_lat = ((dataset.latitude >= min_lat) & (dataset.latitude <= max_lat)) - masked_lon_lat_ds = dataset.where((mask_lon & mask_lat), drop=True) + attributes = PRODUCT_CONFIGS[self.product_type].dataset_name - if masked_lon_lat_ds.sizes['scanline'] == 0 or masked_lon_lat_ds.sizes['ground_pixel'] == 0: - raise ValueError("No data points within Taiwan region") + # 如果提供了範圍，進行過濾 + if extract_range is not None: + min_lon, max_lon, min_lat, max_lat = extract_range + mask_lon = (dataset.longitude >= min_lon) & (dataset.longitude <= max_lon) + mask_lat = (dataset.latitude >= min_lat) & (dataset.latitude <= max_lat) + dataset = dataset.where((mask_lon & mask_lat), drop=True) - mask_qa = (masked_lon_lat_ds.qa_value >= self.mask_qc_value) - masked_ds = masked_lon_lat_ds.where(mask_qa) + # 檢查是否有數據 + if dataset.sizes['scanline'] == 0 or dataset.sizes['ground_pixel'] == 0: + raise ValueError(f"No data points within region: {extract_range}") - attributes = self.product_type.dataset_name + # QA 過濾 + mask_qa = (dataset.qa_value >= self.mask_qc_value) + dataset = dataset.where(mask_qa) - if np.all(np.isnan(masked_ds[attributes])): + # 檢查數據有效性 + if np.all(np.isnan(dataset[attributes])): raise ValueError("No valid data points after QA filtering") - lon = masked_ds.longitude[0].values - lat = masked_ds.latitude[0].values - shape = masked_ds.latitude[0].values.shape - var = masked_ds[attributes][0].values + lon = dataset.longitude[0].values + lat = dataset.latitude[0].values + shape = dataset.latitude[0].values.shape + var = dataset[attributes][0].values - self.nc_info.update({ + info_dict = { 'time': f"{time}", 'shape': f"{shape}", 'latitude': f'{np.nanmin(lat):.2f} to {np.nanmax(lat):.2f}', 'longitude': f'{np.nanmin(lon):.2f} to {np.nanmax(lon):.2f}', - }) + } + + if hasattr(self, 'nc_info'): + self.nc_info.update(info_dict) return lon, lat, var def process_each_data(self, - file_class: ProductClassInput, - file_type: ProductTypeInput, + file_class: ClassInput, + file_type: TypeInput, start_date: str, end_date: str, - use_taiwan_mask: bool = False, ): """ 處理指定日期範圍內的衛星數據 @@ -133,66 +108,57 @@ def process_each_data(self, file_type: ProductTypeInput start_date (str): 開始日期 (YYYY-MM-DD) end_date (str): 結束日期 (YYYY-MM-DD) - use_taiwan_mask (bool): 是否使用台灣遮罩 """ - self.product_type = ProductType[file_type] - - def process_single_file(file_path, output_dir, figure_dir): + def process_single_file(file_path, output_dir): """處理單一數據檔案""" - print(file_path) - ds = xr.open_dataset(file_path, group='PRODUCT', engine='netcdf4') - breakpoint() + ds = xr.open_dataset(file_path, engine='netcdf4', group='PRODUCT') + # 確保輸出目錄存在 output_dir.mkdir(parents=True, exist_ok=True) - ds.to_netcdf(f'{output_dir / file_path.stem}.nc') + output_path = output_dir / file_path.name + + if not output_path.exists(): + ds.to_netcdf(output_path) if ds is None: return try: - # 0. 紀錄 nc 檔訊息 + # 1. 紀錄 nc 檔訊息 self.nc_info = {'file_name': file_path.name} - # 1. 提取數據和信息 - lon, lat, var = self.extract_data(ds, use_taiwan_mask) - - # 2. 顯示檔案信息 - DisplayManager().display_product_info(self.nc_info) - - # 3. 創建網格並進行插值 - lon_grid, lat_grid = self.create_grid(lon, lat) - var_grid = DataInterpolator.interpolate( - lon, lat, var, - lon_grid, lat_grid, - method=self.interpolation_method - ) - - # 4. 創建插值後的數據集 - interpolated_ds = xr.Dataset( - { - self.product_type.dataset_name: ( - ['time', 'latitude', 'longitude'], - var_grid[np.newaxis, :, :] - ) - }, - coords={ - 'time': ds.time.values[0:1], - 'latitude': np.squeeze(lat_grid[:, 0]), - 'longitude': np.squeeze(lon_grid[0, :]) - } - ) - - # 5. 繪製並保存圖形 - plot_global_var( - data=interpolated_ds, - product_params=self.product_type, - savefig_path=figure_dir / file_path.stem, - map_scale='Taiwan' - ) + # # 2. 提取數據和信息 + # lon, lat, var = self.extract_data(ds, extract_range=FIGURE_BOUNDARY) + # + # # 3. 顯示檔案信息 + # DisplayManager().display_product_info(self.nc_info) + # + # # 4. 創建網格並進行插值 + # lon_grid, lat_grid = self.create_grid(lon, lat) + # var_grid = DataInterpolator.interpolate( + # lon, lat, var, + # lon_grid, lat_grid, + # method=self.interpolation_method + # ) + # + # # 5. 創建插值後的數據集 + # interpolated_ds = xr.Dataset( + # { + # self.product_type.dataset_name: ( + # ['time', 'latitude', 'longitude'], + # var_grid[np.newaxis, :, :] + # ) + # }, + # coords={ + # 'time': ds.time.values[0:1], + # 'latitude': np.squeeze(lat_grid[:, 0]), + # 'longitude': np.squeeze(lon_grid[0, :]) + # } + # ) finally: ds.close() - # breakpoint() + # 主處理流程 start = datetime.strptime(start_date, '%Y-%m-%d') end = datetime.strptime(end_date, '%Y-%m-%d') @@ -201,7 +167,7 @@ def process_single_file(file_path, output_dir, figure_dir): # 按月份逐月處理 while current_date <= end: # 1. 準備當月的目錄路徑 - product_type = ProductType[file_type] + product_type = file_type year = current_date.strftime('%Y') month = current_date.strftime('%m') @@ -210,24 +176,45 @@ def process_single_file(file_path, output_dir, figure_dir): figure_dir = FIGURE_DIR / product_type / year / month # 2. 創建必要的目錄 - for directory in [input_dir, output_dir, figure_dir]: + for directory in [output_dir, figure_dir]: directory.mkdir(parents=True, exist_ok=True) - # 3. 選擇處理類型 離線 | 即時 | 物種 - file_class = '' if file_class == '*' else file_class - file_type = '' if file_type == '*' else file_type - - filter_pattern = f"*{file_class}_L2__{file_type}*.nc" - # breakpoint() - for file_path in input_dir.glob(filter_pattern): - # 檢查檔案日期是否在指定範圍內 + file_pattern = f"*{file_class}_L2__{file_type}*.nc" + + # 3. 處理原始數據（如果存在） + if input_dir.exists(): + for file_path in input_dir.glob(file_pattern): + # 檢查檔案日期是否在指定範圍內 + date_to_check = datetime.strptime(file_path.name[20:28], '%Y%m%d') + if not (start <= date_to_check <= end): + continue + try: + process_single_file(file_path, output_dir) + except Exception as e: + logger.error(f"處理檔案 {file_path.name} 時發生錯誤: {e}") + continue + + # 4. 繪製圖片（使用處理後的數據） + processed_files = list(output_dir.glob(file_pattern)) + if not processed_files: + logger.warning(f"在 {output_dir} 中找不到符合條件的處理後檔案") + continue + + for file_path in processed_files: date_to_check = datetime.strptime(file_path.name[20:28], '%Y%m%d') if not (start <= date_to_check <= end): continue + figure_path = figure_dir / f"{file_path.stem}.png" try: - process_single_file(file_path, output_dir, figure_dir) + plot_global_var( + dataset=file_path, + product_params=PRODUCT_CONFIGS[file_type], + savefig_path=figure_path, + map_scale='Taiwan', + show_stations=True + ) except Exception as e: - logger.error(f"處理檔案 {file_path.name} 時發生錯誤: {e}") + logger.error(f"繪製檔案 {file_path.name} 時發生錯誤: {e}") continue # 4. 移至下個月 @@ -269,9 +256,3 @@ def _save_monthly_average(container, grid, year, month, output_file): # 確保輸出目錄存在 output_file.parent.mkdir(parents=True, exist_ok=True) ds_result.to_netcdf(output_file) - - logger.info(f"Saved monthly average to {output_file}") - logger.info(f"Final grid shape: {no2_average.shape}") - logger.info(f"Time: {ds_result.time.values}") - logger.info(f"Longitude range: {grid[0][0].min():.2f} to {grid[0][0].max():.2f}") - logger.info(f"Latitude range: {grid[1][:, 0].min():.2f} to {grid[1][:, 0].max():.2f}") \ No newline at end of file diff --git a/src/utils/catalog.py b/src/utils/catalog.py deleted file mode 100644 index 616c011..0000000 --- a/src/utils/catalog.py +++ /dev/null @@ -1,199 +0,0 @@ -from dataclasses import dataclass -from typing import Literal -from enum import Enum - - -@dataclass -class ProductConfig: - """產品配置""" - display_name: str # 顯示名稱（帶下標） - dataset_name: str # 數據集名稱 - vmin: float | None # 最小值 - vmax: float | None # 最大值 - units: str # 單位 - title: str # 標題 - cmap: str = 'viridis' # 預設色階 - - -class ProductType(str, Enum): - """Available product types""" - # Ozone products - O3 = 'O3' - O3____ = 'O3____' - O3_TCL = 'O3_TCL' # No data ?? - O3__PR = 'O3__PR' # start from 2022 - # Other gas products - CH4___ = 'CH4___' - CO____ = 'CO____' - NO2___ = 'NO2___' - HCHO__ = 'HCHO__' - SO2___ = 'SO2___' - # Cloud and aerosol products - CLOUD_ = 'CLOUD_' - FRESCO = 'FRESCO' - AER_LH = 'AER_LH' - AER_AI = 'AER_AI' - - @property - def config(self) -> ProductConfig: - """獲取產品配置""" - configs = { - 'O3____': ProductConfig( - display_name='O\u2083', - dataset_name='ozone_total_vertical_column', - vmin=0.1, - vmax=0.2, - units=f'O$_3$ Total Vertical Column (mol/m$^2$)', - title=f'O$_3$ Total Vertical Column', - cmap='RdBu_r' - ), - 'O3_TCL': ProductConfig( - display_name='O\u2083_TCL', - dataset_name='ozone_tropospheric_vertical_column', - vmin=0.1, - vmax=0.2, - units=f'O$_3$ Tropospheric Column (mol/m$^2$)', - title=f'O$_3$ Tropospheric Column', - cmap='RdBu_r' - ), - 'O3__PR': ProductConfig( - display_name='O\u2083_PR', - dataset_name='ozone_profile', - vmin=0.1, - vmax=0.2, - units=f'O$_3$ Profile (mol/m$^2$)', - title=f'O$_3$ Profile', - cmap='RdBu_r' - ), - 'CH4___': ProductConfig( - display_name='CH\u2084', - dataset_name='methane_mixing_ratio', - vmin=None, - vmax=None, - units=f'CH$_4$ Methane Mixing Ratio', - title=f'CH$_4$ Methane Mixing Ratio', - cmap='RdBu_r' - ), - 'CO____': ProductConfig( - display_name='CO', - dataset_name='carbonmonoxide_total_column', - vmin=None, - vmax=None, - units=f'CO Total Column (mol/m$^2$)', - title=f'CO Total Column', - cmap='RdBu_r' - ), - 'NO2___': ProductConfig( - display_name='NO\u2082', - dataset_name='nitrogendioxide_tropospheric_column', - vmin=-2e-4, - vmax=2e-4, - units=f'NO$_2$ Tropospheric Column (mol/m$^2$)', - title=f'NO$_2$ Tropospheric Column', - cmap='RdBu_r' - ), - 'HCHO__': ProductConfig( - display_name='HCHO', - dataset_name='formaldehyde_tropospheric_vertical_column', - vmin=-4e-4, - vmax=4e-4, - units=f'HCHO Tropospheric Vertical Column (mol/m$^2$)', - title=f'HCHO Tropospheric Vertical Column', - cmap='RdBu_r' - ), - 'SO2___': ProductConfig( - display_name='SO\u2082', - dataset_name='sulfurdioxide_total_vertical_column', - vmin=-4e-3, - vmax=4e-3, - units=f'SO$_2$ Total Vertical Column (mol/m$^2$)', - title=f'SO$_2$ Total Vertical Column', - cmap='RdBu_r' - ), - 'AER_AI': ProductConfig( - display_name='Aerosol Index', - dataset_name='aerosol_index_340_380', - vmin=None, - vmax=None, - units=f'Aerosol Index', - title=f'Aerosol Index', - cmap='viridis' - ) - - } - return configs.get(self.value, ProductConfig( - display_name=self.value, - dataset_name='unknown_dataset', - vmin=0, - vmax=100, - units='', - title='unknown_dataset', - cmap='RdBu_r' - )) - - @property - def display_name(self) -> str: - return self.config.display_name - - @property - def dataset_name(self) -> str: - return self.config.dataset_name - - @property - def vmin(self) -> float: - return self.config.vmin - - @property - def vmax(self) -> float: - return self.config.vmax - - @property - def units(self) -> str: - return self.config.units - - @property - def title(self) -> str: - return self.config.title - - @property - def cmap(self) -> str: - return self.config.cmap - - -class ProductLevel(str, Enum): - Level0 = 'L0__' - Level1B = 'L1B_' - Level2 = 'L2__' - - -class ProductClass(str, Enum): - """Available processing classes: - - NRTI: near-real time processing - - OFFL: offline processing - - RPRO: reprocessing - - ALL(*): match all classes - Other classes (mainly for testing): - - TEST: internal testing - - OGCA: on-ground calibration - - GSOV: ground segment validation - - OPER: operational processing - """ - TEST = 'TEST' # internal testing - OGCA = 'OGCA' # on-ground calibration - GSOV = 'GSOV' # ground segment overall validation - OPER = 'OPER' # operational processing - NRTI = 'NRTI' # near-real time processing - OFFL = 'OFFL' # offline processing - RPRO = 'RPRO' # reprocessing - ALL = '*' # match all - - -ProductTypeInput = Literal[ - 'O3', 'O3____', 'O3_TCL', 'O3__PR', # Ozone products - 'CH4___', 'CO____', 'NO2___', 'HCHO__', 'SO2___', # Other gas products - 'CLOUD_', 'FRESCO', 'AER_LH', 'AER_AI', # Cloud and aerosol products -] - -ProductClassInput = Literal[ - 'NRTI', 'OFFL', 'RPRO', '*' -] diff --git a/src/utils/logger.py b/src/utils/logger.py deleted file mode 100644 index f548628..0000000 --- a/src/utils/logger.py +++ /dev/null @@ -1,26 +0,0 @@ -"""日誌配置""" -import logging -from datetime import datetime -from pathlib import Path -from src.config.settings import LOGS_DIR - - -def setup_logging(): - """設置日誌配置""" - # 確保日誌目錄存在 - log_dir = Path(LOGS_DIR) - log_dir.mkdir(parents=True, exist_ok=True) - - # 創建日誌檔案路徑 - log_file = log_dir / f"Satellite_S5P_{datetime.now().strftime('%Y%m')}.log" - - # 配置基本設定 - logging.basicConfig( - level=logging.INFO, - format='%(asctime)s - %(name).10s - %(message)s', - datefmt='%Y-%m-%d %H:%M:%S', - handlers=[ - logging.FileHandler(log_file, encoding='utf-8'), - logging.StreamHandler() # 同時輸出到控制台 - ] - ) diff --git a/src/visualization/gif_nc.py b/src/visualization/gif_nc.py index ec4c4ab..844a199 100644 --- a/src/visualization/gif_nc.py +++ b/src/visualization/gif_nc.py @@ -1,14 +1,13 @@ import imageio from datetime import datetime -from pathlib import Path import re from PIL import Image import numpy as np from src.config.settings import FIGURE_DIR -from src.utils.catalog import ProductTypeInput +from src.config.catalog import TypeInput -def animate_data(file_type: ProductTypeInput, start_date, end_date, fps=1, resize=None, **kwargs): +def animate_data(file_type: TypeInput, start_date, end_date, fps=1, resize=None, **kwargs): """ 將 Sentinel 衛星圖片製作成 GIF 動畫 diff --git a/src/visualization/plot_nc.py b/src/visualization/plot_nc.py index add2996..33fe3e2 100644 --- a/src/visualization/plot_nc.py +++ b/src/visualization/plot_nc.py @@ -1,24 +1,230 @@ -import cartopy.crs as ccrs -import numpy as np +import logging +import geopandas as gpd import xarray as xr import matplotlib.pyplot as plt -import cartopy.mpl.gridliner as gridliner -import time -import geopandas as gpd - -from scipy.signal import convolve2d -from netCDF4 import Dataset -from matplotlib.ticker import ScalarFormatter +import cartopy.crs as ccrs +import cartopy.feature as cfeature +import numpy as np +from typing import Literal from pathlib import Path -from src.processing.taiwan_frame import TaiwanFrame +from shapely.geometry import Point +from scipy.signal import convolve2d +from matplotlib.ticker import ScalarFormatter, FixedLocator + +from src.config.settings import FIGURE_BOUNDARY +from src.config.richer import DisplayManager + +plt.rcParams['mathtext.fontset'] = 'custom' +plt.rcParams['mathtext.rm'] = 'Times New Roman' +plt.rcParams['mathtext.it'] = 'Times New Roman: italic' +plt.rcParams['mathtext.bf'] = 'Times New Roman: bold' +plt.rcParams['mathtext.default'] = 'regular' +plt.rcParams['font.family'] = 'Times New Roman' +plt.rcParams['font.weight'] = 'normal' +plt.rcParams['font.size'] = 16 -taiwan_counties = gpd.read_file(Path(__file__).parent / "taiwan/COUNTY_MOI_1090820.shp") -station = gpd.read_file(Path(__file__).parent / "stations/空氣品質監測站位置圖_121_10704.shp") +plt.rcParams['axes.titlesize'] = 'large' +plt.rcParams['axes.titleweight'] = 'bold' +plt.rcParams['axes.labelweight'] = 'bold' -# @setFigure -def platecarree_plot(dataset, zoom=True, path=None, **kwargs): +logger = logging.getLogger(__name__) + + +def smooth_kernel(data, kernel_size=5): + kernel = np.ones((kernel_size, kernel_size)) + return convolve2d(data, kernel, mode='same', boundary='wrap') / np.sum(kernel) + + +def plot_stations(ax, stations: list[str], label_offset: tuple[float, float] = (-0.2, 0)): + """繪製測站標記和標籤 + + Args: + ax: matplotlib axes + stations: 要標記的測站名稱列表 + label_offset: (x偏移, y偏移)，用於調整標籤位置 + """ + # 讀取測站資料 + station_data = gpd.read_file( + Path(__file__).parents[2] / "data/shapefiles/stations/空氣品質監測站位置圖_121_10704.shp" + ) + # 創建測站的 GeoDataFrame + station_geometry = [Point(xy) for xy in zip(station_data['TWD97Lon'], station_data['TWD97Lat'])] + + # all station + # geodata = gpd.GeoDataFrame(station_data, crs=ccrs.PlateCarree(), geometry=station_geometry) + # geodata.plot(ax=ax, color='gray', markersize=10) + + # matplotlib 預設顏色和標記循環 + colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] # 預設顏色循環 + markers = ['o', 's', '^', 'v', 'D', '<', '>', 'p', '*'] # 標記循環 + + # 過濾出要標記的測站 + target_stations = station_data[station_data['SiteName'].isin(stations)] + + legend_labels = [] + legend_handles = [] + + for i, (_, row) in enumerate(target_stations.iterrows()): + # 循環使用顏色和標記 + color = colors[i % len(colors)] + marker = markers[i % len(markers)] + + # 畫點並保存標記物件 + marker_obj = ax.scatter(row['TWD97Lon'], row['TWD97Lat'], + marker=marker, + color=color, + s=15, + transform=ccrs.PlateCarree(), + label=row['SiteEngNam']) + + legend_labels.append(row['SiteEngNam']) + legend_handles.append(marker_obj) + + # 添加圖例 + ax.legend(handles=legend_handles, + labels=legend_labels, + loc='upper right', + bbox_to_anchor=(1, 1), + borderaxespad=0.) + + +def plot_global_var(dataset: Path | str, + product_params, + show_info: bool = True, + savefig_path=None, + map_scale: Literal['global', 'Taiwan'] = 'global', + show_stations: bool = False, + mark_stations: list = ['古亭', '楠梓', '鳳山'], + ): + """ + 在全球地圖上繪製 var 分布圖 + """ + try: + # 判斷輸入類型並適當處理 + from netCDF4 import Dataset + if isinstance(dataset, (str, Path)): + if 'PRODUCT' in Dataset(dataset, 'r').groups: + ds = xr.open_dataset(dataset, engine='netcdf4', group='PRODUCT') + else: + ds = xr.open_dataset(dataset, engine='netcdf4') + elif isinstance(dataset, xr.Dataset): + ds = dataset + else: + raise NotImplementedError + + if show_info: + lon = ds.longitude[0].values + lat = ds.latitude[0].values + var = ds[product_params.dataset_name][0].values + + nc_info = {'file_name': dataset.name if isinstance(dataset, Path) else '', + 'time': np.datetime64(ds.time.values[0], 'D'), + 'shape': ds.latitude[0].values.shape, + 'latitude': f'{np.nanmin(lat):.2f} to {np.nanmax(lat):.2f}', + 'longitude': f'{np.nanmin(lon):.2f} to {np.nanmax(lon):.2f}', + } + + DisplayManager().display_product_info(nc_info) + + # 創建圖形和投影 + fig = plt.figure(figsize=(12, 8) if map_scale == 'global' else (8, 8)) + ax = plt.axes(projection=ccrs.PlateCarree()) + + # 設定全球範圍 + if map_scale == 'global': + ax.set_global() + else: + ax.set_extent(FIGURE_BOUNDARY, crs=ccrs.PlateCarree()) + + # 繪製數據 + dataset = ds[product_params.dataset_name][0] + dataset.data = smooth_kernel(dataset.data, kernel_size=3) + + # 方法1：使用 ScalarFormatter + formatter = ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((-2, 2)) # 可以調整使用科學記號的範圍 + + im = dataset.plot( + ax=ax, + x='longitude', y='latitude', + cmap='RdBu_r', + transform=ccrs.PlateCarree(), + robust=True, # 自動處理極端值 + vmin=product_params.vmin, + vmax=product_params.vmax, + cbar_kwargs={ + 'label': product_params.units, + 'fraction': 0.04, # colorbar 的寬度 (預設是 0.15) + 'pad': 0.04, # colorbar 和圖之間的間距 + 'aspect': 30, # colorbar 的長寬比，增加這個值會讓 colorbar 變長 + 'shrink': 1, + 'format': formatter, + 'extend': 'neither', + } + ) + + # plot = dataset.plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=False, cmap='jet') + + # 如果是台灣範圍且需要顯示測站 + if map_scale == 'Taiwan': + # ax.add_feature(cfeature.COASTLINE.with_scale('10m')) + + # 讀取縣市和測站資料並添加縣市邊界 + taiwan_counties = gpd.read_file(Path(__file__).parents[2] / "data/shapefiles/taiwan/COUNTY_MOI_1090820.shp") + ax.add_geometries(taiwan_counties['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none') + + if show_stations and mark_stations: + plot_stations(ax, mark_stations) + + + else: + # 添加地圖特徵 + ax.add_feature(cfeature.BORDERS.with_scale('10m'), linestyle=':') + ax.add_feature(cfeature.COASTLINE.with_scale('10m')) + ax.add_feature(cfeature.LAND.with_scale('10m'), alpha=0.1) + ax.add_feature(cfeature.OCEAN.with_scale('10m'), alpha=0.1) + + # 設定網格線 + gl = ax.gridlines(draw_labels=True, linestyle='--', alpha=0.7) + gl.top_labels = False + gl.right_labels = False + gl.xlocator = FixedLocator([119, 120, 121, 122, 123]) # 設定經度刻度 + + # 用矩形標記數據範圍 + lon_min, lon_max = float(np.nanmin(ds.longitude)), float(np.nanmax(ds.longitude)) + lat_min, lat_max = float(np.nanmin(ds.latitude)), float(np.nanmax(ds.latitude)) + rect = plt.Rectangle( + (lon_min, lat_min), + lon_max - lon_min, + lat_max - lat_min, + fill=False, + color='red', + transform=ccrs.PlateCarree(), + linewidth=2 + ) + ax.add_patch(rect) + + # 設定標題 + time_str = np.datetime64(dataset.time.values, 'D').astype(str) + plt.title(f'{product_params.title} {time_str}', pad=20, fontdict={'weight': 'bold', 'fontsize': 24}) + + plt.tight_layout() + plt.show() + + if savefig_path is not None: + fig.savefig(savefig_path, dpi=600) + + ds.close() + + except Exception as e: + logger.error(f"繪圖時發生錯誤: {str(e)}") + raise + + +def platecarree_plot(dataset, product_params, zoom=True, path=None, **kwargs): fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={'projection': ccrs.PlateCarree()}) if zoom: @@ -27,16 +233,12 @@ def platecarree_plot(dataset, zoom=True, path=None, **kwargs): gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5) gl.xlocator = plt.FixedLocator([119, 120, 121, 122, 123]) gl.top_labels = gl.right_labels = False - # ax.coastlines(resolution='10m', color='gray', linewidth=1) - ax.add_geometries(taiwan_counties['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=1) else: ax.set_global() gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5) gl.top_labels = gl.right_labels = False - # ax.coastlines(resolution='10m', color='black', linewidth=0.5) - ax.add_geometries(taiwan_counties['geometry'], resolution='10m', crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none') - var = 'nitrogendioxide_tropospheric_column' + var = product_params.dataset_name plot = dataset[var][0].plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=False, cmap='jet', vmin=0, vmax=1.4e-4) @@ -44,9 +246,9 @@ def platecarree_plot(dataset, zoom=True, path=None, **kwargs): cbar.set_label(r'$\bf NO_{2}\ mole/m^2$') # 设置colorbar刻度标签格式为科学记数法 - cbar.formatter = ScalarFormatter(useMathText=True, useOffset=True) - cbar.formatter.set_powerlimits((-2, 2)) # 设置指数的显示范围 - cbar.update_ticks() + formatter = ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((-2, 2)) # 可以調整使用科學記號的範圍 plt.title(kwargs.get('title')) if path is not None: @@ -54,11 +256,11 @@ def platecarree_plot(dataset, zoom=True, path=None, **kwargs): plt.show() -def orthographic_plot(dataset): +def orthographic_plot(dataset, product_params): projection = ccrs.Orthographic(120, 25) fig, ax = plt.subplots(figsize=(8, 6), subplot_kw={'projection': projection}) - var = 'nitrogendioxide_tropospheric_column' + var = product_params.dataset_name vmin = dataset[var][0].min() dataset[var][0].plot.pcolormesh(ax=ax, x='longitude', y='latitude', @@ -68,61 +270,15 @@ def orthographic_plot(dataset): ax.set_global() ax.coastlines(resolution='10m') + plt.show() -def smooth_kernel(data, kernel_size=5): - """ - - :param data: - :param kernel_size: - :return: - """ - kernel = np.ones((kernel_size, kernel_size)) - - return convolve2d(data, kernel, mode='same', boundary='wrap') / np.sum(kernel) - - -if __name__ == '__main__': - year = '2022' # 改年度即可使用 - folder_path = Path(__file__).parent / "ncData" - output_folder = Path(__file__).parent / "Figure" - - lon_coordinate, lat_coordinate = TaiwanFrame().frame() - container = [] - - # 遍历文件夹中符合条件的文件 - for file_path in folder_path.glob(f"NO2_{year}*.nc"): - dataset = xr.open_dataset(file_path) - - # 提取数据并进行平滑处理 - traget_data = dataset.nitrogendioxide_tropospheric_column[0].data - container.append(traget_data) - - dataset.nitrogendioxide_tropospheric_column[0] = smooth_kernel(traget_data) - - # 生成图像文件并保存到输出文件夹 - # print('plot: ' + file_path.name) - # output_path = output_folder / f"{file_path.stem}.png" - # platecarree_plot(dataset, path=output_path, - # title=rf'$\bf {file_path.stem[4:8]}-{file_path.stem[8:10]}\ \ NO_{2}$') - - - # 计算年平均值 - no2_year_average = np.nanmean(container, axis=0) - original_shape = no2_year_average.shape - - ds_time = np.array(np.datetime64(year, 'ns')) - ds_result = xr.Dataset( - {'nitrogendioxide_tropospheric_column': (['latitude', 'longitude'], no2_year_average.reshape(*original_shape))}, - coords={'latitude': lat_coordinate[:, 0], - 'longitude': lon_coordinate[0]}) - - ds_result = ds_result.expand_dims(time=[ds_time]) +if __name__ == "__main__": + file_group = '/Users/chanchihyu/Sentinel-5P/raw/NO2___/2024/03/S5P_OFFL_L2__NO2____20240314T031823_20240314T045953_33252_03_020600_20240315T192700.nc' + file_ungroup = '/Users/chanchihyu/Sentinel-5P/processed/NO2___/2024/01/S5P_OFFL_L2__NO2____20240110T045402_20240110T063532_32345_03_020600_20240111T211523.nc' + ds = xr.open_dataset(file_ungroup) - # 畫年平均 - year_ds = ds_result.copy(deep=True) - traget_data = year_ds.nitrogendioxide_tropospheric_column[0].data - year_ds['nitrogendioxide_tropospheric_column'][0] = smooth_kernel(traget_data) + from netCDF4 import Dataset + nc = Dataset(file_ungroup, 'r').groups - output_path = output_folder / f"NO2_{year}.png" - platecarree_plot(year_ds, path=output_path, title=rf'$\bf {year}\ \ NO_{2}$') + # plot_global_var(file) diff --git a/src/visualization/read_shp.py b/src/visualization/read_shp.py deleted file mode 100644 index 3b5446d..0000000 --- a/src/visualization/read_shp.py +++ /dev/null @@ -1,23 +0,0 @@ -import matplotlib.pyplot as plt -import cartopy.crs as ccrs -import pandas as pd -import geopandas as gpd -from pathlib import Path -from shapely.geometry import Point, Polygon - -taiwan_counties = gpd.read_file(Path(__file__).parents[2] / "data/shapefiles/taiwan/COUNTY_MOI_1090820.shp") -station = gpd.read_file(Path(__file__).parents[2] / "data/shapefiles/stations/空氣品質監測站位置圖_121_10704.shp") - -geometry = [Point(xy) for xy in zip(station['TWD97Lon'], station['TWD97Lat'])] -geodata = gpd.GeoDataFrame(station, crs=ccrs.PlateCarree(), geometry=geometry) - -# 创建地图 -fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()}) -ax.set_extent([119, 123, 21, 26], crs=ccrs.PlateCarree()) - -# 添加县市边界 -ax.add_geometries(taiwan_counties['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none') -ax.add_geometries(station['geometry'], crs=ccrs.PlateCarree(), edgecolor='Red', facecolor='none') -geodata.plot(ax=ax, color='red', markersize=5) - -plt.show() diff --git a/src/visualization/sample_plot_nc.py b/src/visualization/sample_plot_nc.py deleted file mode 100644 index b57f2df..0000000 --- a/src/visualization/sample_plot_nc.py +++ /dev/null @@ -1,126 +0,0 @@ -import xarray as xr -import matplotlib.pyplot as plt -import cartopy.crs as ccrs -import cartopy.feature as cfeature -import numpy as np -from typing import Literal -from pathlib import Path -import logging -from matplotlib.ticker import ScalarFormatter, FuncFormatter, FixedLocator -from src.config.settings import FILTER_BOUNDARY, FIGURE_BOUNDARY - - -plt.rcParams['mathtext.fontset'] = 'custom' -plt.rcParams['mathtext.rm'] = 'Times New Roman' -plt.rcParams['mathtext.it'] = 'Times New Roman: italic' -plt.rcParams['mathtext.bf'] = 'Times New Roman: bold' -plt.rcParams['mathtext.default'] = 'regular' -plt.rcParams['font.family'] = 'Times New Roman' -plt.rcParams['font.weight'] = 'normal' -plt.rcParams['font.size'] = 16 - -plt.rcParams['axes.titlesize'] = 'large' -plt.rcParams['axes.titleweight'] = 'bold' -plt.rcParams['axes.labelweight'] = 'bold' - - -logger = logging.getLogger(__name__) - - -def plot_global_var(data, - product_params, - savefig_path=None, - map_scale: Literal['global', 'Taiwan'] = 'global' - ): - """ - 在全球地圖上繪製 var 分布圖 - """ - try: - # 判斷輸入類型並適當處理 - ds = xr.open_dataset(data) if isinstance(data, (str, Path)) else data - - # 創建圖形和投影 - fig = plt.figure(figsize=(12, 8) if map_scale == 'global' else (8, 8)) - ax = plt.axes(projection=ccrs.PlateCarree()) - - # 設定全球範圍 - if map_scale == 'global': - ax.set_global() - else: - ax.set_extent(FIGURE_BOUNDARY, crs=ccrs.PlateCarree()) - - # 繪製數據 - data = ds[product_params.dataset_name][0] - - # 方法1：使用 ScalarFormatter - formatter = ScalarFormatter(useMathText=True) - formatter.set_scientific(True) - formatter.set_powerlimits((-2, 2)) # 可以調整使用科學記號的範圍 - - im = data.plot( - ax=ax, - cmap='RdBu_r', - transform=ccrs.PlateCarree(), - robust=True, # 自動處理極端值 - vmin=product_params.vmin, - vmax=product_params.vmax, - cbar_kwargs={ - 'label': product_params.units, - 'fraction': 0.04, # colorbar 的寬度 (預設是 0.15) - 'pad': 0.04, # colorbar 和圖之間的間距 - 'aspect': 30, # colorbar 的長寬比，增加這個值會讓 colorbar 變長 - 'shrink': 1, - 'format': formatter, - 'extend': 'neither', - } - ) - - # 添加地圖特徵 - ax.add_feature(cfeature.BORDERS.with_scale('10m'), linestyle=':') - ax.add_feature(cfeature.COASTLINE.with_scale('10m')) - ax.add_feature(cfeature.LAND.with_scale('10m'), alpha=0.1) - ax.add_feature(cfeature.OCEAN.with_scale('10m'), alpha=0.1) - - # 設定網格線 - gl = ax.gridlines(draw_labels=True, linestyle='--', alpha=0.7) - gl.top_labels = False - gl.right_labels = False - gl.xlocator = FixedLocator([119, 120, 121, 122, 123]) # 設定經度刻度 - - # 用矩形標記數據範圍 - lon_min, lon_max = float(np.nanmin(ds.longitude)), float(np.nanmax(ds.longitude)) - lat_min, lat_max = float(np.nanmin(ds.latitude)), float(np.nanmax(ds.latitude)) - rect = plt.Rectangle( - (lon_min, lat_min), - lon_max - lon_min, - lat_max - lat_min, - fill=False, - color='red', - transform=ccrs.PlateCarree(), - linewidth=2 - ) - ax.add_patch(rect) - - # 設定標題 - time_str = np.datetime64(data.time.values, 'D').astype(str) - plt.title(f'{product_params.title} {time_str}', pad=20, fontdict={'weight': 'bold', 'fontsize': 24}) - - plt.tight_layout() - plt.show() - - if savefig_path is not None: - fig.savefig(savefig_path, dpi=600) - - ds.close() - - except Exception as e: - logger.error(f"繪圖時發生錯誤: {str(e)}") - raise - - -# 主程式 -if __name__ == "__main__": - file_list = ["/Users/chanchihyu/Sentinel_data/raw/2024/04/S5P_OFFL_L2__NO2____20240409T051555_20240409T065725_33622_03_020600_20240410T213619.nc"] - - for file in file_list: - plot_global_var(file, 'O3____')