2

I'm using rasterio sample module, list_of_coords are just 3 random points but only 2nd is within the raster bounds:

list_of_coords = [(754.2,4248548.6), (754222.6,4248548.6), (54.9,4248548.4)]
sample = np.array(list(rasterio.sample.sample_gen(raster, list_of_coords))).flatten() 

output:

[  0 896   0]

It works great but as you can see if coordinates are out of raster image it gives value 0. Is there any way to let user know that coords that they put in the list are out of raster bounds? 0 can be also a value for existing within raster bounds point so simple loop:

for idx, element in enumerate(sample):
    if element == 0:
        print(f"coords {a[idx]} out of raster")

is not a good solution. Here is what I came up with so far:

Knowing basic info about geographic coordinate system and raster bounds we could write down some "rules". With raster.bounds I got bbox for my raster and I wrote a bit better loop:

for idx, element in enumerate(a):
    if element[0] > 0 and band2.bounds.right > 0 and element[0] > band2.bounds.right\
       or element[0] > 0 and band2.bounds.left > 0 and element[0] < band2.bounds.left: #more conditions
       print(f"coords {a[idx]} out of raster")

output (correct):

coords (754.6, 4248548.6) out of raster
coords (54.6, 4248548.6) out of raster

Problem is - to cover all posibilities I need to write in this loop way more conditions, is tere a better way to let user know that given point is out of raster?

Edyficjum
  • 81
  • 6
  • Do you want things to stop? You can say `raise ValueError(f"coords {a[idx]} out of raster")`. – Tim Roberts Jun 08 '21 at 19:28
  • @TimRoberts No I'm working with spatial data and I want to know if there is an efficient way (better than my solution) to give user an information that coordinates given in the list are out of raster bounds. So for example, my raster bounds are (x1 = 2, x2=6, y1 = 3, y2 = 7) and user gives coords for point (100,100). That point is not in raster bounds and I'm trying to find a way to inform user about that. Rasterio sample module gives back value 0 in case like that, but that is not good (because if raster image has values 0-255) 0 is in that range. – Edyficjum Jun 08 '21 at 19:40

2 Answers2

3

rasterio.sample.sample_gen provide a masked argument. When True it yield Masked arrays according to bounding box of raster dataset.

>>> import rasterio
>>> ds = rasterio.open("raster.tif")
>>> ds.bounds
BoundingBox(left=-0.0001388888888888889, bottom=40.999861111111116, right=1.000138888888889, top=42.00013888888889)
>>> # Inside bbox
>>> next(rasterio.sample.sample_gen(ds, ((0.5, 41.5), ), masked=True))
masked_array(data=[130],
             mask=False,       # <= No mask (ndim=0)
       fill_value=999999,
            dtype=int16)
>>> # Outside bbox
>>> next(rasterio.sample.sample_gen(ds, ((0, 0), ), masked=True))
masked_array(data=[0],
             mask=[False],     # <= Mask ndim=1
       fill_value=999999,
            dtype=int16)

And them convert to python list with None when coordinates are outside raster bounds:

>>> [None if x.mask.ndim == 1 and not x.mask[0] else x[0]
...  for x in rasterio.sample.sample_gen(ds, ((0.5, 41.5), (0, 0)), masked=True)]
[130, None]
Balaïtous
  • 826
  • 6
  • 9
0

The shortest piece of code I can think of. It's similar to how the sample function checks for masking https://github.com/mapbox/rasterio/blob/master/rasterio/sample.py#L46

def sample_or_error(raster_file, x, y):
    with rasterio.open(raster_file) as src:
        row, col = src.index(x, y)
        if any([row < 0, col < 0, row >= src.height, col >= src.width]):
            raise ValueError(f"({lon}, {lat}) is out of bounds from {src.bounds}")
Scott Staniewicz
  • 712
  • 6
  • 14