Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def sieve_strata(config, image):
""" First create a band in memory that's that's just 1s and 0s """
src_ds = gdal.Open( image, gdal.GA_ReadOnly )
srcband = src_ds.GetRasterBand(1)
srcarray = srcband.ReadAsArray()
mem_rast = save_raster_memory(srcarray, image)
mem_band = mem_rast.GetRasterBand(1)
#Now the code behind gdal_sieve.py
maskband = None
drv = gdal.GetDriverByName('GTiff')
#Need to output sieved file
dst_path = config['postprocessing']['sieve']['sieve_file']
dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_sieve_mask.tif'
dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
srcband.DataType )
wkt = src_ds.GetProjection()
if wkt != '':
dst_ds.SetProjection( wkt )
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dstband = dst_ds.GetRasterBand(1)
# Parameters
prog_func = None
def do_proximity(config, image, srcarray, label, _class):
""" Get proximity of each pixel to to certain value """
#First create a band in memory that's that's just 1s and 0s
src_ds = gdal.Open( image, gdal.GA_ReadOnly )
srcband = src_ds.GetRasterBand(1)
mem_rast = save_raster_memory(srcarray, image)
mem_band = mem_rast.GetRasterBand(1)
#Now the code behind gdal_sieve.py
maskband = None
drv = gdal.GetDriverByName('GTiff')
#Need to output sieved file
dst_path = config['postprocessing']['deg_class']['prox_dir']
dst_filename = dst_path + '/' + image.split('/')[-1].split('.')[0] + '_' + _class + '.tif'
dst_ds = drv.Create( dst_filename,src_ds.RasterXSize, src_ds.RasterYSize,1,
srcband.DataType )
wkt = src_ds.GetProjection()
if wkt != '':
dst_ds.SetProjection( wkt )
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dstband = dst_ds.GetRasterBand(1)
# Parameters
import unittest
import pytest
import glob
import copy
import pyrate.core.config as cf
from pyrate import conv2tif, prepifg
from pyrate.core import gdal_python
import os
import gdal
import numpy as np
import osr
NO_DATA_VALUE = 0
driver = gdal.GetDriverByName( 'GTiff' )
# sample gdal dataset
sample_gdal_filename = 'sample_gdal_dataset.tif'
sample_gdal_dataset = driver.Create(sample_gdal_filename, 5, 5, 1, gdal.GDT_Float32)
srs = osr.SpatialReference()
wkt_projection = srs.ExportToWkt()
sample_gdal_dataset.SetProjection(wkt_projection)
sample_gdal_band = sample_gdal_dataset.GetRasterBand(1)
sample_gdal_band.SetNoDataValue(NO_DATA_VALUE)
sample_gdal_band.WriteArray(np.arange(25).reshape(5, 5))
print("sample_gdal_band: \n", sample_gdal_band.ReadAsArray())
# coherence mask dataset
coherence_mask_filename = 'coherence_mask_dataset.tif'
def convert_raster(input_raster, output_raster, output_format = "GTiff"):
src_ds = gdal.Open(input_raster) #Open existing dataset
driver = gdal.GetDriverByName(output_format) #Open output format driver, see gdal_translate --formats for list
dst_ds = driver.CreateCopy(output_raster, src_ds, 0) #Output to new format
dst_ds = None #Properly close the datasets to flush to disk
src_ds = None
for i in range(layerDefinition.GetFieldCount()):
print(layerDefinition.GetFieldDefn(i).GetName())
# The raster file to be created and receive the rasterized shapefile
outrastername = shapefileshortname + '.tif'
outraster = shapefilefilepath+os.sep+ outrastername
outcsv = shapefilefilepath+os.sep+shapefileshortname+'_lithokey.csv'
print("Full name of out raster is: "+outraster)
# Create the destination data source
inGridSize=float(raster_resolution)
xMin, xMax, yMin, yMax = daLayer.GetExtent()
xRes = int((xMax - xMin) / inGridSize)
yRes = int((yMax - yMin) / inGridSize)
rasterDS = gdal.GetDriverByName('GTiff').Create(outraster, xRes, yRes, 1, gdal.GDT_Byte)
# Define spatial reference
NoDataVal = -9999
rasterDS.SetProjection(daLayer.GetSpatialRef().ExportToWkt())
rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize))
rBand = rasterDS.GetRasterBand(1)
rBand.SetNoDataValue(NoDataVal)
rBand.Fill(NoDataVal)
# Rasterize
gdal.RasterizeLayer(rasterDS, [1], daLayer, options = ["ATTRIBUTE=GEOL_CODE"])
# Make a key for the bedrock
geol_dict = dict()
for feature in daLayer:
ID = feature.GetField(geol_field)
if ct is not None:
for i in range(min(256,ct.GetCount())):
entry = ct.GetColorEntry(i)
for c in range(4):
lookup[c][i] = entry[c]
# ----------------------------------------------------------------------------
# Create the working file.
if format == 'GTiff':
tif_filename = dst_filename
else:
tif_filename = 'temp.tif'
gtiff_driver = gdal.GetDriverByName( 'GTiff' )
tif_ds = gtiff_driver.Create( tif_filename,
src_ds.RasterXSize, src_ds.RasterYSize, out_bands )
# ----------------------------------------------------------------------------
# We should copy projection information and so forth at this point.
tif_ds.SetProjection( src_ds.GetProjection() )
tif_ds.SetGeoTransform( src_ds.GetGeoTransform() )
# ----------------------------------------------------------------------------
# Do the processing one scanline at a time.
gdal.TermProgress( 0.0 )
for iY in range(src_ds.RasterYSize):
def output_ds(out_array, img_params, d_type=GDT_Unknown, fn='result.tif'):
"""
Helper function.
Writes new data-set into disk
and saves output arrays in the data-set.
"""
# create output raster data-set
cols = img_params[0]
rows = img_params[1]
bands = 1 # ndvi image needs only one band
gt = img_params[3]
proj = img_params[4]
driver = gdal.GetDriverByName('GTiff')
driver.Register()
out_ras = driver.Create(fn, cols, rows, bands, d_type)
out_ras.SetGeoTransform(gt)
out_ras.SetProjection(proj)
out_band = out_ras.GetRasterBand(1)
out_band.WriteArray(out_array)
out_band.SetNoDataValue(-99)
out_band.FlushCache()
out_band.GetStatistics(0, 1)
return
# if the tile not exists or cannot be opened, create a nan array with the right projection
except:
if Horizontal==TilesHorizontal[0] and Vertical==TilesVertical[0]:
x1 = (TilesHorizontal[0] - 19) * 4800 * Distance
x4 = (TilesVertical[0] - 9) * 4800 * -1 * Distance
geo = [x1, Distance, 0.0, x4, 0.0, -Distance]
geo_t=tuple(geo)
proj='PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
data=np.ones((4800,4800)) * (-9999)
countYdata=(TilesVertical[1] - TilesVertical[0] + 2) - countY
DataTot[(countYdata - 1) * 4800:countYdata * 4800,(countX - 1) * 4800:countX * 4800] = data
# Make geotiff file
name2 = os.path.join(output_folder, 'Merged.tif')
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(name2, DataTot.shape[1], DataTot.shape[0], 1, gdal.GDT_Float32, ['COMPRESS=LZW'])
try:
dst_ds.SetProjection(proj)
except:
proj='PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not specified (based on custom spheroid)",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
x1 = (TilesHorizontal[0] - 18) * 4800 * Distance
x4 = (TilesVertical[0] - 9) * 4800 * -1 * Distance
geo = [x1, Distance, 0.0, x4, 0.0, -Distance]
geo_t = tuple(geo)
dst_ds.SetProjection(proj)
dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
dst_ds.SetGeoTransform(geo_t)
dst_ds.GetRasterBand(1).WriteArray(DataTot*0.0001)
dst_ds = None
sds = None
# End logging, print blank line for clarity
print('')
bands = src_image.RasterCount
if bands > 1:
for i in range(bands):
outBand = output_dataset.GetRasterBand(i + 1)
outBand.SetNoDataValue(nodata_values[i])
outBand.WriteArray(clip[i])
else:
outBand = output_dataset.GetRasterBand(1)
outBand.SetNoDataValue(nodata_values[0])
outBand.WriteArray(clip)
if dst_img_file:
output_driver = gdal.GetDriverByName('GTiff')
output_driver.CreateCopy(
dst_img_file, output_dataset, False)
"""
self = Interferogram
suffi = {'GTiff':'tif',
'ENVI':'bin'}
try:
suffix = suffi[format]
except:
suffix = ''
if nparray == None:
nparray = load_half(self, half=half)
if outfile == None:
outfile = self.Name + '.' + suffix
try:
driver = gdal.GetDriverByName(format)
nd = driver.Register()
except:
print("unrecognized format, check 'gdalinfo --formats'")
# DATA
rows,cols = nparray.shape
#note order: cols, rows, nbands, dtype (default GDT_Byte=uint8)
outDataset = driver.Create(outfile, cols, rows, 1, gdal.GDT_Float32)
outBand = outDataset.GetRasterBand(1)
if nanval:
print('setting nans to {}'.format(nanval))
outBand.SetNoDataValue(nanval) #or np.NaN?
nparray[np.isnan(nparray)] = nanval
#order: array, xoffset, yoffset (integers)... if successful result=0
result = outBand.WriteArray(nparray, 0, 0)