3

I'm using rasterio's sample function to get values from a (single-layer) raster based on a list of x,y locations:

for val in myraster.sample(locations):
    print(val)

This works fine for the first couple of locations until I get the following error:

IndexError                                Traceback (most recent call last)
<ipython-input-49-9f9c2dcab22c> in <module>
      1 # Sample at these projected locations
----> 2 for val in myraster.sample(locations):
      3     print(val)

~/miniconda3/envs/UHI-NOAA/lib/python3.6/site-packages/rasterio/sample.py in sample_gen(dataset, xy, indexes)
     21             window = Window(col_off, row_off, 1, 1)
     22             data = read(indexes, window=window, masked=False)
---> 23             yield data[:, 0, 0]

IndexError: index 0 is out of bounds for axis 1 with size 0

So I'm assuming there is no data at that point, but I'm trying to understand how to fix that, so here's my current thinking:

  • Is there a way to also print out the location that rasterio is currently looking up?
  • I have also tried to skip the faulty lookup, but a try/except block doesn't help either – the sample method simply stops at that one faulty lookup and skips the rest.

Any leads on how to fix this?

Edit:

I'm on rasterio version 1.0.22, this is what my data looks like:

myraster.bounds

BoundingBox(left=-17619594.54744353, bottom=-6484970.538131511, right=17877405.45255647, top=8750529.46186849)

myraster.meta

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': None,
 'width': 141988,
 'height': 60942,
 'count': 1,
 'crs': CRS.from_wkt('PROJCS["World_Mollweide",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mollweide"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],UNIT["Meter",1.0]]'),
 'transform': Affine(250.0, 0.0, -17619594.54744353,
        0.0, -250.0, 8750529.46186849)}

locations[0:5]

[(-6021233.438259625, 2104298.5232515177),
 (-6020901.044539964, 2106315.744995676),
 (5226731.125781265, 3092960.1175713995),
 (5214407.124651566, 3083679.52532036),
 (5168501.71150224, 2985740.157238809)]

The coordinates in the locations list have been been projected to Mollweide from lat/lon coordinates with pyproj:

inProj = Proj(init='epsg:4326')   # lat/lon 
outProj = Proj('+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs', preserve_flags=True)   # Mollweide

locations = []

# go through the list of stations and project each one to Mollweide
for index, station in stations.iterrows():
    locations.append((transform(inProj,outProj,station['lon'],station['lat'])))

There are actually some locations that are North and South of the raster's bounding box; but IMO that shouldn't produce this error anyway:

xs, ys = zip(*locations)

print(min(xs) > myraster.bounds.left)
True

print(max(xs) < myraster.bounds.right)
True
​

print(min(ys) > myraster.bounds.bottom)
False

print(max(ys) < myraster.bounds.top)
False

  • I just tried with way out of bounds coordinates and just go back a valid "nodata" value. Please edit your question and specify: (1) rasterio version, (2) sample locations coordinates and (3) myraster.bounds and myraster.meta. – user2856 May 22 at 0:47
  • I've made some edits with details. – Carsten May 22 at 6:50
0

Removing the locations that were outside of the raster's bounding box eventually did the trick, even though it seems to me that they should really be producing nodata values (instead of causing an error).

  • Possible bug. I may be misunderstanding the code but it looks like sample is only checking if the coordinates are less than the left and bottom bounds (sample.py#L18). – user2856 May 23 at 9:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.