Merge branch 'master' of gitea:/Sola/model_info_2d

Conflicts:
	model_info_2d.py

__init__.py

@@ -1,5 +1,5 @@
 import cartopy.crs as ccrs
 import numpy as np
-from .model_info_2d import model_info_2d, from_wrf
+from .model_info_2d import model_info_2d, from_wrf, from_ctl
 
-__all__ = ["model_info_2d", "ccrs", "np", "from_wrf"]
+__all__ = ["model_info_2d", "ccrs", "np", "from_wrf", "from_ctl"]

model_info_2d.py

@@ -1,5 +1,7 @@
 import numpy as np
 import cartopy.crs as ccrs
+from typing import Union
+from .proj_info import proj_LC, proj_MERC
 
 class model_info_2d(object):
     """
@@ -20,6 +22,9 @@ class model_info_2d(object):
             type        : str               = None,
             globe       : ccrs.Globe        = None,
             debug       : int               = 0,
+            center      : list              = None,
+            rotate_deg  : Union[int, float] = 0,
+            rotate_poi  : list              = None,
     ) -> None:
         
         """
@@ -60,6 +65,15 @@ class model_info_2d(object):
                 感谢韩雨阳的帮助, 指出了两个差异的问题所在
             2023-12-28 15:42:11 Sola v0.0.5 增加了加密网格的功能
             2023-12-28 15:54:53 Sola v0.0.6 增加了获取绘图范围的功能, 并使其接受浮点数输入
+            2024-07-22 20:36:52 Sola v0.0.7 增加了判断坐标（坐标数组）是否在模式网格内的功能
+            2024-12-18 10:11:55 Sola v0.0.8 增加了与墨卡托投影相关的计算内容
+            2025-04-06 16:39:26 Sola v0.0.9 增加提供网格中心坐标计算网格的功能（优先级低于左下角坐标）
+            2025-04-06 16:45:22 Sola v0.0.10 增加坐标旋转功能
+                修改的关键在于：
+                1. 在将经纬度转化为网格的时候, 围绕中心对网格进行偏移旋转, 需要增加一步后处理
+                2. 在将网格转化为经纬度的时候, 需要先将输入的网格ID旋转回去, 再计算其经纬度
+                设计的网格旋转函数需要保证旋转前后中心位置不变，各网格相对位置不变即可
+                注意, 这里输入的左下角坐标与通过中心计算的左下角坐标均为旋转前的
         测试记录:
             2022-09-28 16:28:10 Sola v2 新的简化网格生成方法测试完成, 结果与旧版一致
             2022-09-28 18:27:59 Sola v2 测试了使用proj_LC投影的相关方法, 网格与WRF一致
@@ -76,20 +90,36 @@ class model_info_2d(object):
             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:
+                if not center is None:
+                    center_x, center_y = self.projection.transform_point(center[0], center[1], ccrs.PlateCarree())
+                else:
+                    center_x, center_y = 0, 0
                 zero_lon, zero_lat = ccrs.PlateCarree().transform_point(
-                    -self.dx*(self.nx-1)/2, -self.dy*(self.ny-1)/2, self.projection)
+                    center_x-self.dx*(self.nx-1)/2, center_y-self.dy*(self.ny-1)/2, self.projection)
-                self.lowerleft = [zero_lon, zero_lat]
+                self.lowerleft = [zero_lon, zero_lat] # 旋转前的左下角坐标
             else:
                 if len(lowerleft) == 2:
-                    self.lowerleft = lowerleft # 左下角坐标(经纬度)
+                    self.lowerleft = lowerleft # 旋转前的左下角坐标(经纬度)
                 else:
+                    # 这是考虑输入的左下角坐标不是经纬度, 而是某个投影系下的坐标位置, 所以先将其转化为经纬度
                     zero_lon, zero_lat = ccrs.PlateCarree().transform_point(\
                         lowerleft[0], lowerleft[1], lowerleft[2])
