Растерио: "фигуры выходят за границы растра" с действительными данными

Я довольно долго застрял в проблеме с растерио. Я пытаюсь прочитать файл netcdf и замаскировать его с помощью шейп-файла. Моя проблема в том, что он работает и всегда возвращает ошибку:

 /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. '

Это происходит с файлами форм, которые, как проверено, находятся в файле netcdf. Я предполагаю, что в этом коде есть небольшая ошибка, но я не уверен, что может быть причиной этого. Ниже приведен код, который работает.

 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:
         #based on https://github.com/mapbox/rasterio/blob/master/rasterio/mask.py
         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

     # Weirdly rasterio's mask operation sets, but doesn't respect the
     # scale_factor or add_offset
     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()

0 ответов

Другие вопросы по тегам