1

I have been stuck for quite a while on a issue with rasterio. I am attempting to read a netcdf file and mask it with a shapefile. My issue is that it runs and always returns the error:

 /usr/local/lib/python3.6/site-packages/rasterio/mask.py:88: 
 UserWarning: shapes are outside bounds of raster. Are they in 
 different coordinate reference systems?
 warnings.warn('shapes are outside bounds of raster. '

This happens with shape files which are verified to be within the netcdf file. I am assuming there is a small bug in my code but am not sure what could be causing this to happen. Below is the code which is running.

 import os
 import logging
 import sys
 import math
 import json

 from netCDF4 import Dataset
 import numpy as np
 import rasterio
 from rasterio.mask import mask as rio_mask
 from rasterio.mask import raster_geometry_mask as rio_getmask
 from shapely.wkt import loads, dumps
 from shapely.affinity import translate
 from shapely.geometry import mapping, shape
 import shapefile

 def polygon_to_mask(nc, fname, poly, variable):
     """Generates a numpy mask from a polygon"""
     nclons = nc.variables['lon'][:]
     if np.any(nclons > 180):
         poly = translate(poly, xoff=180)

     dst_name = 'NETCDF:"{}":{}'.format(fname, variable)
     with rasterio.open(dst_name, 'r', driver='NetCDF') as raster:
         if raster.transform == rasterio.Affine.identity():
             raise Exception("Unable to determine projection parameters for GDAL "
                         "dataset {}".format(dst_name))
         mask, out_transform, window = rio_getmask(raster, [mapping(poly)], all_touched=True)
     return mask, out_transform, window

 def generatePolygon(reader):
     fields = reader.fields[1:]
     field_names = [field[0] for field in fields]
     buffer = []
     for sr in reader.shapeRecords():
         record = sr.record
         record = [r.decode('utf-8', 'ignore') if isinstance(r, bytes)
                   else r for r in record]
         atr = dict(zip(field_names, record))
         geom = sr.shape.__geo_interface__
         buffer.append(dict(type="Feature", geometry=geom, properties=atr))

     return shape(buffer[0]['geometry'])

 def polygon_to_masked_array(nc, fname, poly, variable):
     """Applies a polygon mask to a variable read from a netCDF file,
     in addition to any masks specified in the file itself (_FillValue)
     Returns a numpy masked array with every time slice masked"""

     mask, out_transform, window  = polygon_to_mask(nc, fname, poly, variable)

     dst_name = 'NETCDF:"{}":{}'.format(fname, variable)
     with rasterio.open(dst_name, 'r', driver='NetCDF') as raster:
         height, width = mask.shape
         out_shape = (raster.count, height, width)

         array = raster.read(window=window, out_shape=out_shape, masked=True)
         array.mask = array.mask | mask

     var = nc.variables[variable]
     scale_factor = getattr(var, 'scale_factor', 1.0)
     add_offset = getattr(var, 'add_offset', 0.0)

     return array * scale_factor + add_offset

 def main():
     reader = shapefile.Reader(os.environ.get('shape'))
     polygon = generatePolygon(reader)
     fname = os.environ.get('file')
     variable = os.environ.get('variable')
     nc = Dataset(fname)
     test = polygon_to_masked_array(nc, fname, polygon, variable)
     print(test)

 main()
  • Ryan, I don't see anything obviously wrong with your code. Can you give me a few more data points? 1. the bounds and CRS of your netcdf dataset, 2. the bounds and CRS of your shapefile, 3. the version of Rasterio that you're using. – sgillies Nov 29 '17 at 7:23
1

Ryan, this is a rasterio bug that is fixed in version 1.0a12 (on PyPI now). pip install --pre -U rasterio.

| improve this answer | |

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.