-                    self.lowerleft = [zero_lon, zero_lat]
+                    self.lowerleft = [zero_lon, zero_lat] # 旋转前的左下角经纬度
             self.lowerleft_projxy = self.projection.transform_point(
                 self.lowerleft[0], self.lowerleft[1],
                 ccrs.PlateCarree()
-            ) # 计算投影下的坐标
+            ) # 计算投影下的xy坐标
+            self.rotate = 0 if rotate_deg is None else np.deg2rad(rotate_deg) # 计算旋转的弧度(输入是角度)
+            if rotate_poi is None:
+                # 如果没有给定围绕旋转的点位, 则围绕网格中心进行旋转, 注意这里是 (x, y), 而不是 (ix, iy)
+                # 注意需要考虑如果指定的网格中心和投影中心不一致的情况
+                if not center is None:
+                    self.rotate_poi_x, self.rotate_poi_y = center_x, center_y
+                else:
+                    self.rotate_poi_x, self.rotate_poi_y = self.lowerleft_projxy[0] + (self.nx - 1)*self.dx, self.lowerleft_projxy[1] + (self.ny - 1)*self.dy
+            else:
+                # 如果
+                self.rotate_poi_x, self.rotate_poi_y = self.projection.transform_point(*rotate_poi, ccrs.PlateCarree())
         finally:
             if debug > 0:
                 print(f"{self.__dict__}")
@@ -100,6 +130,7 @@ class model_info_2d(object):
         2022-09-28 11:05:09 Sola 更新为识别传入的对象类型, 判断是否可迭代
         2022-09-28 15:21:07 Sola 增加对proj是否包含相应方法的识别
         2022-09-28 18:25:24 Sola 修正正常情况下未输出ix, iy的bug
+        2025-04-06 20:33:29 Sola 加入坐标旋转的判断
         """
         # 如果是可迭代对象, 则丢给对应的功能处理
         if hasattr(original_x, '__iter__'):
@@ -116,10 +147,10 @@ class model_info_2d(object):
                 # 调用proj的方法计算经纬度
                 ix, iy = self.projection.grid_id_float(lon, lat)
             else: # 如果投影方法本身不具备计算网格ID的方法, 那就手动计算网格
-                x, y = self.projection.transform_point(
+                x, y = self.projection.transform_point(original_x, original_y, original_proj)
-                    original_x, original_y, original_proj)
                 ix = (x - self.lowerleft_projxy[0])/self.dx
                 iy = (y - self.lowerleft_projxy[1])/self.dy
+            ix, iy = self.rotate_grid_revise(ix, iy)
         return ix, iy
 
     def grid_id(self, original_x, original_y, original_proj=ccrs.PlateCarree()):
@@ -145,6 +176,7 @@ class model_info_2d(object):
         2022-09-28 16:40:27 Sola 增加将输入数组转化为numpy数组的功能, 防止传入列表
         2022-10-19 18:52:25 Sola 修正了除错距离的bug
         2023-03-18 15:39:06 Sola 在计算前, 先将数组展开到1维, 返回时折叠
+        2025-04-06 20:33:12 Sola 加入坐标旋转的判断
         注意事项:
             当前存在一个bug, 输入的投影必须是cartopy的投影, 否则无法计算经纬度,
             但是是否有必要在自己写的proj中加入该功能? 需要考虑
@@ -170,6 +202,7 @@ class model_info_2d(object):
             ix_array = ((ix_array - self.lowerleft_projxy[0])/ self.dx).T
             iy_array = ((iy_array - self.lowerleft_projxy[1])/ self.dy).T
         ix_array, iy_array = fold_array(ix_array, iy_array, shape)
+        ix_array, iy_array = self.rotate_grid_revise(ix_array, iy_array)
         return ix_array, iy_array
     
     def grid_ids(self, original_x_array, original_y_array,
@@ -189,13 +222,16 @@ class model_info_2d(object):
         通过网格id获取经纬度坐标
         2022-09-28 16:03:27 Sola 增加判断传入的是数值还是数组的功能
         2022-09-28 16:05:07 Sola 增加判断proj是否有计算网格的功能
+        2025-04-06 20:32:55 Sola 加入坐标旋转的判断
         """
         if hasattr(ix, '__iter__'): # 如果传入的是可迭代对象, 则调用相应功能
             lon, lat = self.grid_lonlats(ix, iy)
         else: # 如果不是, 则由本函数继续运算
+            ix, iy = self.rotate_grid(ix, iy)
             if hasattr(self.projection, 'grid_lonlat'): # 如果投影本身可以计算
                 lon, lat = self.projection.grid_lonlat(ix, iy) # 计算网格对应经纬度
             else: # 如果投影不能根据网格ID计算经纬度, 则手动计算
