20241217-add merc

model_info_2d.py

@@ -17,7 +17,9 @@ class model_info_2d(object):
             nt          : int               = None,
             dt          : float             = None,
             var_list    : list              = None,
-            type        : str               = None
+            type        : str               = None,
+            globe       : ccrs.Globe        = None,
+            debug       : int               = 0,
     ) -> None:
         
         """
@@ -35,6 +37,8 @@ class model_info_2d(object):
             dt                  : 每个模式输出文件的时间间隔(小时)
             var_list            : 模式包含的变量列表
             type                : 模式的类型(只是一个标记)
+            globe               : 地球形状设定
+            debug               : 设置打印的信息
         更新记录:
             2022-08-20 22:08:27 Sola 编写源代码
             2022-08-20 22:08:33 Sola 添加注释
@@ -52,6 +56,9 @@ class model_info_2d(object):
             2023-03-18 16:22:17 Sola v5 增加支持获取加密网格的方法, 用于超采样清单
             2023-03-19 21:53:51 Sola v0.0.2 加入了默认的网格(经纬度网格), 以方便了解功能
             2023-04-29 18:54:06 Sola v0.0.3 加入了从WRF读取数据, 以及输出cartopy.crs的功能
+            2023-09-07 10:42:59 Sola v0.0.4 设定了默认的地球形状, 以修正默认投影与模式的偏差, 加入globe参数
+                感谢韩雨阳的帮助, 指出了两个差异的问题所在
+            2024-07-22 20:36:52 Sola v0.0.5 增加了判断坐标（坐标数组）是否在模式网格内的功能
         测试记录:
             2022-09-28 16:28:10 Sola v2 新的简化网格生成方法测试完成, 结果与旧版一致
             2022-09-28 18:27:59 Sola v2 测试了使用proj_LC投影的相关方法, 网格与WRF一致
@@ -66,6 +73,7 @@ class model_info_2d(object):
             self.dt = 1 if dt is None else dt # 时间间隔(小时)
             self.nt = 1 if nt is None else nt # 每个文件中包含多少时间点
             self.var_list = [] if var_list is None else var_list # 变量列表
+            self.globe = ccrs.Globe(ellipse="sphere", semimajor_axis=6370000, semiminor_axis=6370000) if globe is None else globe
             if lowerleft is None:
                 zero_lon, zero_lat = ccrs.PlateCarree().transform_point(
                     -self.dx*(self.nx-1)/2, -self.dy*(self.ny-1)/2, self.projection)
@@ -82,7 +90,8 @@ class model_info_2d(object):
                 ccrs.PlateCarree()
             ) # 计算投影下的坐标
         finally:
-            print(f"{self.__dict__}")
+            if debug > 0:
+                print(f"{self.__dict__}")
 
     def grid_id_float(self, original_x, original_y, original_proj=ccrs.PlateCarree()):
         """
@@ -214,16 +223,18 @@ class model_info_2d(object):
     def get_grid(self, type=None):
         """
         范围模式所有网格的经纬度坐标
+        type: None | "corner" | "edge" | "corner_2d"
         2023-03-14 10:05:43 Sola 更新边界宽度的功能及边缘网格的功能
                                     获取的边缘网格从左下角开始顺时针排序(左优先)
         2023-03-14 10:30:23 Sola 经过测试, 代码可以正常运行
         2023-03-18 15:40:20 Sola 删除边界宽度的功能(没有用了)
+        2024-08-02 18:01:48 Sola 添加生成边界经纬度的功能
         """
         # 获取网格信息, 下标从0开始
         ys, xs = np.meshgrid(range(self.ny), range(self.nx), indexing='ij')
         if type is None:
             xlon, xlat = self.grid_lonlats(xs, ys) # 从网格信息获取经纬度信息
-        elif type.lower() in ["corner", "c"]: # 四角的网格
+        elif type.lower() in ["corner", "c"]: # 四角的网格 (4, ny, nx)
             result = []
             result.append(self.grid_lonlats(xs - 0.5, ys - 0.5))
             result.append(self.grid_lonlats(xs - 0.5, ys + 0.5))
@@ -239,9 +250,12 @@ class model_info_2d(object):
             result.append(self.grid_lonlats(xs, ys - 0.5))
             xlon = np.array([x[0] for x in result])
             xlat = np.array([x[1] for x in result])
+        elif type.lower() in ["corner_2d", "c2d"]: # 四角网络 (ny + 1, nx + 1)
+            ys, xs = np.meshgrid(range(self.ny+1), range(self.nx+1), indexing='ij')
+            xlon, xlat = self.grid_lonlats(xs-0.5, ys-0.5)
         return xlon, xlat
 
