本文整理汇总了Python中rasterio.warp.calculate_default_transform方法的典型用法代码示例。如果您正苦于以下问题:Python warp.calculate_default_transform方法的具体用法?Python warp.calculate_default_transform怎么用?Python warp.calculate_default_transform使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类rasterio.warp
的用法示例。
在下文中一共展示了warp.calculate_default_transform方法的11个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: get_max_zoom
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_max_zoom(self, snap=0.5, max_z=23):
"""Calculate raster max zoom level."""
dst_affine, w, h = calculate_default_transform(
self.crs,
"epsg:3857",
self.meta["width"],
self.meta["height"],
*self.crs_bounds
)
res_max = max(abs(dst_affine[0]), abs(dst_affine[4]))
tgt_z = max_z
mpp = 0.0
# loop through the pyramid to file the closest z level
for z in range(1, max_z):
mpp = _meters_per_pixel(z, 0)
if (mpp - ((mpp / 2) * snap)) < res_max:
tgt_z = z
break
return tgt_z
示例2: ConvertRaster2LatLong
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def ConvertRaster2LatLong(InputRasterFile,OutputRasterFile):
"""
Convert a raster to lat long WGS1984 EPSG:4326 coordinates for global plotting
MDH
"""
# import modules
import rasterio
from rasterio.warp import reproject, calculate_default_transform as cdt, Resampling
# read the source raster
with rasterio.open(InputRasterFile) as src:
#get input coordinate system
Input_CRS = src.crs
# define the output coordinate system
Output_CRS = {'init': "epsg:4326"}
# set up the transform
Affine, Width, Height = cdt(Input_CRS,Output_CRS,src.width,src.height,*src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': Output_CRS,
'transform': Affine,
'affine': Affine,
'width': Width,
'height': Height
})
with rasterio.open(OutputRasterFile, 'w', **kwargs) as dst:
for i in range(1, src.count+1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.affine,
src_crs=src.crs,
dst_transform=Affine,
dst_crs=Output_CRS,
resampling=Resampling.bilinear)
示例3: get_max_zoom
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_max_zoom(src_dst, lat=0.0, tilesize=256):
"""
Calculate raster max zoom level.
Parameters
----------
src: rasterio.io.DatasetReader
Rasterio io.DatasetReader object
lat: float, optional
Center latitude of the dataset. This is only needed in case you want to
apply latitude correction factor to ensure consitent maximum zoom level
(default: 0.0).
tilesize: int, optional
Mercator tile size (default: 256).
Returns
-------
max_zoom: int
Max zoom level.
"""
dst_affine, w, h = calculate_default_transform(
src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
)
native_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
# Correction factor for web-mercator projection latitude distortion
latitude_correction_factor = math.cos(math.radians(lat))
corrected_resolution = native_resolution * latitude_correction_factor
max_zoom = zoom_for_pixelsize(corrected_resolution, tilesize=tilesize)
return max_zoom
示例4: _get_vrt
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def _get_vrt(src: DatasetReader, rs_method: int) -> WarpedVRT:
from terracotta.drivers.raster_base import RasterDriver
target_crs = RasterDriver._TARGET_CRS
vrt_transform, vrt_width, vrt_height = calculate_default_transform(
src.crs, target_crs, src.width, src.height, *src.bounds
)
vrt = WarpedVRT(
src, crs=target_crs, resampling=rs_method, transform=vrt_transform,
width=vrt_width, height=vrt_height
)
return vrt
示例5: get_min_zoom
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_min_zoom(self, snap=0.5, max_z=23):
"""Calculate raster min zoom level."""
dst_affine, w, h = calculate_default_transform(
self.crs,
"epsg:3857",
self.meta["width"],
self.meta["height"],
*self.crs_bounds
)
res_max = max(abs(dst_affine[0]), abs(dst_affine[4]))
max_decim = self.overiew_levels[-1]
resolution = max_decim * res_max
tgt_z = 0
mpp = 0.0
# loop through the pyramid to file the closest z level
for z in list(range(0, 24))[::-1]:
mpp = _meters_per_pixel(z, 0)
tgt_z = z
if (mpp - ((mpp / 2) * snap)) > resolution:
break
return tgt_z
示例6: reprojectedRaster
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def reprojectedRaster(rasterFn,ref_vectorFn):
dst_crs=gpd.read_file(ref_vectorFn).crs
print(dst_crs) #{'init': 'epsg:4326'}
dst_raster_projected=os.path.join(dataFp_1,r"svf_dstRasterProjected_b.tif")
a_T = datetime.datetime.now()
# dst_crs='EPSG:4326'
with rasterio.open(rasterFn) as src:
transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height,
# 'compress': "LZW",
'dtype':rasterio.float32,
})
# print(src.count)
with rasterio.open(dst_raster_projected, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest
)
b_T = datetime.datetime.now()
print("reprojected time span:", b_T-a_T)
#根据Polgyon统计raster栅格信息
示例7: reprojectedRaster
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def reprojectedRaster(rasterFn,ref_vectorFn,dst_raster_projected):
dst_crs=gpd.read_file(ref_vectorFn).crs
print(dst_crs) #{'init': 'epsg:4326'}
a_T = datetime.datetime.now()
# dst_crs='EPSG:4326'
with rasterio.open(rasterFn) as src:
transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height,
# 'compress': "LZW",
'dtype':rasterio.uint8, #rasterio.float32
})
# print(src.count)
with rasterio.open(dst_raster_projected, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest
)
b_T = datetime.datetime.now()
print("reprojected time span:", b_T-a_T)
#根据Polgyon统计raster栅格信息
示例8: get_zooms
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_zooms(
src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
ensure_global_max_zoom: bool = False,
tilesize: int = 256,
) -> Tuple[int, int]:
"""
Calculate raster min/max mercator zoom level.
Parameters
----------
src_dst: rasterio.io.DatasetReader
Rasterio io.DatasetReader object
ensure_global_max_zoom: bool, optional
Apply latitude correction factor to ensure max_zoom equality for global
datasets covering different latitudes (default: False).
tilesize: int, optional
Mercator tile size (default: 256).
Returns
-------
min_zoom, max_zoom: Tuple
Min/Max Mercator zoom levels.
"""
bounds = transform_bounds(
src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21
)
center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
lat = center[1] if ensure_global_max_zoom else 0
dst_affine, w, h = calculate_default_transform(
src_dst.crs,
constants.WEB_MERCATOR_CRS,
src_dst.width,
src_dst.height,
*src_dst.bounds,
)
mercator_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
# Correction factor for web-mercator projection latitude scale change
latitude_correction_factor = math.cos(math.radians(lat))
adjusted_resolution = mercator_resolution * latitude_correction_factor
max_zoom = zoom_for_pixelsize(adjusted_resolution, tilesize=tilesize)
ovr_resolution = adjusted_resolution * max(h, w) / tilesize
min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)
return (min_zoom, max_zoom)
示例9: get_overview_level
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_overview_level(
src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
bounds: Tuple[float, float, float, float],
height: int,
width: int,
dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> int:
"""
Return the overview level corresponding to the tile resolution.
Freely adapted from https://github.com/OSGeo/gdal/blob/41993f127e6e1669fbd9e944744b7c9b2bd6c400/gdal/apps/gdalwarp_lib.cpp#L2293-L2362
Attributes
----------
src_dst : rasterio.io.DatasetReader
Rasterio io.DatasetReader object
bounds : list
Bounds (left, bottom, right, top) in target crs ("dst_crs").
height : int
Output height.
width : int
Output width.
dst_crs: CRS or str, optional
Target coordinate reference system (default "epsg:3857").
Returns
-------
ovr_idx: Int or None
Overview level
"""
dst_transform, _, _ = calculate_default_transform(
src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
)
src_res = dst_transform.a
# Compute what the "natural" output resolution
# (in pixels) would be for this input dataset
vrt_transform = from_bounds(*bounds, width, height)
target_res = vrt_transform.a
ovr_idx = -1
if target_res > src_res:
res = [src_res * decim for decim in src_dst.overviews(1)]
for ovr_idx in range(ovr_idx, len(res) - 1):
ovrRes = src_res if ovr_idx < 0 else res[ovr_idx]
nextRes = res[ovr_idx + 1]
if (ovrRes < target_res) and (nextRes > target_res):
break
if abs(ovrRes - target_res) < 1e-1:
break
else:
ovr_idx = len(res) - 1
return ovr_idx
示例10: get_vrt_transform
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def get_vrt_transform(
src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
bounds: Tuple[float, float, float, float],
height: Optional[int] = None,
width: Optional[int] = None,
dst_crs: CRS = constants.WEB_MERCATOR_CRS,
) -> Tuple[Affine, int, int]:
"""
Calculate VRT transform.
Attributes
----------
src_dst : rasterio.io.DatasetReader
Rasterio io.DatasetReader object
bounds : list
Bounds (left, bottom, right, top) in target crs ("dst_crs").
height : int, optional
Desired output height of the array for the bounds.
width : int, optional
Desired output width of the array for the bounds.
dst_crs: CRS or str, optional
Target coordinate reference system (default "epsg:3857").
Returns
-------
vrt_transform: Affine
Output affine transformation matrix
vrt_width, vrt_height: int
Output dimensions
"""
dst_transform, _, _ = calculate_default_transform(
src_dst.crs, dst_crs, src_dst.width, src_dst.height, *src_dst.bounds
)
w, s, e, n = bounds
if not height or not width:
vrt_width = math.ceil((e - w) / dst_transform.a)
vrt_height = math.ceil((s - n) / dst_transform.e)
vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)
return vrt_transform, vrt_width, vrt_height
tile_transform = from_bounds(w, s, e, n, width, height)
w_res = (
tile_transform.a
if abs(tile_transform.a) < abs(dst_transform.a)
else dst_transform.a
)
h_res = (
tile_transform.e
if abs(tile_transform.e) < abs(dst_transform.e)
else dst_transform.e
)
vrt_width = math.ceil((e - w) / w_res)
vrt_height = math.ceil((s - n) / h_res)
vrt_transform = from_bounds(w, s, e, n, vrt_width, vrt_height)
return vrt_transform, vrt_width, vrt_height
示例11: _reproject
# 需要导入模块: from rasterio import warp [as 别名]
# 或者: from rasterio.warp import calculate_default_transform [as 别名]
def _reproject(input_data, input_type, input_crs, target_crs, dest_path,
resampling_method='bicubic'):
input_crs = _check_crs(input_crs)
target_crs = _check_crs(target_crs)
if input_type == 'vector':
output = input_data.to_crs(target_crs)
if dest_path is not None:
output.to_file(dest_path, driver='GeoJSON')
elif input_type == 'raster':
if isinstance(input_data, rasterio.DatasetReader):
transform, width, height = calculate_default_transform(
input_crs.to_wkt("WKT1_GDAL"), target_crs.to_wkt("WKT1_GDAL"),
input_data.width, input_data.height, *input_data.bounds
)
kwargs = input_data.meta.copy()
kwargs.update({'crs': target_crs.to_wkt("WKT1_GDAL"),
'transform': transform,
'width': width,
'height': height})
if dest_path is not None:
with rasterio.open(dest_path, 'w', **kwargs) as dst:
for band_idx in range(1, input_data.count + 1):
rasterio.warp.reproject(
source=rasterio.band(input_data, band_idx),
destination=rasterio.band(dst, band_idx),
src_transform=input_data.transform,
src_crs=input_data.crs,
dst_transform=transform,
dst_crs=target_crs.to_wkt("WKT1_GDAL"),
resampling=getattr(Resampling, resampling_method)
)
output = rasterio.open(dest_path)
input_data.close()
else:
output = np.zeros(shape=(height, width, input_data.count))
for band_idx in range(1, input_data.count + 1):
rasterio.warp.reproject(
source=rasterio.band(input_data, band_idx),
destination=output[:, :, band_idx-1],
src_transform=input_data.transform,
src_crs=input_data.crs,
dst_transform=transform,
dst_crs=target_crs,
resampling=getattr(Resampling, resampling_method)
)
elif isinstance(input_data, gdal.Dataset):
if dest_path is not None:
gdal.Warp(dest_path, input_data,
dstSRS='EPSG:' + str(target_crs.to_epsg()))
output = gdal.Open(dest_path)
else:
raise ValueError('An output path must be provided for '
'reprojecting GDAL datasets.')
return output