+                # 这里则是根据网格计算了在给定投影下的坐标XY，然后将其转化为经纬度
                 x = self.lowerleft_projxy[0] + ix * self.dx
                 y = self.lowerleft_projxy[1] + iy * self.dy
                 lon, lat = ccrs.PlateCarree().transform_point(x, y, self.projection)
@@ -208,8 +244,10 @@ class model_info_2d(object):
         2022-09-28 16:08:38 Sola 简化原本的网格计算, 使用转置的方式代替判断返回数组长度
         2022-09-28 16:40:27 Sola 增加将输入数组转化为numpy数组的功能, 防止传入列表
         2023-03-18 15:39:06 Sola 在计算前, 先将数组展开到1维, 返回时折叠
+        2025-04-06 20:33:56 Sola 加入坐标旋转的判断
         """
         ix_array, iy_array, shape = flat_array(np.array(ix_array), np.array(iy_array))
+        ix_array, iy_array = self.rotate_grid(ix_array, iy_array)
         if hasattr(self.projection, 'grid_lonlats'):
             lon_array, lat_array = self.projection.grid_lonlats(ix_array, iy_array)
         else:
@@ -224,16 +262,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))
@@ -249,6 +289,9 @@ 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, flat=False):
@@ -278,7 +321,6 @@ class model_info_2d(object):
         """
         获取用于绘图的地图投影, 目前只支持兰伯特投影
         """
-        from proj_info import proj_LC
         if type(self.projection) is proj_LC:
             proj = ccrs.LambertConformal(
                 central_longitude   = self.projection.stdlon,
@@ -288,14 +330,23 @@ class model_info_2d(object):
                 ],
                 globe = self.globe
             )
+        elif type(self.projection) is proj_MERC:
+            proj = ccrs.Mercator(
+                central_longitude=self.projection.stdlon,
+                globe = self.globe
+            )
+        elif self.projection.__class__.__base__ is ccrs.Projection:
+            proj = self.projection
+        else:
+            proj = ccrs.PlateCarree(globe = self.globe)
         return proj
     
     def get_extent(
         self,
-        cx      : float,
+        cx      : float = None,
-        cy      : float,
+        cy      : float = None,
-        dx      : float,
+        dx      : float = None,
-        dy      : float,
+        dy      : float = None,
         ratio   : float = 1
     ) -> list:
         """
