pySTEPS / pysteps

Python framework for short-term ensemble prediction systems.
https://pysteps.github.io/
BSD 3-Clause "New" or "Revised" License
441 stars 160 forks source link

When running the MRMS data, the following error occurs #361

Closed HY-10 closed 3 months ago

HY-10 commented 3 months ago

Hello, I have a question regarding processing MRMS data. When running the MRMS data, the following error occurs, how should I adjust the parameters?

rainrate, metadata = dimension.aggregate_fields_space(rainrate, metadata, 2000) ValueError: space_window does not equally split R

dnerini commented 3 months ago

hi @HY-10, the method is basically complaining that one of your spatial dimensions (which are assumed to be the last two dimensions of your input) is not a multiple of 2000, see the code for details: https://github.com/pySTEPS/pysteps/blob/master/pysteps/utils/dimension.py#L194

In essence, you have to make sure that the shape of your x and y dimensions, multiplied by the ypixelsize and xpixelsize respectively, is a multiple of 2000.

In your example, what's the shape of rainrate and the content of the metadata dictionnary?

HY-10 commented 3 months ago

I'm sorry, but I still don't understand how to get the appropriate value for the space_window. The MRMS data I'm using has a shape of rainrate.shape=(3,875,1750). Thank you for your guidance.

dnerini commented 3 months ago

what's the resolution of your pixels? can you copy/paste here the content of the metadata dictionary?

HY-10 commented 3 months ago

(3, 875, 1750) {'institution': 'NOAA National Severe Storms Laboratory', 'xpixelsize': 0.04, 'ypixelsize': 0.04, 'unit': 'mm/h', 'accutime': 2.0, 'transform': None, 'zerovalue': 0, 'projection': '+proj=longlat +ellps=IAU76', 'yorigin': 'upper', 'threshold': 0.0125, 'x1': -130.00000000042866, 'x2': -60.000001999571296, 'y1': 20.000000999571306, 'y2': 55.00000000042869, 'cartesian_unit': 'degrees', 'timestamps': array([datetime.datetime(2019, 6, 10, 0, 56), datetime.datetime(2019, 6, 10, 0, 58), datetime.datetime(2019, 6, 10, 1, 0)], dtype=object)}

像素的分辨率是多少?您可以在这里复制/粘贴字典的内容吗?metadata

aperezhortal commented 3 months ago

Hi @HY-10,

Due to the large size of the dataset (3500 x 7000), the importer reduces the original 1kmx1km resolution to 4km x 4km.

Hence, we end with a grid size of (875, 1750). Now, in the metadata, the xpixelsize/ypixelsize is expressed in degrees. That is why you see 0.04 (that is roughly 4km in CONUS).

Now, the aggregate_fields_space can only reduce (or keep) the resolution of the input data. So, let's say you want to downscale the data by 5 (to 20km).

Both 875 and 1750 are divisible by 5. Hence, we use *window=0.045** (remember that the window needs to be in the same units as the pixel size ).

upscaled_rainrate, upscaled_metadata = dimension.aggregate_fields_space(rainrate, metadata, 0.04*5)
print(upscaled_metadata)

{'institution': 'NOAA National Severe Storms Laboratory',
 'xpixelsize': 0.2,
 'ypixelsize': 0.2,
 'unit': 'mm/h',
 'accutime': 2.0,
 'transform': None,
 'zerovalue': 0,
 'projection': '+proj=longlat  +ellps=IAU76',
 'yorigin': 'upper',
 'threshold': 0.0125,
 'x1': -130.00000000042866,
 'x2': -60.000001999571296,
 'y1': 20.000000999571306,
 'y2': 55.00000000042869,
 'cartesian_unit': 'degrees',
 'timestamps': array([datetime.datetime(2019, 6, 10, 0, 0)], dtype=object)}

In which resolution do you need the data? For example, if you need 2km, you need to load the data in the original 1km resolution and then upscale it using a window of pixelsize*2 (0.02).

Let us know if this helps.

Andres

HY-10 commented 3 months ago

Hi @HY-10,

Due to the large size of the dataset (3500 x 7000), the importer reduces the original 1kmx1km resolution to 4km x 4km.

Hence, we end with a grid size of (875, 1750). Now, in the metadata, the xpixelsize/ypixelsize is expressed in degrees. That is why you see 0.04 (that is roughly 4km in CONUS).

Now, the aggregate_fields_space can only reduce (or keep) the resolution of the input data. So, let's say you want to downscale the data by 5 (to 20km).

Both 875 and 1750 are divisible by 5. Hence, we use *window=0.045** (remember that the window needs to be in the same units as the pixel size ).

upscaled_rainrate, upscaled_metadata = dimension.aggregate_fields_space(rainrate, metadata, 0.04*5)
print(upscaled_metadata)

