Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_rasterize_fill_value(basic_geometry, basic_image_2x2):
"""All pixels not covered by shapes should be given fill value."""
default_value = 2
assert np.array_equal(
basic_image_2x2 + 1,
rasterize(
[basic_geometry], out_shape=DEFAULT_SHAPE, fill=1,
default_value=default_value
)
def test_rasterize_missing_out(basic_geometry):
"""If both out and out_shape are missing, should raise exception."""
with pytest.raises(ValueError):
rasterize([basic_geometry], out=None, out_shape=None)
def test_rasterize_geomcollection_no_hole():
"""
Make sure that bug reported in
https://github.com/mapbox/rasterio/issues/1253
does not recur. GeometryCollections are flattened to individual parts,
and should result in no holes where parts overlap.
"""
geomcollection = {'type': 'GeometryCollection', 'geometries': [
{'type': 'Polygon',
'coordinates': (((0, 0), (0, 5), (5, 5), (5, 0), (0, 0)),)},
{'type': 'Polygon',
'coordinates': (((2, 2), (2, 7), (7, 7), (7, 2), (2, 2)),)}
]}
expected = rasterize(geomcollection['geometries'], out_shape=DEFAULT_SHAPE)
assert np.array_equal(
rasterize([geomcollection], out_shape=DEFAULT_SHAPE),
expected
)
def test_rasterize_geometries_symmetric():
"""Make sure that rasterize is symmetric with shapes."""
transform = (1.0, 0.0, 0.0, 0.0, -1.0, 0.0)
truth = np.zeros(DEFAULT_SHAPE, dtype=rasterio.ubyte)
truth[2:5, 2:5] = 1
s = shapes(truth, transform=transform)
result = rasterize(s, out_shape=DEFAULT_SHAPE, transform=transform)
assert np.array_equal(result, truth)
# This potentially creates a large in-memory object so if processing an entire layer it is
# best to explicitly define min/max, especially because its also faster.
if bbox:
x_min, y_min, x_max, y_max = bbox
else:
_bbox, ftrz = min_bbox(ftrz, return_iter=True)
x_min, y_min, x_max, y_max = _bbox
x_delta = x_max - x_min
y_delta = y_max - y_min
cell_size = x_delta / width
height = int(y_delta / cell_size)
if height is 0:
height = 1
output_array = rasterize(
fill=0,
default_value=1,
shapes=(g for g in _geometry_extractor(ftrz)),
out_shape=(height, width),
transform=affine.Affine.from_gdal(*(x_min, cell_size, 0.0, y_max, 0.0, -cell_size)),
all_touched=all_touched,
dtype=rio.uint8
)
# Convert to string dtype and do character replacements
output_array = output_array.astype(np.str_)
if fill is not None and fill != '0':
output_array = np.char.replace(output_array, '0', fill)
if char is not None and fill != '1':
output_array = np.char.replace(output_array, '1', char)
shape = reference_im.shape
if do_transform:
affine_obj = reference_im.transform
# extract geometries and pair them with burn values
if burn_field:
if out_type == 'int':
feature_list = list(zip(df[geom_col],
df[burn_field].astype('uint8')))
else:
feature_list = list(zip(df[geom_col],
df[burn_field].astype('uint8')))
else:
feature_list = list(zip(df[geom_col], [burn_value]*len(df)))
output_arr = features.rasterize(shapes=feature_list, out_shape=shape,
transform=affine_obj)
if out_file:
meta = reference_im.meta.copy()
meta.update(count=1)
if out_type == 'int':
meta.update(dtype='uint8')
with rasterio.open(out_file, 'w', **meta) as dst:
dst.write(output_arr, indexes=1)
return output_arr
tgdf[args.lookup_field][tgdf[field] == key] = dictionary[key]
# Get the template raster
rst = rasterio.open(rst_fn)
# copy and update the metadata from the input raster for the output
meta = rst.meta.copy()
meta.update(compress='lzw')
with rasterio.open(out_fn, 'w+', **meta) as out:
out_arr = out.read(1)
# this is where we create a generator of geom, value pairs to use in rasterizing
shapes = ((geom,value) for geom, value in zip(tgdf.geometry, tgdf[args.lookup_field]))
burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)
out.write_band(1, burned)
else:
out_fn = this_dir+args.fname_prefix+".tif"
print("I am going to proceed using unique values")
print("This will relsult in a new shapefile")
# creating the rasterization column
tgdf['RZ_key'] = ''
#Hosts the equivalences
eqdic = {}
# Fill the null values
tgdf[field].fillna(-5, inplace=True)
# add the center polygons to adj_polygons
adj_polygons.extend([center_polygon])
adj_polygons_class.extend([class_int])
with rasterio.open(sub_image_path) as src:
transform = src.transform
burn_out = np.zeros((src.height, src.width))
# rasterize the shapes
burn_shapes = [(item_shape, item_class_int) for (item_shape, item_class_int) in
zip(adj_polygons, adj_polygons_class)]
#
out_label = rasterize(burn_shapes, out=burn_out, transform=transform,
fill=0, all_touched=False, dtype=rasterio.uint8)
# test: save it to disk
kwargs = src.meta
kwargs.update(
dtype=rasterio.uint8,
count=1)
# width=burn_out.shape[1],
# height=burn_out.shape[0],#transform=transform
# remove nodta in the output
if 'nodata' in kwargs.keys():
del kwargs['nodata']
with rasterio.open(save_path, 'w', **kwargs) as dst:
dst.write_band(1, out_label.astype(rasterio.uint8))
bound_poly = ops.transform(lambda x, y: dst_geotransform * (x, y), box(0., 0., ds_width, ds_height, ccw=False))
if not bound_poly.intersects(dst_poly):
return BrdfTileSummary.empty()
ocean_poly = ops.transform(lambda x, y: fid_mask.transform * (x, y), box(0., 0., fid_mask.width, fid_mask.height))
if not ocean_poly.intersects(dst_poly):
return BrdfTileSummary.empty()
# read ocean mask file for correspoing tile window
# land=1, ocean=0
bound_poly_coords = list(bound_poly.exterior.coords)[:4]
ocean_mask, _ = read_subset(fid_mask, *bound_poly_coords)
ocean_mask = ocean_mask.astype(bool)
# inside=1, outside=0
roi_mask = rasterize([(dst_poly, 1)], fill=0, out_shape=(ds_height, ds_width), transform=dst_geotransform)
roi_mask = roi_mask.astype(bool)
# both ocean_mask and mask shape should be same
if ocean_mask.shape != roi_mask.shape:
raise ValueError('ocean mask and ROI mask do not have the same shape')
if roi_mask.shape != ds.shape:
raise ValueError('BRDF dataset and ROI mask do not have the same shape')
roi_mask = roi_mask & ocean_mask
def layer_sum(param):
layer = ds[param][:, :]
common_mask = roi_mask & (layer != ds.attrs['_FillValue'])
layer = layer.astype('float32')
layer[~common_mask] = np.nan
layer = ds.attrs['scale_factor'] * (layer - ds.attrs['add_offset'])