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_utm(self):
with create_tmp_geotiff() as (tmp_file, expected):
with xr.open_rasterio(tmp_file) as rioda:
assert_allclose(rioda, expected)
assert rioda.attrs["scales"] == (1.0, 1.0, 1.0)
assert rioda.attrs["offsets"] == (0.0, 0.0, 0.0)
assert rioda.attrs["descriptions"] == ("d1", "d2", "d3")
assert rioda.attrs["units"] == ("u1", "u2", "u3")
assert isinstance(rioda.attrs["crs"], str)
assert isinstance(rioda.attrs["res"], tuple)
assert isinstance(rioda.attrs["is_tiled"], np.uint8)
assert isinstance(rioda.rio._cached_transform(), Affine)
np.testing.assert_array_equal(
rioda.attrs["nodatavals"], [np.NaN, np.NaN, np.NaN]
)
# Check no parse coords
with xr.open_rasterio(tmp_file, parse_coordinates=False) as rioda:
assert "x" not in rioda.coords
def test_rasterio_environment(self):
with create_tmp_geotiff() as (tmp_file, expected):
# Should fail with error since suffix not allowed
with pytest.raises(Exception):
with rasterio.Env(GDAL_SKIP="GTiff"):
with xr.open_rasterio(tmp_file) as actual:
assert_allclose(actual, expected)
def test_no_mftime(self):
# rasterio can accept "filename" urguments that are actually urls,
# including paths to remote files.
# In issue #1816, we found that these caused dask to break, because
# the modification time was used to determine the dask token. This
# tests ensure we can still chunk such files when reading with
# rasterio.
with create_tmp_geotiff(
8, 10, 3, transform_args=[1, 2, 0.5, 2.0], crs="+proj=latlong"
) as (tmp_file, expected):
with mock.patch("os.path.getmtime", side_effect=OSError):
with xr.open_rasterio(tmp_file, chunks=(1, 2, 2)) as actual:
assert isinstance(actual.data, da.Array)
assert_allclose(actual, expected)
def readImage(image_folder_fp):
print("Reading Sentinel image from: %s" % (image_folder_fp))
### Rather than figuring out what the filepath inside SAFE folder is, this is just finding the red and nir files with correct endings
for subdir, dirs, files in os.walk(image_folder_fp):
for file in files:
if file.endswith("_B04_10m.jp2"):
red_fp = os.path.join(subdir,file)
if file.endswith("_B08_10m.jp2"):
nir_fp = os.path.join(subdir,file)
### Read the red and nir band files to xarray and with the chunk-option to dask
red = xr.open_rasterio(red_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
nir = xr.open_rasterio(nir_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
### Scale the image values back to real reflectance values
red = red /10000
nir = nir /10000
return red,nir
def read(self):
"""Read the image."""
dataset = rasterio.open(self.finfo['filename'])
# Create area definition
if hasattr(dataset, 'crs') and dataset.crs is not None:
self.area = utils.get_area_def_from_raster(dataset)
data = xr.open_rasterio(dataset, chunks=(1, CHUNK_SIZE, CHUNK_SIZE))
attrs = data.attrs.copy()
# Rename to Satpy convention
data = data.rename({'band': 'bands'})
# Rename bands to [R, G, B, A], or a subset of those
data['bands'] = BANDS[data.bands.size]
# Mask data if alpha channel is present
try:
data = mask_image_data(data)
except ValueError as err:
logger.warning(err)
data.attrs = attrs
self.file_content['image'] = data
S=6205900.
W=344520.
E=344620.
plot_bbox = np.asarray([[W,N],[E,N],[E,S],[W,S]])
pts, starting_ids, trees = io.load_lidar_data_by_polygon(las_list,plot_bbox,laz_files=False,max_pts_per_tree = 5*10**5)
N_trees = len(trees)
# filter LiDAR data by return classification
pts[np.any((pts[:,4]==3,pts[:,4]==4,pts[:,4]==5),axis=0),4]=1
pts[pts[:,2]<0,2]=0
# Load sensitivity analysis to spatial resolution
resolution_sensitivity = np.load('SEOS_MH_sensitivity_resolution.npy')[()]
density_sensitivity = np.load('SEOS_MH_sensitivity_pulse_density.npy')[()]
# Load 3D maps of canopy density
pad = xr.open_rasterio(pad_file)
# other parameters
max_height = 40.
layer_thickness = 0.5
heights = np.arange(0.,max_height,layer_thickness)+0.5
n_layers = heights.size
plot_width = 100.
sample_res = 5.
kappa = 1.
n_iter = 100
# bin lidar returns
heights,first_return_profile,n_ground_returns = LAD.bin_returns(pts, max_height, layer_thickness)
"""
PART B: Plot canopy profiles
def open(self, path, fmt, with_nodata=False, isel_band=None):
"Open raster and return in requested format"
if fmt is None or fmt.lower() == 'xarray':
xarr = xr.open_rasterio(path, parse_coordinates=True)
if isel_band is not None:
xarr = xarr.isel(band=0)
if with_nodata:
xarr = convert_nodata_to_nans(xarr)
return xarr
raster_data = self.read(path)
if fmt.lower() == 'rasterio':
return raster_data
if fmt.lower() == 'array':
return raster_data.read()
raise NotImplementedError('format %s not recognized' % fmt)
def readImage(image_folder_fp):
print("Reading Sentinel image from: %s" % (image_folder_fp))
### Rather than figuring out what the filepath inside SAFE folder is, this is just finding the red and nir files with correct endings
for subdir, dirs, files in os.walk(image_folder_fp):
for file in files:
if file.endswith("_B04_10m.jp2"):
red_fp = os.path.join(subdir,file)
if file.endswith("_B08_10m.jp2"):
nir_fp = os.path.join(subdir,file)
### Read the red and nir band files to xarray and with the chunk-option to dask
red = xr.open_rasterio(red_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
nir = xr.open_rasterio(nir_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
### Scale the image values back to real reflectance values
red = red /10000
nir = nir /10000
return red,nir
Using rasterio's projection information for more accurate plots.
This example extends :ref:`recipes.rasterio` and plots the image in the
original map projection instead of relying on pcolormesh and a map
transformation.
"""
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import xarray as xr
# Read the data
url = "https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif"
da = xr.open_rasterio(url)
# The data is in UTM projection. We have to set it manually until
# https://github.com/SciTools/cartopy/issues/813 is implemented
crs = ccrs.UTM("18N")
# Plot on a map
ax = plt.subplot(projection=crs)
da.plot.imshow(ax=ax, rgb="band", transform=crs)
ax.coastlines("10m", color="r")
plt.show()