{'institution': 'NOAA National Severe Storms Laboratory',
 'xpixelsize': 0.2,
 'ypixelsize': 0.2,
 'unit': 'mm/h',
 'accutime': 2.0,
 'transform': None,
 'zerovalue': 0,
 'projection': '+proj=longlat  +ellps=IAU76',
 'yorigin': 'upper',
 'threshold': 0.0125,
 'x1': -130.00000000042866,
 'x2': -60.000001999571296,
 'y1': 20.000000999571306,
 'y2': 55.00000000042869,
 'cartesian_unit': 'degrees',
 'timestamps': array([datetime.datetime(2019, 6, 10, 0, 0)], dtype=object)}

In which resolution do you need the data? For example, if you need 2km, you need to load the data in the original 1km resolution and then upscale it using a window of pixelsize*2 (0.02).

Let us know if this helps.

Andres

hi @aperezhortal , Thank you very much! This has been very helpful to me. Now I understand how the value of space_window is obtained.

HY-10 commented 2 months ago

hi@

Hi @HY-10,

Due to the large size of the dataset (3500 x 7000), the importer reduces the original 1kmx1km resolution to 4km x 4km.

Hence, we end with a grid size of (875, 1750). Now, in the metadata, the xpixelsize/ypixelsize is expressed in degrees. That is why you see 0.04 (that is roughly 4km in CONUS).

Now, the aggregate_fields_space can only reduce (or keep) the resolution of the input data. So, let's say you want to downscale the data by 5 (to 20km).

Both 875 and 1750 are divisible by 5. Hence, we use *window=0.045** (remember that the window needs to be in the same units as the pixel size ).

upscaled_rainrate, upscaled_metadata = dimension.aggregate_fields_space(rainrate, metadata, 0.04*5)
print(upscaled_metadata)

{'institution': 'NOAA National Severe Storms Laboratory',
 'xpixelsize': 0.2,
 'ypixelsize': 0.2,
 'unit': 'mm/h',
 'accutime': 2.0,
 'transform': None,
 'zerovalue': 0,
 'projection': '+proj=longlat  +ellps=IAU76',
 'yorigin': 'upper',
 'threshold': 0.0125,
 'x1': -130.00000000042866,
 'x2': -60.000001999571296,
 'y1': 20.000000999571306,
 'y2': 55.00000000042869,
 'cartesian_unit': 'degrees',
 'timestamps': array([datetime.datetime(2019, 6, 10, 0, 0)], dtype=object)}

In which resolution do you need the data? For example, if you need 2km, you need to load the data in the original 1km resolution and then upscale it using a window of pixelsize*2 (0.02).

Let us know if this helps.

Andres

hi @aperezhortal , I want to load the data in the original 1km resolution and then upscale it using a window of pixelsize*2(0.02). So, what should I do

HY-10 commented 2 months ago

嗨,

由于数据集的体积较大 (3500 x 7000),导入器将原来的 1kmx1km 分辨率降低到 4km x 4km。

因此,我们以 (875, 1750) 的网格大小结束。现在,在元数据中,xpixelsize/ypixelsize 以度表示。这就是为什么你会看到 0.04(在 CONUS 中大约是 4 公里)。

现在,aggregate_fields_space只能降低(或保持)输入数据的分辨率。因此,假设您想将数据缩小 5(到 20 公里)。

875 和 1750 都可以被 5 整除。因此,我们使用 *window=0.045**(请记住,窗口需要与像素大小的单位相同)。

upscaled_rainrate, upscaled_metadata = dimension.aggregate_fields_space(rainrate, metadata, 0.04*5)
print(upscaled_metadata)

{'institution': 'NOAA National Severe Storms Laboratory',
 'xpixelsize': 0.2,
 'ypixelsize': 0.2,
 'unit': 'mm/h',
 'accutime': 2.0,
 'transform': None,
 'zerovalue': 0,
 'projection': '+proj=longlat  +ellps=IAU76',
 'yorigin': 'upper',
 'threshold': 0.0125,
 'x1': -130.00000000042866,
 'x2': -60.000001999571296,
 'y1': 20.000000999571306,
 'y2': 55.00000000042869,
 'cartesian_unit': 'degrees',
 'timestamps': array([datetime.datetime(2019, 6, 10, 0, 0)], dtype=object)}

您需要哪种分辨率的数据?例如,如果需要 2km,则需要以原始的 1km 分辨率加载数据,然后使用像素大小*2 (0.02) 的窗口对其进行放大。

如果有帮助,请告诉我们。

安德烈斯

hi @aperezhortal , I want to load the data in the original 1km resolution and then upscale it using a window of pixelsize*2(0.02). So, what should I do?