@@ -306,7 +357,9 @@ class model_info_2d(object):
             dx: 中心点周围x网格数
             dy: 中心点周围y网格数
         """
-        XLON, XLAT = self.get_grid()
+        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)
@@ -324,6 +377,56 @@ class model_info_2d(object):
         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):
+        """
+        用于判断坐标（经纬度）是否在模式网格范围内
+        Update:
+            2025-05-05 00:13:01 Sola 修正使用浮点数计算时的问题
+        """
+        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 - 1) & (0 <= iy) & (iy <= self.ny - 1)
+        return result
+
+    def rotate_xy(self, x, y, rotate_rad=None):
+        rotate_rad = self.rotate if rotate_rad is None else rotate_rad
+        x_new, y_new = rotate_xy(x, y, self.rotate_poi_x, self.rotate_poi_y, rotate_rad)
+        return x_new, y_new
+    
+    def rotate_xy_revise(self, x, y, rotate_rad=None):
+        rotate_rad = self.rotate if rotate_rad is None else rotate_rad
+        x_new, y_new = rotate_xy(x, y, self.rotate_poi_x, self.rotate_poi_y, -rotate_rad)
+        return x_new, y_new
+    
+    def rotate_grid(self, ix, iy, rotate_rad=None):
+        rotate_rad = self.rotate if rotate_rad is None else rotate_rad
+        if np.sum(np.abs(rotate_rad % (np.pi*2)) > 1e-8):
+            x, y = self.lowerleft_projxy[0] + ix*self.dx, self.lowerleft_projxy[1] + iy*self.dy
+            x_new, y_new = self.rotate_xy(x, y, rotate_rad)
+            ix_new, iy_new = (x_new - self.lowerleft_projxy[0])/self.dx, (y_new - self.lowerleft_projxy[1])/self.dy
+        else:
+            ix_new, iy_new = ix, iy
+        return ix_new, iy_new
+    
+    def rotate_grid_revise(self, ix, iy, rotate_rad=None):
+        rotate_rad = self.rotate if rotate_rad is None else rotate_rad
+        if np.sum(np.abs(rotate_rad % (np.pi*2)) > 1e-8):
+            x, y = self.lowerleft_projxy[0] + ix*self.dx, self.lowerleft_projxy[1] + iy*self.dy
+            x_new, y_new = self.rotate_xy_revise(x, y, rotate_rad)
+            ix_new, iy_new = (x_new - self.lowerleft_projxy[0])/self.dx, (y_new - self.lowerleft_projxy[1])/self.dy
+        else:
+            ix_new, iy_new = ix, iy
+        return ix_new, iy_new
+
+
+def rotate_xy(xx, yy, cx, cy, rad):
+    xx_offset = (xx - cx)*np.cos(rad) - (yy - cy)*np.sin(rad)
+    yy_offset = (xx - cx)*np.sin(rad) + (yy - cy)*np.cos(rad)
+    xx_new, yy_new = cx + xx_offset, cy + yy_offset
+    return xx_new, yy_new
 
 def flat_array(
         x : np.ndarray,
@@ -371,20 +474,42 @@ def from_wrf(file: str) -> model_info_2d:
     
     # import need library
     import netCDF4 as nc
-    from proj_info import proj_LC
     # open dataset
     with nc.Dataset(file) as nf:
+        dx, dy = nf.DX, nf.DY
+        nx, ny = nf.dimensions["west_east"].size, nf.dimensions["south_north"].size
+        truelat1 = nf.TRUELAT1
+        stdlon = nf.STAND_LON
+        lat1, lon1 = nf.CEN_LAT, nf.CEN_LON
         if nf.MAP_PROJ == 1: # Lambert proj
-            dx, dy = nf.DX, nf.DY
+            truelat2 = nf.TRUELAT2
-            nx, ny = nf.dimensions["west_east"].size, nf.dimensions["south_north"].size
-            truelat1, truelat2 = nf.TRUELAT1, nf.TRUELAT2
-            stdlon = nf.STAND_LON
-            lat1, lon1 = nf.CEN_LAT, nf.CEN_LON
             # make projection
             proj = proj_LC(dx=dx, dy=dy, truelat1=truelat1, truelat2=truelat2,
                         lat1=lat1, lon1=lon1, stdlon=stdlon, nx=nx, ny=ny)
-            # make model_info
+        elif nf.MAP_PROJ == 3: # Mercator proj
-            model = model_info_2d(proj=proj, nx=nx, ny=ny, dx=dx, dy=dy,
+            proj = proj_MERC(dx=dx, dy=dy, truelat1=truelat1, lat1=lat1,
-                                lowerleft=proj.grid_lonlat(0, 0))
+                             lon1=lon1, stdlon=stdlon, nx=nx, ny=ny)
-
+        elif nf.MAP_PROJ == 6: # lon-lat proj
+            proj = ccrs.PlateCarree
+            dx /= 111177.473
+            dy /= 111177.473
+            
+        # make model_info
+        model = model_info_2d(proj=proj, nx=nx, ny=ny, dx=dx, dy=dy,
+                            lowerleft=proj.grid_lonlat(0, 0))
+    return model
+
+def from_ctl(file: str) -> model_info_2d:
+    with open(file, "r") as f:
+        lines = f.readlines()
+    for line in lines:
+        if "PDEF" in line.upper():
+            if "LCC" in line.upper():
+                _, nx, ny, _, lat1, lon1, knowi, knowj, truelat1, truelat2, stdlon, dx, dy = line.split()
+                nx, ny = int(nx), int(ny)
+                lat1, lon1, knowi, knowj, truelat1, truelat2, stdlon, dx, dy =\
+                    float(lat1), float(lon1), float(knowi), float(knowj), float(truelat1), float(truelat2), float(stdlon), float(dx), float(dy)
+                proj = proj_LC(dx=dx, dy=dy, truelat1=truelat1, truelat2=truelat2, lat1=lat1, lon1=lon1,
+                    knowni=knowi, knownj=knowj, stdlon=stdlon, nx=nx, ny=ny)
+    model = model_info_2d(proj=proj, nx=nx, ny=ny, dx=dx, dy=dy, lowerleft=proj.grid_lonlat(0, 0))
     return model

proj_info.py

@@ -1,11 +1,16 @@
 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 增加变量检测的内容
+    2024-12-17 15:39:45 Sola v3 增加墨卡托投影
 """
 
+# logging.basicConfig(format='[%(asctime)s][%(levelname)s]: %(message)s',
+#                     level=logging.DEBUG, datefmt='%Y-%m-%dT%H:%M:%S %Z')
 EARTH_RADIUS_M = 6370000.
 
 class proj_info(object):
@@ -17,7 +22,7 @@ class proj_info(object):
         nxmax=None, hemi=None, cone=None, polei=None, polej=None, 
         rsw=None, knowni=None, knownj=None, re_m=EARTH_RADIUS_M, 
         init=False, wrap=False, rho0=None, nc=None, bigc=None, comp_ll=False, 
-        gauss_lat=None) -> None:
+        gauss_lat=None, nx=None, ny=None) -> None:
         self.code       = code
         self.lat1       = lat1      # SW latitude (1,1) in degrees (-90->90N) 格点(1, 1)纬度, 西南角, 度
         self.lon1       = lon1      # SW longitude (1,1) in degrees (-180->180E) 格点(1, 1)经度, 西南角, 度
@@ -54,7 +59,63 @@ class proj_info(object):
         self.bigc       = bigc
         self.comp_ll    = comp_ll
         self.gauss_lat  = gauss_lat
+        self.nx         = nx
+        self.ny         = ny
 
+        if self.knowni is None and self.knownj is None:
+            self.knowni = (self.nx + 1) / 2
+            self.knownj = (self.ny + 1) / 2
+        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 # 地球半径除以网格距
+    
+    def grid_id_float(self, lon, lat):
+        """返回以0为开始的下标"""
+        i,j = self.llij(lon, lat)
+        return i - 1, j - 1
+
+    def grid_ids_float(self, lon_array, lat_array):
+        """返回以0为开始的下标数组"""
+        i_array, j_array = self.llij(lon_array, lat_array)
+        return i_array - 1, j_array - 1
+
+    def grid_lonlat(self, ix, iy):
+        """返回对应网格(以0为下标开始)的经纬度"""
+        lon, lat = self.ijll(ix + 1, iy + 1)
+        return lon, lat
+
+    def grid_lonlats(self, ix_array, iy_array):
+        """返回对应网格数组(以0为下标开始)的经纬度数组"""
+        lon_array, lat_array = self.ijll(ix_array + 1, iy_array + 1)
+        return lon_array, lat_array
+    
+    def transform_point(self, lon, lat, proj_useless=None):
+        """返回对应经纬度坐标的网格坐标(m)"""
+        pass
 
 class proj_LC(proj_info):
     """
@@ -86,25 +147,11 @@ class proj_LC(proj_info):
         if truelat1 is None or truelat2 is None or lat1 is None or lon1 is None\
             or nx is None or ny is None or stdlon 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, 
+        super().__init__(code=code, lat1=lat1, lon1=lon1, dx=dx, dy=dy, 
-            stdlon=dummy_stdlon, truelat1=truelat1, truelat2=truelat2, hemi=hemi, 
+            stdlon=stdlon, truelat1=truelat1, truelat2=truelat2, 
-            knowni=knowni, knownj=knownj, re_m=re_m) # 初始化各变量
+            knowni=knowni, knownj=knownj, re_m=re_m, nx=nx, ny=ny) # 初始化各变量
-        self.rebydx = re_m / dx # 地球半径除以网格距
         self.set_lc()           # 计算其他变量
         self.check_init()       # 确认是否所有变量都计算完毕
 
@@ -137,23 +184,7 @@ class proj_LC(proj_info):
         else:
             self.cone = sin(radians(abs(self.truelat1)))
 
-    def llij_lc(self, lon, lat):
+    def llij(self, lon, lat):
-        """通过经纬度计算坐标"""
-        if not self.init:
-            print('[ERROR] proj cannot use!')
-        deltalon = lon - self.stdlon            # 计算经度与中央经线差
-        deltalon = (deltalon + 180) % 360 - 180 # 限定范围 -180~180
-        ctl1r = cos(radians(self.truelat1))     # 剩下就看不懂了
-        rm = self.rebydx * ctl1r / self.cone * (tan(radians(90*self.hemi \
-            - lat)/2)/tan(radians(90*self.hemi - self.truelat1)/2))**self.cone
-        arg = self.cone * radians(deltalon)
-        i = self.polei + self.hemi * rm * sin(arg)
-        j = self.polej - rm * cos(arg)
-        i = self.hemi * i
-        j = self.hemi * j
-        return i, j
-
-    def llijs_lc(self, lon, lat):
         """通过经纬度序列计算坐标"""
         if not self.init:
             print('[ERROR] proj cannot use!')
@@ -169,32 +200,7 @@ class proj_LC(proj_info):
         j = self.hemi * j
         return i, j
 
-    def ijll_lc(self, i, j):
+    def ijll(self, i, j):
-        """通过坐标计算经纬度"""
-        if not self.init:
-            print('[ERROR] proj cannot use!')
-        chi1 = radians(90 - self.hemi * self.truelat1)
-        chi2 = radians(90 - self.hemi * self.truelat2)
-        inew = self.hemi * i
-        jnew = self.hemi * j
-        xx = inew - self.polei
-        yy = self.polej - jnew
-        r2 = xx**2 + yy**2
-        r = sqrt(r2) / self.rebydx
-        if r2 == 0.:
-            lat = self.hemi * 90
-            lon = self.stdlon
-        else:
-            lon = self.stdlon + degrees(atan2(self.hemi*xx, yy))/self.cone
-            lon = (lon + 180) % 360 - 180
-            if chi1 == chi2:
-                chi = 2 * atan((r/tan(chi1))**(1/self.cone) * tan(chi1*0.5))
-            else:
-                chi = 2 * atan((r*self.cone/sin(chi1))**(1/self.cone) * tan(chi1*0.5))
-            lat = (90 - degrees(chi)) * self.hemi
-        return lon, lat
-
-    def ijlls_lc(self, i, j):
         """
         通过坐标计算经纬度
         2022-09-28 18:22:26 Sola 修正计算chi时, xx, yy未进行筛选的问题
