Dimensions error when filtering #700
Comments
|
It seems that applying a filter with only single pressure level works, but when combined with filtering by date, we get this error. When I swapped the order of the filters, I got this error: If you are looking for an alternative approach, filtering with square brackets can be applied this way: plev = st_get_dimension_values(hu, "plev")
time = st_get_dimension_values(hu, "time")
plev_idx = which(plev == units::set_units(850*100, Pa))
time_idx = which(time <= as.POSIXct("1950-06-01 12:00:00", tz = "UTC"))
hu[,,, plev_idx, time_idx]
# stars object with 4 dimensions and 1 attribute
# attribute(s), summary of first 1e+05 cells:
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# Relative Humidity [%] 0.03111257 48.86616 70.07285 63.57896 81.61889 124.086
# dimension(s):
# from to offset delta refsys values x/y
# x 1 512 -0.3516 0.7031 NA NULL [x]
# y 1 256 89.81 -0.7017 NA NULL [y]
# plev 2 2 NA NA udunits [85000,70000) [Pa]
# time 1 152 1950-01-01 12:00:00 UTC 1 days POSIXct NULL |
|
Yes, this solution is working. Thanks!! |
|
I just noticed that it works only without using stars_proxy. |
|
This seems to now work: library(stars)
hu <- read_mdim("hur_day_EC-Earth3_historical_r1i1p1f1_gr_19500101-19501231.nc", proxy = T)
plev = st_get_dimension_values(hu, "plev")
time = st_get_dimension_values(hu, "time")
plev_idx = which(plev == units::set_units(850*100, Pa))
time_idx = which(time <= as.Date("1950-06-01"))
hu <- hu[,,, plev_idx, time_idx, drop = T]
write_mdim(hu, "/tmp/test.nc") |
|
Yes, it worked. Thanks. But as far as I see, it doesn't work with |
Can you give an example where that doesn't work? |
|
The example you pasted here first (I saw it only in my email) is what I searched for, but as you say, I would need it as a proxy. I tried out Another issue I found now is in Thank you for your help. |
|
Yes, but |
hu <- st_set_dimensions(hu, "lon", values = st_get_dimension_values(hu, "lon")-180)now generates an error, as dimensions are being reread and overwritten in > hu <- read_stars("hur_day_EC-Earth3_historical_r1i1p1f1_gr_19500101-19501231.nc", proxy = T)
Warning messages:
1: In CPL_get_metadata(file, NA_character_, options) :
GDAL Message 1: Latitude grid not spaced evenly. Setting projection for grid spacing is within 0.1 degrees threshold.
2: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver), :
GDAL Message 1: Latitude grid not spaced evenly. Setting projection for grid spacing is within 0.1 degrees threshold.
> st_crop(hu, st_bbox(c(xmin = 20, ymin = 25, xmax = 40, ymax = 50))) |> st_as_stars()
stars object with 4 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
Min. 1st Qu. Median Mean 3rd Qu. Max.
Relative Humidity [%] 0.002519234 0.7560223 19.43109 31.13065 58.62669 114.6512
dimension(s):
from to offset delta refsys
x 29 58 -0.3516 0.7031 NA
y 57 93 89.81 -0.7017 NA
plev 1 8 NA NA udunits
time 1 365 1950-01-01 12:00:00 UTC 1 days POSIXct
values x/y
x NULL [x]
y NULL [y]
plev [1e+05,85000) [Pa],...,[1000,-3000) [Pa]
time NULL
Warning messages:
1: In CPL_get_metadata(file, NA_character_, options) :
GDAL Message 1: Latitude grid not spaced evenly. Setting projection for grid spacing is within 0.1 degrees threshold.
2: In CPL_read_gdal(as.character(x), as.character(options), as.character(driver), :
GDAL Message 1: Latitude grid not spaced evenly. Setting projection for grid spacing is within 0.1 degrees threshold. |
|
In this case, which is the best way to change from longitude 360 to 180? I still get the following error when I include the rest of the steps even with With read_stars(): read_mdim() using 360º lon, both errors are the same with or without the steps |
Arguably no, because you loose information (the value of that single level). You can put an
Yes, I'll look into either fixing this or making the error message more helpful.
Tell the data producers to do so? |
|
More seriously,
I managed to do so outside R, with cdo, cdo -sellonlatbox,-180,180,-90,90 hur_day_EC-Earth3_historical_r1i1p1f1_gr_19500101-19501231.nc out.ncthen r = read_mdim("out.nc")
plot(adrop(r[,,,1,1]),axes = TRUE)gives But to be honest, it took me quite some time to get |
|
Ok. Thanks. Terra has the rotate() function for this kind of operation. Regarding telling the data producers to do so, climate projections are run with 360º by default as standard. |
Yes, I discussed this with @Nowosad over lunch, and given that these are all discrete global grids, this whole shifting or rotating should not be needed; the reason it's needed is that much of our software stacks still tend to treat geographic coordinates as Cartesian coordinates. I think that needs to change. (Terra will also not preserve the irregular latitude values.) |
|
I found another possibility, but as proxy_stars, it doesn't work. For me, the grid is regular, but the only difference is in longitude 360º. I thought that I could use |
I get a dimension error when I want to filter a single pressure level from my ncdf. I can filter the time dimension without any issues. Data example can be downloaded here.
The text was updated successfully, but these errors were encountered: