Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def init_from_name(self, filename):
"""
Initialize file_info from filename
filename -- Name of file to read.
Returns 1 on success or 0 if the file can't be opened.
"""
namebbox = filename.split('::')
bbox = None
filename = namebbox[0]
if len(namebbox) > 1:
bbox = map(float, namebbox[1].split(':'))
fh = gdal.Open( filename, gdal.GA_ReadOnly )
if fh is None:
return False
self.filename = filename
self.bands = fh.RasterCount
self.xsize = fh.RasterXSize
self.ysize = fh.RasterYSize
self.projection = fh.GetProjection()
self.gcpprojection = fh.GetGCPProjection()
self.gcps = fh.GetGCPs()
self.geotransform = fh.GetGeoTransform()
if bbox:
self.geotransform = [0.0,0.0,0.0,0.0,0.0,0.0]
self.geotransform[0] = bbox[0]
self.geotransform[1] = (bbox[2] - bbox[0]) / float(self.xsize)
def Tiff2Envi(infile):
#check if shapefile exists
if infile==None:
return
#Get Path and Name of Inputfile
(infilefilepath, infilename) = os.path.split(infile) #get path and filename seperately
(infileshortname, extension) = os.path.splitext(infilename) #get file name without extension
outfile = infilepath + '\\'+ infileshortname + '.img' #create outputfilename
#Open Dataset and Read as Array
dataset = gdal.Open(infile, gdal.GA_ReadOnly)
band = dataset.GetRasterBand(1)
data = band.ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize)
#Driver for output dataset
driver = gdal.GetDriverByName('ENVI')
driver.Register()
#Get Dimensions for Outputdataset
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
datatype = band.DataType
#Create outputfile
print 'Creating ', outfile, ' with ', cols, rows, bands, datatype
outDataset = driver.Create(outfile, cols, rows, bands, datatype)
def lowPassFilter(self,dataFile, sigDataFile, maskFile, Sx, Sy, sig_x, sig_y, iteration=5, theta=0.0):
ds = gdal.Open(dataFile + '.vrt', gdal.GA_ReadOnly)
length = ds.RasterYSize
width = ds.RasterXSize
dataIn = np.memmap(dataFile, dtype=np.float32, mode='r', shape=(length,width))
sigData = np.memmap(sigDataFile, dtype=np.float32, mode='r', shape=(length,width))
mask = np.memmap(maskFile, dtype=np.byte, mode='r', shape=(length,width))
dataF, sig_dataF = iterativeFilter(self,dataIn[:,:], mask[:,:], sigData[:,:], iteration, Sx, Sy, sig_x, sig_y, theta)
filtDataFile = dataFile + ".filt"
sigFiltDataFile = sigDataFile + ".filt"
filtData = np.memmap(filtDataFile, dtype=np.float32, mode='w+', shape=(length,width))
filtData[:,:] = dataF[:,:]
filtData.flush()
sigFilt= np.memmap(sigFiltDataFile, dtype=np.float32, mode='w+', shape=(length,width))
tilebands)
# TODO: Implement more clever walking on the tiles with cache functionality
# probably walk should start with reading of four tiles from top left corner
# Hilbert curve...
children = []
# Read the tiles and write them to query window
for y in range(2 * ty, 2 * ty + 2):
for x in range(2 * tx, 2 * tx + 2):
minx, miny, maxx, maxy = self.tminmax[tz + 1]
if x >= minx and x <= maxx and y >= miny and y <= maxy:
dsquerytile = gdal.Open(
path.join(self.output, str(tz + 1), str(x),
"%s.%s" %
(y, self.tileext)), gdal.GA_ReadOnly)
if (ty == 0 and y == 1) or (ty != 0 and
(y % (2 * ty)) != 0):
tileposy = 0
else:
tileposy = self.tilesize
if tx:
tileposx = x % (2 * tx) * self.tilesize
elif tx == 0 and x == 1:
tileposx = self.tilesize
else:
tileposx = 0
dsquery.WriteRaster(
tileposx,
tileposy,
self.tilesize,
self.tilesize,
num_full_path = requested_files.at[i + 1, 'Full Path']
num_filename = requested_files.at[i + 1, 'Full Path']
zipped_num_data = ZipFile(num_full_path)
zipped_num_full_path = zipped_num_data.infolist()[0].filename
num_data = np.frombuffer(zipped_num_data.open(zipped_num_full_path).read(),
np.dtype('uint8')).reshape(array_shape)
masked_dem_data[(num_data == 1) | (num_data == 2)] = water_value
i += 1
zipped_hgt_data = ZipFile(hgt_full_path)
dem_dataset = gdal.Open(hgt_full_path, gdal.GA_ReadOnly)
dem_data = dem_dataset.ReadAsArray()
masked_dem_data *= dem_data
metadata_dict[label]['WKT'] = dem_dataset.GetProjection()
metadata_dict[label]['GeoTransform'] = dem_dataset.GetGeoTransform()
else:
geo_transform = []
geo_transform.append(extents[0])
geo_transform.append(x_res)
geo_transform.append(0)
geo_transform.append(extents[-1])
def main():
print("For: {}, extract elev data in: {}".format(sys.argv[1], sys.argv[2]))
extr_data = open(sys.argv[2], "w")
#local_max = open(sys.argv[3], "w")
dataset = gdal.Open(sys.argv[1], gdal.GA_ReadOnly)
if not dataset:
print("ERROR: CANNOT OPEN DEM FILE")
elif not extr_data:
print("ERROR: CANNOT OPEN EXTRACT DATA OUTPUT FILE")
else:
# Contains metadata and pointer to channels (i.e. rgb) of DEM
geotransform = dataset.GetGeoTransform()
if geotransform:
num_cols = dataset.RasterXSize
num_rows = dataset.RasterYSize
# Get coordinates in cartesian
ext = get_coords.get_corner_coordinates(
geotransform, num_cols, num_rows
print("Error reading the RPC")
return ERROR
# open the DSM
dsm = gdal.Open(args_dsm, gdal.GA_ReadOnly)
if not dsm:
return ERROR
band = dsm.GetRasterBand(1)
dsmRaster = band.ReadAsArray(
xoff=0, yoff=0,
win_xsize=dsm.RasterXSize, win_ysize=dsm.RasterYSize)
dsm_nodata_value = band.GetNoDataValue()
print("DSM raster shape {}".format(dsmRaster.shape))
if args_dtm:
dtm = gdal.Open(args_dtm, gdal.GA_ReadOnly)
if not dtm:
return ERROR
band = dtm.GetRasterBand(1)
dtmRaster = band.ReadAsArray(
xoff=0, yoff=0,
win_xsize=dtm.RasterXSize, win_ysize=dtm.RasterYSize)
newRaster = numpy.where(dsmRaster != dsm_nodata_value, dsmRaster, dtmRaster)
dsmRaster = newRaster
# apply morphology to denoise the DSM
if (args_denoise_radius > 0):
morph_struct = circ_structure(args_denoise_radius)
dsmRaster = morphology.grey_opening(dsmRaster, structure=morph_struct)
dsmRaster = morphology.grey_closing(dsmRaster, structure=morph_struct)
# create the rectified image
def transfer_metadata(corrected_image_path, original_image_path):
corrected_dataset = gdal.Open(corrected_image_path, gdal.GA_Update)
original_dataset = gdal.Open(original_image_path, gdal.GA_ReadOnly)
corrected_dataset.SetMetadata(original_dataset.GetMetadata())
rpcs = original_dataset.GetMetadata('RPC')
corrected_dataset.SetMetadata(rpcs, 'RPC')
corrected_dataset.SetGeoTransform(original_dataset.GetGeoTransform())
corrected_dataset.SetProjection(original_dataset.GetProjection())
corrected_dataset = None
original_dataset = None
def getSizes(im):
raster = gdal.Open(im,gdal.GA_ReadOnly)
x = raster.RasterXSize
y = raster.RasterYSize
d = raster.RasterCount
return raster,x,y,d
data = np.array(ds.ReadAsArray(), dtype=np.float32)
data[data == ds.GetRasterBand(1).GetNoDataValue()] = np.nan
h5['height'][:,:] = data
# slantRangeDistance
h5['slantRangeDistance'][:,:] = float(h5.attrs['STARTING_RANGE'])
# incidenceAngle
ds = gdal.Open(incAngleFile, gdal.GA_ReadOnly)
data = ds.ReadAsArray()
data[data == ds.GetRasterBand(1).GetNoDataValue()] = np.nan
h5['incidenceAngle'][:,:] = data
# azimuthAngle
if azAngleFile is not None:
ds = gdal.Open(azAngleFile, gdal.GA_ReadOnly)
data = ds.ReadAsArray()
data[data == ds.GetRasterBand(1).GetNoDataValue()] = np.nan
# azimuth angle of the line-of-sight vector:
# ARIA: vector from target to sensor measured from the east in counterclockwise direction
# ISCE: vector from sensor to target measured from the north in counterclockwise direction
# convert ARIA format to ISCE format, which is used in mintpy
data -= 90
h5['azimuthAngle'][:,:] = data
# waterMask
if waterMaskFile is not None:
ds = gdal.Open(waterMaskFile, gdal.GA_ReadOnly)
water_mask = ds.ReadAsArray()
water_mask[water_mask == ds.GetRasterBand(1).GetNoDataValue()] = False
# assign False to invalid pixels based on incAngle data