-    def get_density_grid(self, density=10):
+    def get_density_grid(self, density=10, flat=False):
         """
         获取一个更密的网格, 原先的每个网格均匀返回多个点, 例如返回10*10=100个点
             可用于超采样, 以进行清单的分配等操作, 注意不要设置太大的密度, 否则
@@ -249,12 +263,16 @@ class model_info_2d(object):
         更新记录:
             2023-03-18 16:09:39 Sola 编写源代码
             2023-03-18 16:21:46 Sola 测试功能正常, 从网格到经纬度及反向都正常
+            2023-10-18 16:19:10 Sola 增加了将结果展开成2D的功能
         """
 
         sub_jj, sub_ii, jj, ii = np.meshgrid(range(density), range(density), 
             range(self.ny), range(self.nx), indexing='ij')
         fii = ii - 0.5 + (sub_ii + 0.5)/density
         fjj = jj - 0.5 + (sub_jj + 0.5)/density
+        if flat:
+            fii = np.transpose(fii, (2, 0, 3, 1)).reshape((self.ny*density, self.nx*density))
+            fjj = np.transpose(fjj, (2, 0, 3, 1)).reshape((self.ny*density, self.nx*density))
         xlonf, xlatf = self.grid_lonlats(fii, fjj)
 
         return xlonf, xlatf
@@ -271,8 +289,57 @@ class model_info_2d(object):
                     self.projection.truelat1,
                     self.projection.truelat2,
                 ],
+                globe = self.globe
             )
         return proj
