xarray-contrib / xvec

Vector data cubes for Xarray
https://xvec.readthedocs.io
MIT License
93 stars 9 forks source link

ENH: extract_points from a geospatial raster #36

Closed martinfleis closed 1 year ago

martinfleis commented 1 year ago

Implements a sample_points method designed to generate a subset of the original DataArray/Dataset with N-1 dimensions indexed by an array of points instead of lat/lon (or x/y) coordinates.

Covers the first suggestion from #35

codecov[bot] commented 1 year ago

Codecov Report

Merging #36 (2962704) into main (d2a7dba) will increase coverage by 0.06%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main      #36      +/-   ##
==========================================
+ Coverage   98.98%   99.05%   +0.06%     
==========================================
  Files           3        3              
  Lines         297      317      +20     
==========================================
+ Hits          294      314      +20     
  Misses          3        3              
Impacted Files Coverage Δ
xvec/accessor.py 100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

martinfleis commented 1 year ago

I'd need someone to weigh in on terminology here. If we want to be consistent with stars, we shall use "extract" but there's also rasterio.sample doing the same in Python ecosystem, so for the sake of consistency with that we should use "sample" (as I had until the last commit). Any preference for which to use?

edzer commented 1 year ago

I associate the word sampling with choosing locations, e.g. by random or regular sampling. After you've chosen the locations, you can extract values from the raster data cube at these locations. stars picked the extract terminology from R package raster (now terra), which goes back to 2010.

set.seed(13531)
sample(1:10, 5, replace = TRUE) # pick 5, with replacement
# [1] 10  2  3  2  5
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
st_bbox() |> 
   st_as_sfc() |> 
   st_sample(3, "random") # randomly distributed over the sphere
# Geometry set for 3 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -169.2486 ymin: 1.488203 xmax: -81.12489 ymax: 51.48541
# Geodetic CRS:  WGS 84
# POINT (-81.12489 1.488203)
# POINT (-169.2486 19.72468)
# POINT (-90.56986 51.48541)