Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
Return
----------
dt : datatype
"""
if gdal_dt == gdal.GDT_Byte:
dt = 'uint8'
elif gdal_dt == gdal.GDT_Int16:
dt = 'int16'
elif gdal_dt == gdal.GDT_UInt16:
dt = 'uint16'
elif gdal_dt == gdal.GDT_Int32:
dt = 'int32'
elif gdal_dt == gdal.GDT_UInt32:
dt = 'uint32'
elif gdal_dt == gdal.GDT_Float32:
dt = 'float32'
elif gdal_dt == gdal.GDT_Float64:
dt = 'float64'
elif gdal_dt == gdal.GDT_CInt16 or gdal_dt == gdal.GDT_CInt32 or gdal_dt == gdal.GDT_CFloat32 or gdal_dt == gdal.GDT_CFloat64:
dt = 'complex64'
else:
print('Data type unkown')
# exit()
return dt
def WriteTifOutput(self,outfile,output,geo, prj, fields):
'''Writes the arrays of an output dictionary which keys match the list in fields to a GeoTIFF '''
import gdal
import numpy as np
rows,cols=np.shape(output['H1'])
driver = gdal.GetDriverByName('GTiff')
nbands=len(fields)
ds = driver.Create(outfile, cols, rows, nbands, gdal.GDT_Float32)
ds.SetGeoTransform(geo)
ds.SetProjection(prj)
for i,field in enumerate(fields):
band=ds.GetRasterBand(i+1)
band.SetNoDataValue(0)
band.WriteArray(output[field])
band.FlushCache()
ds.FlushCache()
del ds
def Save_as_MEM(data='', geo='', projection=''):
"""
This function save the array as a memory file
Keyword arguments:
data -- [array], dataset of the geotiff
geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
pixelsize], (geospatial dataset)
projection -- interger, the EPSG code
"""
# save as a geotiff
driver = gdal.GetDriverByName("MEM")
dst_ds = driver.Create('', int(data.shape[1]), int(data.shape[0]), 1,
gdal.GDT_Float32, ['COMPRESS=LZW'])
srse = osr.SpatialReference()
if projection == '':
srse.SetWellKnownGeogCS("WGS84")
else:
srse.SetWellKnownGeogCS(projection)
dst_ds.SetProjection(srse.ExportToWkt())
dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
dst_ds.SetGeoTransform(geo)
dst_ds.GetRasterBand(1).WriteArray(data)
return(dst_ds)
filename = filename.replace(".vrt", ".tif")
filename = os.path.join(self.output_directory, filename)
metadata = get_dataset_metadata(arg)
data = get_dataset_data_with_pq(arg, Ls57Arg25Bands, pqa)
# Calculate TCI Wetness
tci = calculate_tassel_cap_index(data,
coefficients=TCI_COEFFICIENTS[arg.satellite][TasselCapIndex.WETNESS])
_log.info("TCI shape is %s | min = %s | max = %s", numpy.shape(tci), tci.min(), tci.max())
raster_create(filename,
[tci],
metadata.transform, metadata.projection, numpy.nan, gdal.GDT_Float32)
def write_tif(path_tif, array, geotransform, geoprojection, size):
dim_array = array.shape
if len(dim_array) > 2:
depth = dim_array[2]
else:
depth = 1
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(
path_tif, size[0], size[1], depth, gdal.GDT_Float32)
# sets same geotransform as input
outdata.SetGeoTransform(geotransform)
outdata.SetProjection(geoprojection) # sets same projection as input
for i in range(depth):
try:
arr = array[:, :, i]
except:
arr = array[:, :]
arr = cv2.resize(arr, size)
outdata.GetRasterBand(i+1).WriteArray(arr)
# outdata.GetRasterBand(1).SetNoDataValue(-9999)##if you want these values transparent
outdata.FlushCache() # saves to disk!!
products = dc.list_products()
resolution = products.resolution[products.name == 'ls7_ledaps']
lon_dist = resolution.values[0][1]
lat_dist = resolution.values[0][0]
# Rotation
lon_rtn = 0
lat_rtn = 0
geotransform = (ul_lon, lon_dist, lon_rtn, ul_lat, lat_rtn, lat_dist)
dataset_out = frac_coverage_classify(dataset_in)
out_file = (str(min_lon) + '_' + str(min_lat) + '_' + start_date_str + '_' + end_date_str + '_frac_coverage.tif')
utilities.save_to_geotiff(out_file, gdal.GDT_Float32, dataset_out, geotransform, spatial_ref)
"""
gt = [origin[0], pixel_width, 0, origin[1], 0, pixel_height]
# Apply rotation by tweaking geotransform. The data remains the
# same but will appear roated in a viewer e.g. ArcGIS.
if angle:
if rotate_origin:
gt = rotate_transform(gt, angle, origin[0], origin[1])
else:
gt = rotate_transform(gt, angle, center.east, center.north)
filename = '{}_rot{}.tif'\
.format(os.path.splitext(filename)[0], angle)
rows, cols = data.shape
driver = gdal.GetDriverByName('GTiff')
out_raster = driver.Create(filename, cols, rows, 1, gdal.GDT_Float32)
out_raster.SetGeoTransform(gt)
out_band = out_raster.GetRasterBand(1)
out_band.SetNoDataValue(ndv)
out_band.WriteArray(data)
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg_code)
out_raster.SetProjection(srs.ExportToWkt())
out_band.FlushCache()
# output to ascii format
ascii_filename = "{}.asc".format(os.path.splitext(filename)[0])
driver2 = gdal.GetDriverByName('AAIGrid')
driver2.CreateCopy(ascii_filename, out_raster)
return filename, ascii_filename
# Start the gdal driver for GeoTIFF
if outPath == "MEM":
driver = gdal.GetDriverByName("MEM")
driverOpt = []
shape = data.shape
if len(shape) > 2:
ds = driver.Create(outPath, shape[1], shape[0], shape[2], gdal.GDT_Float32, driverOpt)
ds.SetProjection(proj)
ds.SetGeoTransform(geotransform)
for i in range(shape[2]):
ds.GetRasterBand(i+1).WriteArray(data[:, :, i])
ds.GetRasterBand(i+1).SetNoDataValue(noDataValue)
else:
ds = driver.Create(outPath, shape[1], shape[0], 1, gdal.GDT_Float32, driverOpt)
ds.SetProjection(proj)
ds.SetGeoTransform(geotransform)
ds.GetRasterBand(1).WriteArray(data)
ds.GetRasterBand(1).SetNoDataValue(noDataValue)
return ds
def create_output_file(base_filepath, out_filepath, raster_count = 1, dataType = gdal.GDT_Float32, \
imageFormat = 'GTiff', formatOptions = ['COMPRESS=LZW', 'TILED=True', 'BIGTIFF=YES']):
driver = gdal.GetDriverByName(imageFormat)
base_ds = gdal.Open(base_filepath)
x_start, pixel_width, _, y_start, _, pixel_height = base_ds.GetGeoTransform()
x_size = base_ds.RasterXSize
y_size = base_ds.RasterYSize
out_srs = osr.SpatialReference()
out_srs.ImportFromWkt(base_ds.GetProjectionRef())
output_img_ds = driver.Create(out_filepath, x_size, y_size, raster_count, dataType, formatOptions)
output_img_ds.SetGeoTransform((x_start, pixel_width, 0, y_start, 0, pixel_height))
output_img_ds.SetProjection(out_srs.ExportToWkt())
def numpy_to_layer(self,array,name):
array = np.rot90(array)
sy, sx = array.shape
pathname = QgsProject.instance().readPath("./")+'/'+name
driver = gdal.GetDriverByName("GTiff")
dsOut = driver.Create(pathname, sx+1,sy+1,1,gdal.GDT_Float32 ,)
dsOut.SetGeoTransform(self.transform)
dsOut.SetProjection(self.wkt)
bandOut=dsOut.GetRasterBand(1)
gdalnumeric.BandWriteArray(bandOut, array)
bandOut = None
dsOut = None
layer = QgsRasterLayer(pathname,name)
QgsProject.instance().addMapLayer(layer)
def run_calculator(self,string,name):