+    
+    def get_extent(
+        self,
+        cx      : float = None,
+        cy      : float = None,
+        dx      : float = None,
+        dy      : float = None,
+        ratio   : float = 1
+    ) -> list:
+        """
+        用于获取指定数据范围的经纬度坐标
+        参数:
+            cx: 中心点x坐标
+            cy: 中心点y坐标
+            dx: 中心点周围x网格数
+            dy: 中心点周围y网格数
+        """
+        if cx is None:
+            cx, cy, dx, dy = self.nx/2, self.ny/2, self.nx/2, self.ny/2
+        XLON, XLAT = self.get_grid()
+        # ys, ye, xs, xe = np.floor(cy-dy), np.ceil(cy+dy), np.floor(cx-dx), np.ceil(cx+dx)
+        lon_start, _ = self.grid_lonlat(cx-dx*ratio, cy)
+        lon_end, _ = self.grid_lonlat(cx+dx*ratio, cy)
+        _, lat_start = self.grid_lonlat(cx, cy-dy*ratio)
+        _, lat_end = self.grid_lonlat(cx, cy+dy*ratio)
+        # if lon_start > lon_end:
+        #     lon_end += 360
+        # XLON, XLAT = XLON[cy-dy:cy+dy, cx-dx:cx+dx], XLAT[cy-dy:cy+dy, cx-dx:cx+dx]
+        # clon, clat = np.mean(XLON), np.mean(XLAT)
+        # dlon, dlat = (np.max(XLON) - np.min(XLON))/2*ratio, (np.max(XLAT) - np.min(XLAT))/2*ratio
+        # clon, clat = (lon_end + lon_start)/2, (lat_end + lat_start)/2
+        # dlon, dlat = (lon_end - lon_start)/2*ratio, (lat_end - lat_start)/2*ratio
+        # extent = [(clon-dlon+180)%360-180, (clon+dlon+180)%360-180, clat-dlat if clat-dlat>=-90 else -90, clat+dlat if clat+dlat<=90 else 90]
+        constrain_lon = lambda x: (x+180)%360-180
+        constrain_lat = lambda x: min(abs(x), 90) * (1 if x > 0 else -1)
+        extent = [constrain_lon(lon_start), constrain_lon(lon_end), constrain_lat(lat_start), constrain_lat(lat_end)]
+        return extent
+    
+    def is_in_domain(self, origin_x, origin_y, use_float=False):
+        """
+        用于判断坐标（经纬度）是否在模式网格范围内
+        """
+        if use_float:
+            ix, iy = self.grid_id_float(origin_x, origin_y)
+        else:
+            ix, iy = self.grid_id(origin_x, origin_y)
+        result = (0 <= ix) & (ix < self.nx) & (0 <= iy) & (iy < self.ny)
+        return result
 
 def flat_array(
         x : np.ndarray,

proj_info.py

@@ -1,9 +1,11 @@
 import numpy as np
 from math import radians, cos, tan, log10, sin, sqrt, atan2, atan, degrees
+import logging
 
 """
 更新记录:
     2022-09-23 11:59:01 Sola v1 编写源代码, 修正set_lc代码错误的问题
+    2024-12-16 09:54:40 Sola v2 增加变量检测的内容
 """
 
 EARTH_RADIUS_M = 6370000.
@@ -55,6 +57,34 @@ class proj_info(object):
         self.comp_ll    = comp_ll
         self.gauss_lat  = gauss_lat
 
+        if self.lat1:
+            if abs(self.lat1) > 90:
+                logging.error("Latitude of origin corner required as follows: -90N <= lat1 < = 90.N")
+        if self.lon1: # 限制经度范围
+            dummy_lon1 = (self.lon1 + 180) % 360 - 180
+            self.lon1 = dummy_lon1
+        if self.lon0: # 限制中央经线范围
+            dummy_lon0 = (self.lon0 + 180) % 360 - 180
+            self.lon0 = dummy_lon0
+        if self.dx:
+            if self.dx <= 0 and self.code != "PROJ_LATLON":
+                logging.error("Require grid spacing (dx) in meters be positive!")
+        if self.stdlon:
+            dummp_stdlon = (self.stdlon + 180) % 360 - 180
+            self.stdlon = dummp_stdlon
+        if self.truelat1:
+            if abs(self.truelat1) > 90:
+                logging.error("Set true latitude 1 for all projections!")
+        if not self.dy and self.dx: # 设置dy, 如果dy不存在, 则利用dx给定
+            self.dy = self.dx
+        if self.dx:
+            if self.code in ["PROJ_LC", "PROJ_PS", "PROJ_PS_WGS84", "PROJ_A:NERS_NAD83", "PROJ_MERC"]:
+                if self.truelat1 < 0: # 所在半球, 1为北半球, -1为南半球
+                    self.hemi = -1
+                else:
+                    self.hemi = 1
+                self.rebydx = self.re_m / self.dx # 地球半径除以网格距
+        
 
 class proj_LC(proj_info):
     """
@@ -259,6 +289,73 @@ class proj_LC(proj_info):
         cy = (cy1 + cy2) / 2
         return (ix - cx) * self.dx, (iy - cy) * self.dy
 
+class proj_MERC(proj_info):
+    """
+    参考WPS源码中的proj_LC改写, 因为WRF计算得到的网格与cartopy的不同
+    更新记录:
+        2022-09-22 22:07:51 Sola 编写源代码
+    """
+    def __init__(self, code='PROJ_MERC', truelat1=None, lat1=None,
+        lon1=None, knowni=None, knownj=None, stdlon=None, dx=None,
+        dy=None, nx=None, ny=None, re_m=EARTH_RADIUS_M) -> None:
+        """
+        初始化
+        必要参数:
+            code        投影编码
+            truelat1    标准纬线1
+            truelat2    标准纬线2
+            lat1        参考点纬度
+            lon1        参考点经度
+            stdlon      中央经线
+            dx          x方向网格距(m)
+            nx          x方向格点数
+            ny          y方向格点数
+        可选参数:
+            knowni      参考点x方向坐标, 默认为网格中心
+            knownj      参考点y方向坐标, 默认为网格中心
+            dy          y方向网格距(m), 默认与dx一致
+            re_m        地球半径, 默认为6370000
+        """
+        if truelat1 is None is None or lat1 is None or lon1 is None\
+            or nx is None or ny is None or dx is None:
+            print('[ERROR] cannot generate proj!')
+        if abs(lat1) > 90 or dx <= 0 or truelat1 > 90:
+            pass
+        dummy_lon1 = (lon1 + 180) % 360 - 180       # 限制经度范围
+        dummy_stdlon = (stdlon + 180) % 360 - 180   # 限制中央经线范围
+        if knowni is None and knownj is None:
+            knowni = (nx + 1) / 2
+            knownj = (ny + 1) / 2
+        if dy is None:          # 设置dy, 如果dy不存在, 则利用dx给定
+            dy = dx
+        if truelat1 < 0:        # 所在半球, 1为北半球, -1为南半球
+            hemi = -1
+        else:
+            hemi = 1
+        if abs(truelat2) > 90:  # 如果标准纬线2超过范围, 则用标准纬线1赋值
+            truelat2 = truelat1
+        super().__init__(code=code, lat1=lat1, lon1=dummy_lon1, dx=dx, dy=dy, 
+            stdlon=dummy_stdlon, truelat1=truelat1, truelat2=truelat2, hemi=hemi, 
+            knowni=knowni, knownj=knownj, re_m=re_m) # 初始化各变量
+        self.rebydx = re_m / dx # 地球半径除以网格距
+        self.set_lc()           # 计算其他变量
+        self.check_init()       # 确认是否所有变量都计算完毕
+    def set_merc(self):
+        clain = np.cos(np.deg2rad(self.truelat1))
+        self.dlon = self.dx / (self.re_m * clain)
+        # 计算原点到迟到的距离，并保存在 self.rsw 变量中
+        self.rsw = 0 if self.lat1 == 0 else np.log(np.tan(0.5*(np.deg2rad(self.lat1+90))))/self.dlon
+    def llij_merc(self, lat, lon):
+        deltalon = lon - self.lon1
+        deltalon = (deltalon + 180) % 360 - 180
+        i = self.knowni + deltalon / np.rad2deg(self.dlon)
+        j = self.knownj + np.log(np.tan(0.5*np.deg2rad(lat + 90)))/self.dlon - self.rsw
+        return i, j
+    def ijll_merc(self, i, j):
+        lat = np.rad2deg(2*np.arctan(np.exp(self.dlon*(self.rsw + j - self.knownj)))) - 90
+        lon = (i - self.knowni) * np.rad2deg(self.dlon) + self.lon1
+        lon = (lon + 180) % 360 - 180
+        return lon, lat
 
 if __name__ == '__main__':
     proj = proj_LC(truelat1=45, truelat2=15, lat1=30, lon1=108, stdlon=108, dx=3000, dy=3000, nx=2025, ny=2025)