Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inconsistent treatment of units in st_as_sf() #652

Closed
fgoerlich opened this issue Oct 4, 2023 · 3 comments
Closed

Inconsistent treatment of units in st_as_sf() #652

fgoerlich opened this issue Oct 4, 2023 · 3 comments

Comments

@fgoerlich
Copy link

Hi,
H1084 is an SpatRaster. Executing

H1084 |>
  crop(filter(MTN50grid, MTN50 == "1084-1086")) |>
  st_as_stars() |>
  st_xy2sfc(as_points = TRUE) |>
  st_as_sf() |>
  rename(H1084 = 1)

produces numeric, without units, output

Simple feature collection with 155 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 607385 ymin: 3180900 xmax: 617162.5 ymax: 3182725
Projected CRS: ETRS89 / UTM zone 28N
First 10 features:
   H1084                 geometry
1    4.1   POINT (615690 3182520)
2    2.8 POINT (615692.5 3182520)
3    4.0   POINT (615685 3182518)
4    3.9 POINT (615687.5 3182518)
5    4.0   POINT (615690 3182518)
6    4.1 POINT (615692.5 3182518)
7    0.5 POINT (615697.5 3182518)
8    4.3 POINT (615682.5 3182515)
9    4.0   POINT (615685 3182515)
10   3.8 POINT (615687.5 3182515)

Without the cropping (from terra) instruction,

H0184 |>
  st_as_stars() |>
  st_xy2sfc(as_points = TRUE) |>
  st_as_sf() |>
  rename(H1084 = 1)

I get numeric, with units, output

Simple feature collection with 238404 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 484497.5 ymin: 4686732 xmax: 512062.5 ymax: 4705362
Projected CRS: ETRS89 / UTM zone 29N
First 10 features:
       H1084                 geometry
1  4.735 [m]   POINT (510505 4705295)
2  5.330 [m] POINT (510507.5 4705295)
3  2.100 [m] POINT (510492.5 4705293)
4  6.137 [m]   POINT (510505 4705293)
5  6.900 [m] POINT (510507.5 4705293)
6  2.723 [m]   POINT (510490 4705290)
7  2.005 [m] POINT (510492.5 4705290)
8  2.789 [m]   POINT (510450 4705288)
9  4.000 [m] POINT (510452.5 4705288)
10 5.750 [m] POINT (510502.5 4705288)

I wonder what makes the difference, since croping produces an SpatRaster also.
PS: I can provide the data if necessary.
Thanks

@edzer
Copy link
Member

edzer commented Oct 4, 2023

Speculating: H1084 is a pretty large file, which terra doesn't read into memory. On convertion to stars, stars gets to read the file and figures out it has units (while, for instance, it is a NetCDF file or so). If terra crops this, it reads a section of the file, and ignores the units. So the question is: is this a stars issue or a terra issue?

@fgoerlich
Copy link
Author

Probably you are rigth.
If I read H1084, which is a raster (tif) object, as star object, I always get units.
stars reads the units, but apparently terra does not. I am not sure if this is a stars issue or a terra issue.

Find the data an small example.

Units.zip

@edzer
Copy link
Member

edzer commented Oct 4, 2023

Confirmed.

@edzer edzer closed this as completed Oct 4, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants