20241218 update proj merc

proj_info.py

@@ -6,8 +6,11 @@ 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):
@@ -19,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)经度, 西南角, 度
@@ -56,7 +59,12 @@ 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")
@@ -84,7 +92,30 @@ class proj_info(object):
                 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):
     """
@@ -116,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, 
-            stdlon=dummy_stdlon, truelat1=truelat1, truelat2=truelat2, hemi=hemi, 
-            knowni=knowni, knownj=knownj, re_m=re_m) # 初始化各变量
-        self.rebydx = re_m / dx # 地球半径除以网格距
+        super().__init__(code=code, lat1=lat1, lon1=lon1, dx=dx, dy=dy, 
+            stdlon=stdlon, truelat1=truelat1, truelat2=truelat2, 
+            knowni=knowni, knownj=knownj, re_m=re_m, nx=nx, ny=ny) # 初始化各变量
         self.set_lc()           # 计算其他变量
         self.check_init()       # 确认是否所有变量都计算完毕
 
@@ -167,23 +184,7 @@ class proj_LC(proj_info):
         else:
             self.cone = sin(radians(abs(self.truelat1)))
 
-    def llij_lc(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):
+    def llij(self, lon, lat):
         """通过经纬度序列计算坐标"""
         if not self.init:
             print('[ERROR] proj cannot use!')
@@ -199,32 +200,7 @@ class proj_LC(proj_info):
         j = self.hemi * j
         return i, j
 
-    def ijll_lc(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):
+    def ijll(self, i, j):
         """
         通过坐标计算经纬度
         2022-09-28 18:22:26 Sola 修正计算chi时, xx, yy未进行筛选的问题
@@ -235,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
-        yy = self.polej - jnew
+        xx = np.array(inew - self.polei)
+        yy = np.array(self.polej - jnew)
         r2 = xx**2 + yy**2
         r = np.sqrt(r2) / self.rebydx
         lon = np.zeros(r2.shape)
@@ -261,39 +237,20 @@ 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)
-        cx, cy1 = self.llij_lc(self.stdlon, self.truelat1)
-        cx, cy2 = self.llij_lc(self.stdlon, self.truelat2)
+        ix, iy = self.llij(lon, lat)
+        cx, cy1 = self.llij(self.stdlon, self.truelat1)
+        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_LC改写, 因为WRF计算得到的网格与cartopy的不同
+    参考WPS源码中的proj_MERC改写, 因为WRF计算得到的网格与cartopy的不同
     更新记录:
-        2022-09-22 22:07:51 Sola 编写源代码
+        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,
@@ -319,54 +276,72 @@ class proj_MERC(proj_info):
         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()       # 确认是否所有变量都计算完毕
+        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 变量中
+        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):
+
+    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_merc(self, i, j):
-        lat = np.rad2deg(2*np.arctan(np.exp(self.dlon*(self.rsw + j - self.knownj)))) - 90
+    
+    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
 
 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)
-    x = np.arange(1, 2000, 10)
-    y = np.arange(1, 2000, 10)
-    print(proj.ijlls_lc(x, y))
-    print(proj.ijll_lc(x[0], y[0]), proj.ijll_lc(x[-1], y[-1]))
-    x, y = proj.llij_lc(0, 90)
-    lon, lat = proj.ijll_lc(x, y)
-    x, y = proj.llij_lc(lon, lat)
-    lon1, lat1 = proj.ijll_lc(x, y)
-    print(lon, lat, x, y, lon1, lat1)
+    """
+    主要用于测试投影的基本功能，坐标前后转换是否正常，具体测试项如下：
+    1. 构建投影
+    2. 坐标点位转换以及转换结果是否一致
+    3. 坐标序列转换以及转换结果是否一致
+    """
+    logging.info("test proj_LC")
+    try:
+        proj = proj_LC(truelat1=45, truelat2=15, 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((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.")