Possible floating point issue when merging on x,y

[gschlauch]

In overlay_weights.R, I believe I am encountering a floating point issue when merging on "x" and "y" in line 181.

Using

1. grid = the ERA5-Land grid,
2. secondary_raster = Landsat 2010 population estimates, and
3. polygons = 2010 US County shapefile,

w_merged only contains 43,567 rows with non-missing values for "weight". However, when I round "x" and "y" to 7 decimal points in area_weights and weights_dt before merging, this number jumps up to 120,754.

A possible fix could be to round "x" and "y" in area_weights and weights_dt just before merging in line 181.

[tliddell4]

Thank you for pointing this out - this one I need to think about for a bit. I've never worked with ERA5-Land data but apparently the grid layout it is based upon is completely different from that of the ERA5 data.

I imagine you already know all of this, but for the sake of clarity (and solidifying my own newfound knowledge), ERA5 data is based on a rectangular, planar grid which is then "wrapped" around the sphere of the earth, which means each cell is uniform in degrees of lon/lat, but not in area. ERA5-Land data is based on what's called an octahedral reduced gaussian grid, a layout of (somewhat) equally spaced points on the sphere of the earth, which, when unwrapped, do not form a rectangular grid with distinguishable columns. I've pasted a helpful diagram from https://sites.ecmwf.int/docs/atlas/design/grid/ below which illustrates this.

stagg was designed with the rectangular raster grid in mind, and so incorporating other grid types may take a while and may even be outside the scope of the package. I'm not sure if the terra package supports octahedral grids, though if you know the answer to this I would love to hear your thoughts. I'll talk with the team and figure out if we think we could eventually support different grid layouts (I apologize if the argument name grid made it sound like we could), but in the mean time if you wish to use the package with ERA5-land data, I think the best thing to do would be to reproject to a standard rectangular raster grid.

[gschlauch]

Thanks for the response! I didn't know that about the native ERA5-Land grid.

Fortunately, it appears that the ERA5-Land data you can download from the CDS (and therefore through KrigR) have been converted to a rectangular 0.1° × 0.1° lat-lon grid with CRS WGS84. That's what it says on this page at least: "
0.1° x 0.1°; Native resolution is 9 km." I also downloaded a sample file from here with the following options selected and confirmed it is rectangular, 0.1° x 0.1°: 2m temperature, Jan 1 2024, 00:00, Sub-region extraction (default -180 180 -90 90), NetCDF4.

I wasn't having any floating point issues merging the ERA5 data in overlay_weights.R, but I don't think it would hurt to round x and y to 7-8 decimal points before merging to avoid floating point issues with other grids. But I'm not an expert in this area, so feel free to correct me!

[tliddell4]

Oh yes sorry I got caught up considering the gaussian grid (good to know about the ERA5-land data btw) but didn't respond to your original suggestion. I think that's a great idea and I'll let you know when we implement it!

[tliddell4]

Sorry, for the purpose of issue reproduction, was your original post describing what happened when using a gaussian grid or the CDS rectangular grid? If you originally were having trouble with just the rectangular grid, does your more recent post where you say you weren't having floating point issues imply that the NAs were due to something else?

(btw, in working through your example I caught and was able to fix a separate floating point issue when doing 360 %% xres on the landscan data, so thank you for pointing this out to me).

[gschlauch]

The rectangular grid! I was saying that the floating point arose when using the ERA5-Land rectangular grid and Landscan population estimate data, but not the ERA5 rectangular grid and Landscan population estimate data. Does that answer your question?

[gschlauch]

Adding another comment: I think rounding the lat/lons to 6 decimal points is a good number. I tried with 8 decimal points in one application (which I can't remember now) but still had a floating point issue when merging.

[yougsanghvi]

Bumping this issue and adding a little more to it - I've also been using ERA5 Land and has the same issue. It would be nice to have the "tolerance" feature which we have in staggregate in overlay_weights too. However, I wanted to flag that this happened after I used secondary_weights() to create the weights I used in overlay_weights(). I would assume the secondary_weights function should re-project to the grid in the same way but it seems some rounding / floating point mismatches happened? I'm assuming that cause the x-y coords of cropland_10 below don't match those of r_2m_land but when I round both then the match runs just fine. I've added the code snippet below! Let me know if I am doing something wrong on my end too!

secondary weights

cropland_10 <- stagg::secondary_weights(
secondary_raster = cropland_world_2015_era5,
grid = r_2m_land #,
# extent = st_bbox(shapefile_nigeria)
)

county weights

county_weights_modified <- stagg::overlay_weights(
shapefile_nigeria,
polygon_id_col = "GID_2",
grid = r_2m_land,
secondary_weights = cropland_10
)