How to raise an error that coordinates are out of raster bounds (rasterio.sample)?

Viewed 368

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?

2 Answers

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]

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}")
Related