Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"Either gdal2tiles with parameter -p 'raster' or use another GIS software for georeference e.g. gdal_transform -gcp / -a_ullr / -a_srs")
in_srs_code = self.in_srs.GetAttrValue("AUTHORITY", 0)
in_ds_srs = osr.SpatialReference()
res = in_ds_srs.ImportFromWkt(in_ds.GetProjection())
if res != 0 and in_srs_code is None:
print "ERROR! The inumpyut file %s has no SRS associated and no SRS has been defined in inumpyut (-s parameter)" % _inumpyut
exit(1)
if self.in_srs:
if in_ds_srs.ExportToProj4() != self.out_srs.ExportToProj4():
if (self.in_srs.ExportToProj4() != self.out_srs.ExportToProj4()) or (in_ds.GetGCPCount() != 0):
print "WARNING! Inumpyut file %s has a SR different from EPSG:4326 (WGS84). This can make the processing significantly slow." % _inumpyut
# Generation of VRT dataset in tile projection, default 'nearest neighbour' warping
out_ds = gdal.AutoCreateWarpedVRT( in_ds, self.in_srs_wkt, self.out_srs.ExportToWkt() )
# TODO: HIGH PRIORITY: Correction of AutoCreateWarpedVRT according the max zoomlevel for correct direct warping!!!
if self.options.verbose:
print("Warping of the raster by AutoCreateWarpedVRT (result saved into 'tiles.vrt')")
out_ds.GetDriver().CreateCopy("%s.vrt" % _inumpyut, out_ds)
inumpyut_or_vrt = "%s.vrt" % _inumpyut
# Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!!
else:
self.error("Inumpyut file has unknown SRS.", "Use --s_srs ESPG:xyz (or similar) to provide source reference system." )
if out_ds and self.options.verbose:
print("Projected file:", "tiles.vrt", "( %sP x %sL - %s bands)" % (out_ds.RasterXSize, out_ds.RasterYSize, out_ds.RasterCount))
def calc_zoom(self, zoom_parm):
'determine and set a list of zoom levels to generate'
#----------------------------
# check raster parameters to find default zoom range
ld('automatic zoom levels')
# modify target srs to allow charts crossing meridian 180
shifted_srs = self.shift_srs()
t_ds = gdal.AutoCreateWarpedVRT(self.src_ds, None, txt2wkt(shifted_srs))
geotr = t_ds.GetGeoTransform()
res = (geotr[1], geotr[5])
max_zoom = max(self.res2zoom_xy(res))
# calculate min_zoom
ul_c = (geotr[0], geotr[3])
lr_c = gdal.ApplyGeoTransform(geotr, t_ds.RasterXSize, t_ds.RasterYSize)
wh = (lr_c[0]-ul_c[0], ul_c[1]-lr_c[1])
ld('ul_c, lr_c, wh', ul_c, lr_c, wh)
min_zoom = min(self.res2zoom_xy([wh[i]/abs(self.tile_dim[i]) for i in (0, 1)]))
self.set_zoom_range(zoom_parm, (min_zoom, max_zoom))
"""
# VRT to be warped
src_vrt = self.copy()
# if destination GCPs are given: create and add fake GCPs to src
dst_wkt = src_vrt._set_fake_gcps(dst_srs, dst_gcps, skip_gcps)
if resize_only:
src_vrt._set_geotransform_for_resize()
else:
src_vrt._set_gcps_geolocation_geotransform()
# create Warped VRT GDAL Dataset
warped_dataset = gdal.AutoCreateWarpedVRT(src_vrt.dataset, None, dst_wkt, resample_alg)
# check if Warped VRT was created
if warped_dataset is None:
raise NansatGDALError('Cannot create warpedVRT!')
# create VRT object from Warped VRT GDAL Dataset
warped_vrt = VRT.copy_dataset(warped_dataset)
# set x/y size, geo_transform, block_size
warped_vrt._update_warped_vrt_xml(x_size, y_size, geo_transform, block_size, working_data_type)
# apply thin-spline-transformation option
if self.tps:
warped_vrt.write_xml(warped_vrt.xml.replace('GCPTransformer', 'TPSTransformer'))
# if given, add dst GCPs
def calc_zoom(self, zoom_parm):
'determine and set a list of zoom levels to generate'
#----------------------------
# check raster parameters to find default zoom range
ld('automatic zoom levels')
# modify target srs to allow charts crossing meridian 180
shifted_srs = self.shift_srs()
t_ds = gdal.AutoCreateWarpedVRT(self.src_ds, None, txt2wkt(shifted_srs))
geotr = t_ds.GetGeoTransform()
res = (geotr[1], geotr[5])
max_zoom = max(self.res2zoom_xy(res))
# calculate min_zoom
ul_c = (geotr[0], geotr[3])
lr_c = gdal.ApplyGeoTransform(geotr, t_ds.RasterXSize, t_ds.RasterYSize)
wh = (lr_c[0]-ul_c[0], ul_c[1]-lr_c[1])
ld('ul_c, lr_c, wh', ul_c, lr_c, wh)
min_zoom = min(self.res2zoom_xy([wh[i]/abs(self.tile_dim[i]) for i in (0, 1)]))
self.set_zoom_range(zoom_parm, (min_zoom, max_zoom))
1.0,
) and self.in_ds.GetGCPCount() == 0:
self.error("There is no georeference - neither affine transformation (worldfile) nor GCPs. You can generate only 'raster' profile tiles."
,
"Either gdal2tiles with parameter -p 'raster' or use another GIS software for georeference e.g. gdal_transform -gcp / -a_ullr / -a_srs"
)
if self.in_srs:
if self.in_srs.ExportToProj4() \
!= self.out_srs.ExportToProj4() \
or self.in_ds.GetGCPCount() != 0:
# Generation of VRT dataset in tile projection, default 'nearest neighbour' warping
self.out_ds = gdal.AutoCreateWarpedVRT(self.in_ds,
self.in_srs_wkt, self.out_srs.ExportToWkt())
# TODO: HIGH PRIORITY: Correction of AutoCreateWarpedVRT according the max zoomlevel for correct direct warping!!!
if self.options.verbose:
print "Warping of the raster by AutoCreateWarpedVRT (result saved into 'tiles.vrt')"
self.out_ds.GetDriver().CreateCopy('tiles.vrt',
self.out_ds)
# Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!!
# Correction of AutoCreateWarpedVRT for NODATA values
if self.in_nodata != []:
import tempfile
tempfilename = tempfile.mktemp('-gdal2tiles.vrt'
# calculate zoom range
self.calc_zoom(zoom_parm)
self.max_zoom = self.zoom_range[0]
# shift target SRS to avoid crossing 180 meridian
shifted_srs = self.shift_srs(self.max_zoom)
shift_x = GdalTransformer(SRC_SRS=shifted_srs, DST_SRS=self.proj_srs).transform_point((0, 0))[0]
if shift_x != 0:
self.proj_srs = shifted_srs
self.proj2geog = GdalTransformer(SRC_SRS=self.proj_srs, DST_SRS=self.geog_srs)
self.pix_origin = (self.pix_origin[0]-shift_x, self.pix_origin[1])
self.tile_origin = (self.tile_origin[0]-shift_x, self.tile_origin[1])
ld('new_srs', shifted_srs, 'shift_x', shift_x, 'pix_origin', self.pix_origin)
# get corners at the target SRS
target_ds = gdal.AutoCreateWarpedVRT(self.src_ds, None, txt2wkt(shifted_srs))
target_bounds = GdalTransformer(target_ds).transform([
(0, 0),
(target_ds.RasterXSize, target_ds.RasterYSize)])
# self.bounds are set to a world raster, now clip to the max tileset area
self.bounds = ((target_bounds[0][0],
min(self.bounds[0][1], target_bounds[0][1])),
(target_bounds[1][0],
max(self.bounds[1][1], target_bounds[1][1])))
ld('target raster')
ld('Upper Left', self.bounds[0], target_bounds[0], self.proj2geog.transform([self.bounds[0], target_bounds[0]]))
ld('Lower Right', self.bounds[1], target_bounds[1], self.proj2geog.transform([self.bounds[1], target_bounds[1]]))
# orig_ul = GdalTransformer(SRC_SRS=self.geog_srs, DST_SRS=self.srs).transform_point(
# self.proj2geog.transform_point(target_bounds[0]))
# ld(orig_ul[0]-target_bounds[0][0], orig_ul)
if (self.in_ds.GetGeoTransform() ==
(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)) and (
self.in_ds.GetGCPCount() == 0):
self.error(
"There is no georeference - neither affine transformation (worldfile) nor GCPs. You can generate only 'raster' profile tiles.",
"Either gdal2tiles with parameter -p 'raster' or use another GIS software for georeference e.g. gdal_transform -gcp / -a_ullr / -a_srs")
if self.in_srs:
if (self.in_srs.ExportToProj4() !=
self.out_srs.ExportToProj4()) or (
self.in_ds.GetGCPCount() != 0):
# Generation of VRT dataset in tile projection, default 'nearest neighbour' warping
self.out_ds = gdal.AutoCreateWarpedVRT(
self.in_ds, self.in_srs_wkt,
self.out_srs.ExportToWkt())
# TODO: HIGH PRIORITY: Correction of AutoCreateWarpedVRT according the max zoomlevel for correct direct warping!!!
if self.options.verbose:
print(
"Warping of the raster by AutoCreateWarpedVRT (result saved into 'tiles.vrt')")
self.out_ds.GetDriver().CreateCopy("tiles.vrt",
self.out_ds)
# Note: self.in_srs and self.in_srs_wkt contain still the non-warped reference system!!!
# Correction of AutoCreateWarpedVRT for NODATA values
if self.in_nodata != []:
tempfilename = mktemp('-gdal2tiles.vrt')
"""
# Open source dataset
src_ds = get_dataset(src)
# Define target SRS
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(int(epsg))
dst_wkt = dst_srs.ExportToWkt()
# Resampling might be passed as a string
if not isinstance(resampling, int):
resampling = getattr(gdal, resampling)
# Call AutoCreateWarpedVRT() to fetch default values
# for target raster dimensions and geotransform
reprojected_ds = gdal.AutoCreateWarpedVRT(src_ds,
None,
dst_wkt,
resampling,
error_threshold)
# Create the final warped raster
if dst:
gdal.GetDriverByName('GTiff').CreateCopy(dst, reprojected_ds)
return reprojected_ds