Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def _update_warped_vrt_xml(self, x_size, y_size, geo_transform, block_size, working_data_type):
"""Update rasterXsize, rasterYsize, geotransform, block_size and working_data_type"""
node0 = Node.create(str(self.xml))
node0.replaceAttribute('rasterXSize', str(x_size))
node0.replaceAttribute('rasterYSize', str(y_size))
if geo_transform is not None:
invGeotransform = gdal.InvGeoTransform(geo_transform)
# convert proper string style and set to the GeoTransform element
node0.node('GeoTransform').value = str(geo_transform).strip('()')
node0.node('DstGeoTransform').value = str(geo_transform).strip('()')
node0.node('DstInvGeoTransform').value = (
str(invGeotransform[1]).strip('()'))
if node0.node('SrcGeoLocTransformer'):
node0.node('BlockXSize').value = str(x_size)
node0.node('BlockYSize').value = str(y_size)
if block_size is not None:
node0.node('BlockXSize').value = str(block_size)
node0.node('BlockYSize').value = str(block_size)
if working_data_type is not None:
node0.node('WorkingDataType').value = working_data_type
def create_warped_vrt(self, top_left_coord, res, size):
#----------------------------
# generate warp transform
src_geotr = self.src_ds.GetGeoTransform()
src_proj = txt2proj4(self.src_ds.GetProjection())
gcp_proj = None
if not self.options.tps and src_geotr and src_geotr != (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
src_igeotr = gdal.InvGeoTransform(src_geotr)
src_transform = '%s\n%s' % (warp_src_geotr % src_geotr, warp_src_igeotr % src_igeotr)
else:
gcps = self.src_ds.GetGCPs()
assert gcps, 'Neither geotransform, nor gpcs are in the source file %s' % self.src
gcp_lst = [(g.Id, g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY, g.GCPZ) for g in gcps]
ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
if src_proj and gcp_proj != src_proj:
coords = GdalTransformer(
SRC_SRS=gcp_proj,
DST_SRS=src_proj
).transform([g[3:6] for g in gcp_lst])
gcp_lst = [tuple(p[:3] + c) for p, c in zip(gcp_lst, coords)]
ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
if src_proj and gcp_proj != src_proj:
coords = GdalTransformer(SRC_SRS=gcp_proj, DST_SRS=src_proj).transform([g[3:6] for g in gcp_lst])
gcp_lst = [tuple(p[:3]+c) for p, c in zip(gcp_lst, coords)]
gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
#src_transform = warp_src_gcp_transformer % (0, gcp_txt)
src_transform = warp_src_tps_transformer % gcp_txt
res = self.zoom2res(zoom)
#ul_ll, lr_ll = self.coords2longlat([ul_c, lr_c])
ld('max_zoom', zoom, 'size', dst_xsize, dst_ysize, '-tr', res[0], res[1], '-te', ul_c[0], lr_c[1], lr_c[0], ul_c[1], '-t_srs', self.proj_srs)
dst_geotr = ( ul_c[0], res[0], 0.0,
ul_c[1], 0.0, res[1] )
ok, dst_igeotr = gdal.InvGeoTransform(dst_geotr)
assert ok
dst_transform = '%s\n%s' % (warp_dst_geotr % dst_geotr, warp_dst_igeotr % dst_igeotr)
# generate warp options
warp_options = []
def w_option(name, value): # warp options template
return ' <option name="%s">%s</option>' % (name, value)
warp_options.append(w_option('INIT_DEST', 'NO_DATA'))
# generate cut line
if self.options.cut or self.options.cutline:
cut_wkt = self.get_cutline()
else:
cut_wkt = None
if cut_wkt:
def latlon_to_xy(lat, lon, transform):
"""
Convert lat/lon to x/y for raster
NOTE: No projection done - assumes raster has native lat/lon projection
:param lat: latitude
:param lon: longitude
:param transform: GDAL GeoTransform
:return: x, y pair
"""
# Get the reverse direction GeoTransform
_, transform = gdal.InvGeoTransform(transform)
'''
ulx, uly = transform[0], transform[3]
psx, psy = transform[1], transform[5]
x = int(math.floor(ulx + psx * lon))
y = int(math.floor(uly + psy * lat))
'''
# This is to fix to make nearest rounding before casting to integer and apply gdal utility
x, y= gdal.ApplyGeoTransform(transform, lon, lat)
return int(round(x)), int(round(y))
"""Initializes the Tile Information.
Inputs:
tile_path: path the tile file.
"""
self.tile_path = tile_path
ds = gdal.Open(tile_path)
self.width = ds.RasterXSize
self.height = ds.RasterYSize
# Gets the gdal geo transformation tuples
gdal_version = osgeo.gdal.__version__
self._txf = ds.GetGeoTransform()
if gdal_version[0] == '1':
self._inv_txf = gdal.InvGeoTransform(self._txf)[1]
else:
self._inv_txf = gdal.InvGeoTransform(self._txf)
# Gets the transformation from lat/lon to coordinates
wgs84_ref = osr.SpatialReference()
wgs84_ref.ImportFromEPSG(4326) # WGS84
sref = osr.SpatialReference()
sref.ImportFromWkt(ds.GetProjection())
self._transform = osr.CoordinateTransformation(wgs84_ref, sref)
inv_transform = osr.CoordinateTransformation(sref, wgs84_ref)
# Find a loose lat/lon bounding box for quick check without
# having to do full coordinates transformation
corners = []
for x in [0, self.width]:
for y in [0, self.height]:
corners.append([self._txf[0] + self._txf[1] * x + self._txf[2] * y,
self._txf[3] + self._txf[4] * x + self._txf[5] * y])
def LoadTileForLatLng(self, lat, lng):
latf = int(math.ceil(lat))
lngf = int(math.floor(lng))
k = '%s.%s' % (latf, lngf)
if k in self.tile_cache:
return
filename = self.latlng_file[k]
print 'Loading tile %s from %s' % (k, filename)
dataset = gdal.Open(filename)
# Store the inverse geo transform which will map (lat, lng) to array indices
# in this tile.
tx = dataset.GetGeoTransform()
self.txf[k] = gdal.InvGeoTransform(tx)
self.tile_cache[k] = dataset.ReadAsArray().astype(numpy.float)
self.tile_lru[k] = time.clock()
# evict a tile if there are more than tile_lru_size elements in the cache
if len(self.tile_lru) > self.tile_lru_size:
mint = 0
mink = ''
for k in self.tile_lru.keys():
if self.tile_lru[k] > mint:
mint = self.tile_lru[k]
mink = k
print 'Evicting tile %s' % k
del self.tile_cache[k]
del self.tile_lru[k]
).transform([g[3:6] for g in gcp_lst])
gcp_lst = [tuple(p[:3] + c) for p, c in zip(gcp_lst, coords)]
gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
#src_transform = warp_src_gcp_transformer % (0, gcp_txt)
src_transform = warp_src_tps_transformer % gcp_txt
#tl_ll, br_ll = self.coords2longlat([tl_c, br_c])
#~ ld('create_target_dataset', 'max_zoom', self.max_zoom, 'size', size[0], size[1], '-tr', res[0], res[1], '-te', tl_c[0], br_c[1], br_c[0], tl_c[1], '-t_srs', self.proj_srs)
left, top = top_left_coord
dst_geotr = ( left, res[0], 0.0,
top, 0.0, -res[1] )
dst_igeotr = gdal.InvGeoTransform(dst_geotr)
dst_transform = '%s\n%s' % (warp_dst_geotr % dst_geotr, warp_dst_igeotr % dst_igeotr)
# generate warp options
warp_options = []
def w_option(name, value): # warp options template
return ' <option name="%s">%s</option>' % (name, value)
warp_options.append(w_option('INIT_DEST', 'NO_DATA'))
# generate cut line
if self.options.cut or self.options.cutline:
cut_wkt = self.get_cutline()
else:
cut_wkt = None
if cut_wkt:
warp_options.append(w_option('CUTLINE', cut_wkt))
# base zoom level raster size
dst_xsize = lr_pix[0]-ul_pix[0]
dst_ysize = lr_pix[1]-ul_pix[1]
ld('target Upper Left', self.bounds[0], ul_c, self.proj2geog.transform([self.bounds[0], ul_c]))
ld('target Lower Right', self.bounds[1], lr_c, self.proj2geog.transform([self.bounds[1], lr_c]))
# create VRT for base image warp
# generate warp transform
src_geotr = self.src_ds.GetGeoTransform()
src_proj = txt2proj4(self.src_ds.GetProjection())
gcp_proj = None
if not self.options.tps and src_geotr and src_geotr != (0.0, 1.0, 0.0, 0.0, 0.0, 1.0):
ok, src_igeotr = gdal.InvGeoTransform(src_geotr)
assert ok
src_transform = '%s\n%s' % (warp_src_geotr % src_geotr, warp_src_igeotr % src_igeotr)
else:
gcps = self.src_ds.GetGCPs()
assert gcps, 'Neither geotransform, nor gpcs are in the source file %s' % self.src
gcp_lst = [(g.Id, g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY, g.GCPZ) for g in gcps]
ld('src_proj', self.src_ds.GetProjection(), 'gcp_proj', self.src_ds.GetGCPProjection())
gcp_proj = txt2proj4(self.src_ds.GetGCPProjection())
if src_proj and gcp_proj != src_proj:
coords = GdalTransformer(SRC_SRS=gcp_proj, DST_SRS=src_proj).transform([g[3:6] for g in gcp_lst])
gcp_lst = [tuple(p[:3]+c) for p, c in zip(gcp_lst, coords)]
gcp_txt = '\n'.join((gcp_templ % g for g in gcp_lst))
#src_transform = warp_src_gcp_transformer % (0, gcp_txt)
src_transform = warp_src_tps_transformer % gcp_txt