v0.0.13 改进了旋转的输入, 增加了读取 proj4 string 的功能

model_info_2d.py

@@ -23,8 +23,9 @@ class model_info_2d(object):
             globe       : ccrs.Globe        = None,
             debug       : int               = 0,
             center      : list              = None,
-            rotate_deg  : Union[int, float] = 0,
+            rotate_deg  : Union[int, float] = None,
             rotate_poi  : list              = None,
+            wind_dir_rad: float             = None,
     ) -> None:
         
         """
@@ -37,11 +38,11 @@ class model_info_2d(object):
             dx      : x方向网格距离(在目标网格投影下, 例如兰伯特是米, 等经纬是度)
             dy      : y方向网格距离
         可选参数:
-            lowerleft_lonlat    : 左下角坐标(经纬度)
+            lowerleft_lonlat    : 左下角 (0, 0) 位置坐标 (经纬度)
             nt                  : 每个模式输出文件的时间段个数
-            dt                  : 每个模式输出文件的时间间隔(小时)
+            dt                  : 每个模式输出文件的时间间隔 (小时)
             var_list            : 模式包含的变量列表
-            type                : 模式的类型(只是一个标记)
+            type                : 模式的类型 (只是一个标记)
             globe               : 地球形状设定
             debug               : 设置打印的信息
         更新记录:
@@ -74,6 +75,9 @@ class model_info_2d(object):
                 2. 在将网格转化为经纬度的时候, 需要先将输入的网格ID旋转回去, 再计算其经纬度
                 设计的网格旋转函数需要保证旋转前后中心位置不变，各网格相对位置不变即可
                 注意, 这里输入的左下角坐标与通过中心计算的左下角坐标均为旋转前的
+            2025-07-14 15:42:22 Sola v0.0.11 增加select方法, 用于选取某个经纬度范围的数据
+            2025-07-14 23:24:51 Sola v0.0.12 改进了select方法, 并增加了判断是否在某个extent内的功能
+            2025-08-03 21:51:56 Sola v0.0.13 改进了旋转的输入, 增加了读取 proj4 string 的功能
         测试记录:
             2022-09-28 16:28:10 Sola v2 新的简化网格生成方法测试完成, 结果与旧版一致
             2022-09-28 18:27:59 Sola v2 测试了使用proj_LC投影的相关方法, 网格与WRF一致
@@ -109,7 +113,13 @@ class model_info_2d(object):
                 self.lowerleft[0], self.lowerleft[1],
                 ccrs.PlateCarree()
             ) # 计算投影下的xy坐标
-            self.rotate = 0 if rotate_deg is None else np.deg2rad(rotate_deg) # 计算旋转的弧度(输入是角度)
+            if rotate_deg is None:
+                if wind_dir_rad is None:
+                    self.rotate = 0
+                else:
+                    self.rotate = -wind_dir_rad
+            else:
+                self.rotate = np.deg2rad(rotate_deg) # 计算旋转的弧度(输入是角度)
             if rotate_poi is None:
                 # 如果没有给定围绕旋转的点位, 则围绕网格中心进行旋转, 注意这里是 (x, y), 而不是 (ix, iy)
                 # 注意需要考虑如果指定的网格中心和投影中心不一致的情况
@@ -378,7 +388,7 @@ class model_info_2d(object):
         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):
+    def is_in_domain(self, origin_x, origin_y, use_float=False, extent=None):
         """
         用于判断坐标（经纬度）是否在模式网格范围内
         Update:
@@ -388,7 +398,12 @@ class model_info_2d(object):
             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)
+        if extent is None:
+            xs, xe, ys, ye = 0, self.nx - 1, 0, self.ny - 1
+        else:
+            xs, xe, ys, ye = self.get_select_xy_offset(extent)
+            xe, ye = xe - 1, ye - 1
+        result = (xs <= ix) & (ix <= xe) & (ys <= iy) & (iy <= ye)
         return result
 
     def rotate_xy(self, x, y, rotate_rad=None):
@@ -420,6 +435,49 @@ class model_info_2d(object):
         else:
             ix_new, iy_new = ix, iy
         return ix_new, iy_new
+    
+    def get_select_xy_extent(self, extent: list = [-180, 180, -90, 90]):
+        """
+        根据经纬度范围获取坐标范围
+        """
+        nx, ny = self.nx, self.ny
+        lon_s, lon_e, lat_s, lat_e = extent
+        lon_list = np.concatenate([
+            np.linspace(lon_s, lon_e, nx-1),
+            np.linspace(lon_e, lon_e, ny-1),
+            np.linspace(lon_e, lon_s, nx-1),
+            np.linspace(lon_s, lon_s, ny-1)
+        ])
+        lat_list = np.concatenate([
+            np.linspace(lat_s, lat_s, nx-1),
+            np.linspace(lat_s, lat_e, ny-1),
+            np.linspace(lat_e, lat_e, nx-1),
+            np.linspace(lat_e, lat_s, ny-1)
+        ])
+        x_list, y_list = self.grid_id_float(lon_list, lat_list)
+        xs_float, xe_float, ys_float, ye_float = np.min(x_list), np.max(x_list), np.min(y_list), np.max(y_list)
+        return xs_float, xe_float, ys_float, ye_float
+    
+    def get_select_xy_offset(self, extent: list = [-180, 180, -90, 90]):
+        """
+        根据经纬度范围获取坐标偏移量
+        """
+        nx, ny = self.nx, self.ny
+        xs_float, xe_float, ys_float, ye_float = self.get_select_xy_extent(extent)
+        limit_range = lambda x, vmin, vmax: x + (vmin - x)*(x < vmin) - (x - vmax)*(x > vmax) # x \in [xs, xe]
+        xs, xe, ys, ye = round(xs_float), round(xe_float), round(ys_float), round(ye_float)
+        xs, xe = limit_range(xs, 0, nx-1), limit_range(xe, 0, nx-1)+1
+        ys, ye = limit_range(ys, 0, ny-1), limit_range(ye, 0, ny-1)+1
+        return xs, xe, ys, ye
+
+    
+    def select(self, data, extent: list = [-180, 180, -90, 90]):
+        """
+        根据经纬度范围截取数据
+        """
+        xs, xe, ys, ye = self.get_select_xy_offset(extent)
+        data_select = data[..., ys:ye, xs:xe]
+        return data_select
 
 
 def rotate_xy(xx, yy, cx, cy, rad):
@@ -513,3 +571,19 @@ def from_ctl(file: str) -> model_info_2d:
                     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
+
+def from_ncatts(file: str) -> model_info_2d:
+    """
+    proj4 example:
+        +proj=lcc +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +ellps=sphere +a=6370000 +b=6370000
+        +proj=lonlat
+    """
+    import netCDF4 as nc
+    with nc.Dataset(file) as nf:
+        dx, dy, nx, ny, lowerleft = nf.dx, nf.dy, nf.nx, nf.ny, nf.lowerleft
+        proj = ccrs.CRS(getattr(nf, "proj4", "+latlon"))
+    model = model_info_2d(proj=proj, nx=nx, ny=ny, dx=dx, dy=dy, lowerleft=lowerleft)
+    return model
+
+
+

proj_info.py

@@ -24,8 +24,8 @@ class proj_info(object):
         init=False, wrap=False, rho0=None, nc=None, bigc=None, comp_ll=False, 
         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)经度, 西南角, 度
+        self.lat1       = lat1      # SW latitude (1,1) in degrees (-90->90N) 格点(1, 1)纬度, 西南角, 度, 如果没有指定 knowj, 则为网格中心
+        self.lon1       = lon1      # SW longitude (1,1) in degrees (-180->180E) 格点(1, 1)经度, 西南角, 度, 如果没有指定 knowi, 则为网格中心
         self.lat0       = lat0
         self.lon0       = lon0
         self.dx         = dx        # Grid spacing in meters at truelats, used x方向网格距, m