@@ -205,8 +211,8 @@ class proj_LC(proj_info):
         chi2 = np.radians(90 - self.hemi * self.truelat2)
         inew = self.hemi * i
         jnew = self.hemi * j
-        xx = inew - self.polei
+        xx = np.array(inew - self.polei)
-        yy = self.polej - jnew
+        yy = np.array(self.polej - jnew)
         r2 = xx**2 + yy**2
         r = np.sqrt(r2) / self.rebydx
         lon = np.zeros(r2.shape)
@@ -231,45 +237,116 @@ class proj_LC(proj_info):
         else:
             self.init = True
 
-    def grid_id_float(self, lon, lat):
-        """返回以0为开始的下标"""
-        i,j = self.llij_lc(lon, lat)
-        return i - 1, j - 1
-
-    def grid_ids_float(self, lon_array, lat_array):
-        """返回以0为开始的下标数组"""
-        i_array, j_array = self.llijs_lc(lon_array, lat_array)
-        return i_array - 1, j_array - 1
-
-    def grid_lonlat(self, ix, iy):
-        """返回对应网格(以0为下标开始)的经纬度"""
-        lon, lat = self.ijll_lc(ix + 1, iy + 1)
-        return lon, lat
-
-    def grid_lonlats(self, ix_array, iy_array):
-        """返回对应网格数组(以0为下标开始)的经纬度数组"""
-        lon_array, lat_array = self.ijlls_lc(ix_array + 1, iy_array + 1)
-        return lon_array, lat_array
-
     def transform_point(self, lon, lat, proj_useless=None):
         """返回对应经纬度坐标的网格坐标(m)"""
-        ix, iy = self.llij_lc(lon, lat)
+        ix, iy = self.llij(lon, lat)
-        cx, cy1 = self.llij_lc(self.stdlon, self.truelat1)
+        cx, cy1 = self.llij(self.stdlon, self.truelat1)
-        cx, cy2 = self.llij_lc(self.stdlon, self.truelat2)
+        cx, cy2 = self.llij(self.stdlon, self.truelat2)
         cy = (cy1 + cy2) / 2
         return (ix - cx) * self.dx, (iy - cy) * self.dy
 
 
