Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""
Create Cloud Optimized Geotiff.
Parameters
----------
src_path : str or PathLike object
A dataset path or URL. Will be opened in "r" mode.
This script is the rasterio equivalent of
https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/validate_cloud_optimized_geotiff.py
"""
errors = []
details = {}
if not GDALVersion.runtime().at_least("2.2"):
raise Exception("GDAL 2.2 or above required")
with rasterio.open(src_path) as src:
if not src.driver == "GTiff":
raise Exception("The file is not a GeoTIFF")
filelist = [os.path.basename(f) for f in src.files]
src_bname = os.path.basename(src_path)
if len(filelist) > 1 and src_bname + ".ovr" in filelist:
errors.append(
"Overviews found in external .ovr file. They should be internal"
)
if src.width >= 512 or src.height >= 512:
if not src.is_tiled:
errors.append(
"The file is greater than 512xH or 512xW, but is not tiled"
)
def get_projection(realpath, path):
with rasterio.open(os.path.join(str(realpath), str(path))) as img:
left, bottom, right, top = img.bounds
spatial_reference = str(str(getattr(img, 'crs_wkt', None) or img.crs.wkt))
geo_ref_points = {
'ul': {'x': left, 'y': top},
'ur': {'x': right, 'y': top},
'll': {'x': left, 'y': bottom},
'lr': {'x': right, 'y': bottom},
}
return geo_ref_points, spatial_reference
def from_geotiff(path_to_tif):
with rasterio.drivers(CPL_DEBUG=True):
with rasterio.open(path_to_tif) as src:
b, g, r = src.read()
total = np.zeros(r.shape, dtype=rasterio.uint16)
for band in r, g, b:
total += band
total /= 3
return total, src
import pprint
import rasterio
from rasterio.features import shapes
with rasterio.open('tests/data/shade.tif') as src:
image = src.read(1)
# Print the first two shapes...
pprint.pprint(
list(shapes(image))[:2]
)
def read_oneband_image_to_1dArray(image_path,nodata=None, ignore_small=None):
if os.path.isfile(image_path) is False:
raise IOError("error, file not exist: " + image_path)
with rasterio.open(image_path) as img_obj:
# read the all bands (only have one band)
indexes = img_obj.indexes
if len(indexes) != 1:
raise IOError('error, only support one band')
data = img_obj.read(indexes)
data_1d = data.flatten() # convert to one 1d, row first.
if nodata is not None:
data_1d = data_1d[data_1d != nodata]
if ignore_small is not None:
data_1d = data_1d[data_1d >= ignore_small ]
return data_1d
def get_projection(path):
with rasterio.open(str(path)) as img:
left, bottom, right, top = img.bounds
return {
'spatial_reference': str(str(getattr(img, 'crs_wkt', None) or img.crs.wkt)),
'geo_ref_points': {
'ul': {
'x': left,
'y': top
},
'ur': {
'x': right,
'y': top
},
'll': {
'x': left,
'y': bottom
},
valid_samples: array-like
Numpy array of extracted raster values, typically 2d
valid_coordinates: 2d array-like
2D numpy array of xy coordinates of extracted values
"""
# set the seed
np.random.seed(seed=random_state)
# determine number of GDAL rasters to sample
n_features = len(rasters)
# open rasters
dataset = [0] * n_features
for n in range(n_features):
dataset[n] = rasterio.open(rasters[n], mode='r')
# create np array to store randomly sampled data
# we are starting with zero initial rows because data will be appended,
# and number of columns are equal to n_features
valid_samples = np.zeros((0, n_features))
valid_coordinates = np.zeros((0, 2))
# loop until target number of samples is satified
satisfied = False
n = size
while satisfied is False:
# generate random row and column indices
Xsample = np.random.choice(range(0, dataset[0].shape[1]), n)
Ysample = np.random.choice(range(0, dataset[0].shape[0]), n)
def estimateModel():
# Open image files of predicted data and test data
with rasterio.open(predicted_image_output_path, 'r') as prediction_dataset:
with rasterio.open(test_image_path, 'r') as test_labels_dataset:
#Find out the overlappin area of two images.
#Because of tiling the prediction image is slightly smaller than the original clip.
left = max(prediction_dataset.bounds.left,test_labels_dataset.bounds.left)
bottom = max(prediction_dataset.bounds.bottom,test_labels_dataset.bounds.bottom)
right = min(prediction_dataset.bounds.right,test_labels_dataset.bounds.right)
top = min(prediction_dataset.bounds.top,test_labels_dataset.bounds.top)
common_bbox = {
"type": "Polygon",
"coordinates": [[
[left, bottom],
[left, top],
[right, top],
[right, bottom],
crs = rasterio.crs.CRS.from_epsg(4326) # WGS 84,
x_left_top = xll_lon
y_left_top = yll_lon + res*height
profile = {'driver': 'GTiff', 'nodata': no_data, 'width': width, 'height': height, 'count': 1,
'crs': crs,
'transform': (x_left_top, res, 0.0, y_left_top , 0.0, -res),
}
# And then change the band count to 1, set the
# dtype to uint8, and specify LZW compression.
profile.update(
dtype=rasterio.float32,
count=1,
compress='lzw')
with rasterio.open(output, 'w', **profile) as dst:
dst.write(array_2d.astype(rasterio.float32), 1)
def get_model_coords(self):
"""Get LFP model coordinates for 1D and 2D mesh via BMI. The LFP model
should be initialized first in order to access the variables."""
logger.info('Getting LFP model coordinates.')
i_ind, j_ind = np.where(np.logical_and(self.get_var('SGCwidth') > 0., self.get_var('DEM') != -9999))
print i_ind.shape, j_ind.shape
fn_map = join(self.model_data_dir,
self.model_config['dummysection']['DEMfile'])
if not isfile(fn_map):
raise IOError('landmask file not found')
self._fn_landmask = fn_map
with rasterio.open(fn_map, 'r') as ds:
self.grid_index = ds.index
self.model_grid_res = ds.res
self.model_grid_bounds = ds.bounds
self.model_grid_shape = ds.shape
self.model_grid_transform = ds.transform
list_x_coords, list_y_coords = self.model_grid_transform * (i_ind, j_ind)
self.model_1d_coords = zip(list_x_coords, list_y_coords)
self.model_1d_indices = zip(i_ind, j_ind)
pass