diff --git a/cinrad/io/level3.py b/cinrad/io/level3.py index e8996f2..3d939b8 100644 --- a/cinrad/io/level3.py +++ b/cinrad/io/level3.py @@ -208,17 +208,13 @@ def __init__(self, file: Any, product: Optional[str] = None): header["y_grid_num"][0], header["z_grid_num"][0], ) + self.heights = header["height"][0][:zdim] dtype = header["m_data_type"][0] data_size = int(xdim) * int(ydim) * int(zdim) * self.size_conv[dtype] bittype = self.dtype_conv[dtype] data_body = np.frombuffer(f.read(data_size), bittype).astype(int) # Convert data to i4 to avoid overflow in later calculation - if zdim == 1: - # 2D data - out = data_body.reshape(ydim, xdim) - else: - # 3D data - out = data_body.reshape(zdim, ydim, xdim) + out = data_body.reshape(zdim, ydim, xdim) self.data_time = datetime.datetime( header["year"][0], header["month"][0], @@ -227,15 +223,14 @@ def __init__(self, file: Any, product: Optional[str] = None): header["minute"][0], ) # TODO: Recognize correct product name - product_name = ( + self.product_name = ( (b"".join(header["data_name"]).decode("gbk", "ignore").replace("\x00", "")) if not product else product ) - self.product_name = product_name - for pname in ["CR", "3DREF"]: - if product_name.startswith(pname): - self.product_name = pname + for pname in ["CR", "3DREF", "反射率"]: + if pname in self.product_name: + self.product_name = "CR" start_lon = header["start_lon"][0] start_lat = header["start_lat"][0] center_lon = header["center_lon"][0] @@ -246,13 +241,13 @@ def __init__(self, file: Any, product: Optional[str] = None): # y_reso = header['y_reso'][0] self.lon = np.linspace(start_lon, end_lon, xdim) # For shape compatibility self.lat = np.linspace(start_lat, end_lat, ydim) - if self.product_name in ["CR", "3DREF"]: + if self.product_name == "CR": self.data = (np.ma.masked_equal(out, 0) - 66) / 2 else: # Leave data unchanged because the scale and offset are unclear self.data = np.ma.masked_equal(out, 0) - def get_data(self, level: int = 0) -> Dataset: + def get_data(self) -> Dataset: r""" Get radar data with extra information @@ -262,16 +257,13 @@ def get_data(self, level: int = 0) -> Dataset: Returns: xarray.Dataset: Data. """ - dtype = self.product_name - if self.data.ndim == 2: - ret = self.data - else: - ret = self.data[level] - if self.product_name == "3DREF": - dtype = "CR" - da = DataArray(ret, coords=[self.lat, self.lon], dims=["latitude", "longitude"]) + da = DataArray( + self.data, + coords=[self.heights, self.lat, self.lon], + dims=["height", "latitude", "longitude"], + ) ds = Dataset( - {dtype: da}, + {self.product_name: da}, attrs={ "scan_time": self.data_time.strftime("%Y-%m-%d %H:%M:%S"), "site_code": "SWAN", @@ -281,6 +273,8 @@ def get_data(self, level: int = 0) -> Dataset: "elevation": 0, }, ) + if len(self.heights) == 1: + ds = ds.squeeze("height") return ds @@ -518,13 +512,17 @@ def _parse_header(self): if header["magic_number"] != 0x4D545352: raise RadarDecodeError("Invalid standard data") site_config = np.frombuffer(self.f.read(128), SDD_site) - self.code = site_config["site_code"][0].decode().replace("\x00", "") + self.code = ( + site_config["site_code"][0].decode("ascii", "ignore").replace("\x00", "") + ) self.geo = geo = dict() geo["lat"] = site_config["Latitude"] geo["lon"] = site_config["Longitude"] geo["height"] = site_config["ground_height"] task = np.frombuffer(self.f.read(256), SDD_task) - self.task_name = task["task_name"][0].decode().replace("\x00", "") + self.task_name = ( + task["task_name"][0].decode("ascii", "ignore").replace("\x00", "") + ) cut_num = task["cut_number"][0] self.scan_config = np.frombuffer(self.f.read(256 * cut_num), SDD_cut) ph = np.frombuffer(self.f.read(128), SDD_pheader) @@ -572,7 +570,7 @@ def _parse_radial_fmt(self): az[az > 360] -= 360 azi = az * deg2rad # self.azi = np.deg2rad(azi) - dist = np.arange(start_range + reso, end_range + reso, reso) + dist = np.arange(start_range // reso + 1, end_range // reso + 1, 1) * reso lon, lat = get_coordinate( dist, azi, self.params["elevation"], self.stationlon, self.stationlat ) @@ -1061,6 +1059,14 @@ class MocMosaic(object): 解析中国气象局探测中心-天气雷达拼图系统V3.0-产品 """ + dtype_conv = [ + ("QREF", "REF"), + ("CREF", "CR"), + ("CRF", "CR"), + ("CAP", "REF"), + ("UCR", "CR"), + ] + def __init__(self, file): self.f = prepare_file(file) radartype = self._check_ftype() @@ -1151,6 +1157,10 @@ def _parse_single(self): r = np.reshape(out_data, (len(out_data), self.ny, self.nx)) ret = np.flipud(r) # Data store reverse in y axis data = (np.ma.masked_less(ret, 0)) / self.scale + for oname, cname in self.dtype_conv: + if self.dtype == oname: + self.dtype = cname + break da = DataArray( data, coords=[heights, self.lat, self.lon], @@ -1170,6 +1180,8 @@ def _parse_single(self): "task": "VCP{}".format(self.vcp), }, ) + if len(heights) == 1: + ds = ds.squeeze("height") self._dataset = ds def _parse_mosaic(self): @@ -1198,7 +1210,7 @@ def _parse_mosaic(self): ) self.lon = np.linspace(edge_w, edge_e, nx) self.lat = np.linspace(edge_s, edge_n, ny) - databody = self.f.read(nx * ny * 2) + databody = self.f.read() if compress == 0: databody = databody elif compress == 1: @@ -1210,6 +1222,10 @@ def _parse_mosaic(self): out = np.frombuffer(databody, "i2").astype(int).reshape(ny, nx) out = np.flipud(out) self.data = (np.ma.masked_less(out, 0)) / scale + for oname, cname in self.dtype_conv: + if self.dtype == oname: + self.dtype = cname + break da = DataArray( self.data, coords=[self.lat, self.lon], dims=["latitude", "longitude"] )