+class proj_MERC(proj_info):
+    """
+    参考WPS源码中的proj_MERC改写, 因为WRF计算得到的网格与cartopy的不同
+    更新记录:
+        2024-12-17 15:12: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
+            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!')
+        super().__init__(code=code, lat1=lat1, lon1=lon1, dx=dx, dy=dy, nx=nx, ny=ny,
+            stdlon=stdlon, truelat1=truelat1, knowni=knowni, knownj=knownj, re_m=re_m) # 初始化各变量
+        self.set_merc()           # 计算其他变量
+
+    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(self, lon, lat):
+        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(self, i, j):
+        lat = 2*np.rad2deg(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
+    
+    def transform_point(self, lon, lat, proj_useless=None):
+        """返回对应经纬度坐标的网格坐标(m)"""
+        ix, iy = self.llij(lon, lat)
+        cx, cy = self.llij(self.stdlon, self.truelat1)
+        return (ix - cx) * self.dx, (iy - cy) * self.dy
+
 if __name__ == '__main__':
-    proj = proj_LC(truelat1=45, truelat2=15, lat1=30, lon1=108, stdlon=108, dx=3000, dy=3000, nx=2025, ny=2025)
+    """
-    proj.llij_lc(108, 30)
+    主要用于测试投影的基本功能，坐标前后转换是否正常，具体测试项如下：
-    proj.ijll_lc(1013, 1013)
+    1. 构建投影
-    x = np.arange(1, 2000, 10)
+    2. 坐标点位转换以及转换结果是否一致
-    y = np.arange(1, 2000, 10)
+    3. 坐标序列转换以及转换结果是否一致
-    print(proj.ijlls_lc(x, y))
+    """
-    print(proj.ijll_lc(x[0], y[0]), proj.ijll_lc(x[-1], y[-1]))
+    logging.info("test proj_LC")
-    x, y = proj.llij_lc(0, 90)
+    try:
-    lon, lat = proj.ijll_lc(x, y)
+        proj = proj_LC(truelat1=45, truelat2=15, lat1=30, lon1=108, stdlon=108, dx=3000, dy=3000, nx=2025, ny=2025)
-    x, y = proj.llij_lc(lon, lat)
+        x0, y0 = 0, 0
-    lon1, lat1 = proj.ijll_lc(x, y)
+        lon, lat = proj.grid_lonlat(x0, y0)
-    print(lon, lat, x, y, lon1, lat1)
+        x1, y1 = proj.grid_id_float(lon, lat)
+        if not ((x0 - x1)**2 + (y0 - y1)**2)**0.5 < 1e-5:
+            logging.error((x0, y0), (x1, y1))
+            raise ValueError("point convert error!")
+        x0_array, y0_array = np.arange(1, 2000, 10), np.arange(1, 2000, 10)
+        lon_array, lat_array = proj.grid_lonlat(x0_array, y0_array)
+        x1_array, y1_array = proj.grid_id_float(lon_array, lat_array)
+        if not np.max(((x0_array - x1_array)**2 + (y0_array - y1_array)**2)**0.5) < 1e-5:
+            raise ValueError("array convert error!")
+    except Exception as e:
+        logging.error("proj_LC not pass!")
+        logging.error(e)
+    else:
+        logging.info("proj_LC pass.")
+    
+    logging.info("test proj_MERC")
+    try:
+        proj = proj_MERC(truelat1=30, lat1=30, lon1=108, stdlon=108, dx=3000, dy=3000, nx=2025, ny=2025)
+        x0, y0 = 0, 0
+        lon, lat = proj.grid_lonlat(x0, y0)
+        x1, y1 = proj.grid_id_float(lon, lat)
+        if not ((x0 - x1)**2 + (y0 - y1)**2)**0.5 < 1e-5:
+            logging.error(f"{(x0, y0)}, {(x1, y1)}")
+            raise ValueError("point convert error!")
+        x0_array, y0_array = np.arange(1, 2000, 10), np.arange(1, 2000, 10)
+        lon_array, lat_array = proj.grid_lonlat(x0_array, y0_array)
+        x1_array, y1_array = proj.grid_id_float(lon_array, lat_array)
+        if not np.max(((x0_array - x1_array)**2 + (y0_array - y1_array)**2)**0.5) < 1e-5:
+            raise ValueError("array convert error!")
+    except Exception as e:
+        logging.error("proj_MERC not pass!")
+        logging.error(e)
+    else:
+        logging.info("proj_MERC pass.")