Python osgeo.gdal.GDT_Byte() Examples

The following are code examples for showing how to use osgeo.gdal.GDT_Byte(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 votes vote down vote up
def writeRaster(arr, out_path, template):
    no_data = 0
    # First of all, gather some information from the template file
    data = gdal.Open(template)
    [cols, rows] = arr.shape
    trans = data.GetGeoTransform()
    proj = data.GetProjection()
    # nodatav = 0 #data.GetNoDataValue()
    # Create the file, using the information from the template file
    outdriver = gdal.GetDriverByName("GTiff")
    # http://www.gdal.org/gdal_8h.html
    # GDT_Byte = 1, GDT_UInt16 = 2, GDT_UInt32 = 4, GDT_Int32 = 5, GDT_Float32 = 6,
    outdata   = outdriver.Create(str(out_path), rows, cols, 1, gdal.GDT_UInt32)
    # Write the array to the file, which is the original array in this example
    outdata.GetRasterBand(1).WriteArray(arr)
    # Set a no data value if required
    outdata.GetRasterBand(1).SetNoDataValue(no_data)
    # Georeference the image
    outdata.SetGeoTransform(trans)
    # Write projection information
    outdata.SetProjection(proj)
    return arr


# raster to vector 
Example 2
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 
Example 3
Project: Neural-Network-Remote-Sensing-Training-Framework   Author: Thimm   File: Predict_B.py    GNU General Public License v3.0 6 votes vote down vote up
def getHDF5(path):
    f = h5py.File(path, 'r')
    data = f['/data']
    label = np.array(f['/label'], dtype=np.int)
    newdata = np.empty((900,900, 4), dtype=np.float32)
    mean_pixel = np.array([102.93, 111.36, 116.52, 116.52], dtype=np.float32)
    newdata[...] =  data[()].transpose([2,3,1,0])[:,:,:,0] + mean_pixel
    print newdata[:,:,1:4].shape
    im = Image.fromarray(np.uint8(newdata[:,:,1:4]))
    im.save('Orig.jpg', 'JPEG')
    dst_ds = gdal.GetDriverByName('GTiff').Create('myGeoTIFF.tif', newdata.shape[0], newdata.shape[1], newdata.shape[2], gdal.GDT_Byte)
    print newdata
    dst_ds.GetRasterBand(2).WriteArray(newdata[:,:,0])   # write r-band to the raster
    dst_ds.GetRasterBand(1).WriteArray(newdata[:,:,1])   # write g-band to the raster
    dst_ds.GetRasterBand(0).WriteArray(newdata[:,:,2])
    dst_ds.FlushCache()                     # write to disk
    dst_ds = None
    return 
Example 4
Project: surfclass   Author: Kortforsyningen   File: rasterio.py    MIT License 6 votes vote down vote up
def gdaltype_to_creationoptions(gdaltype):
    """Gets GDAL creation options suitable for a GDAL datatype.

    Args:
        gdaltype (int): GDAL datatype value. For instance `osgeo.gdal.GDT_Byte`

    Raises:
        NotImplementedError: If given GDAL data type is not supported.

    Returns:
        list of str: List of GDAL creation options.

    """
    if gdaltype in [
        gdal.GDT_Byte,
        gdal.GDT_Int16,
        gdal.GDT_Int32,
        gdal.GDT_UInt16,
        gdal.GDT_UInt32,
    ]:
        return gdal_int_options
    if gdaltype in [gdal.GDT_Float32, gdal.GDT_Float64]:
        return gdal_float_options
    raise NotImplementedError() 
Example 5
Project: coast-def   Author: zdb999   File: geo.py    GNU General Public License v3.0 6 votes vote down vote up
def make_extent_layer(dem_path, water_points, height, out_path = "flood_extent"):

  # Open DEM layer

  raster, img, projection, transform = utils.import_dem(dem_path)

  # TODO: implement water points, change this

  points = water_points

  # Get xy coorinates

  h,w = np.shape(img)

  # Make output file
  driver = gdal.GetDriverByName('GTiff')
  outRaster = driver.Create(out_path, w, h, 1, gdal.GDT_Byte)
  outband = outRaster.GetRasterBand(1)
  outRaster.SetGeoTransform(transform)
  outRaster.SetProjection(projection)

  # Find extent and write it to file
  outband.WriteArray(flood_extent(img, points, 10))
  outband.FlushCache() 
Example 6
Project: fire-data-processing   Author: ANU-WALD   File: update_flammability_mosaic.py    GNU Affero General Public License v3.0 6 votes vote down vote up
def update_flammability_mosaic(date, n_band, dst, tmp, comp):
            
    flam = compose_mosaic(date, n_band, "flammability", gdal.GDT_Float32)
    anom = compose_mosaic(date, n_band, "anomaly", gdal.GDT_Float32)
    qmask = compose_mosaic(date, n_band, "quality_mask", gdal.GDT_Byte)
    
    tmp_file = os.path.join(tmp, uuid.uuid4().hex + ".nc")
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    pack_flammability_mosaic(d, flam, anom, qmask, tmp_file)

    if not os.path.isfile(dst):
        shutil.move(tmp_file, dst)
    else:
        tmp_file2 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
        os.system("cdo mergetime {} {} {}".format(tmp_file, dst, tmp_file2))
        os.remove(tmp_file)
        if comp:
            # In case the UNLIMITED time dimension is causing problems. nccopy has: 
            # [-u] convert unlimited dimensions to fixed size in output
            tmp_file3 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
            os.system("nccopy -d 4 -c 'time/4,y/240,x/240' {} {}".format(tmp_file2, tmp_file3))
            os.remove(tmp_file2)
            shutil.move(tmp_file3, dst)
        else: 
            shutil.move(tmp_file2, dst) 
Example 7
Project: fire-data-processing   Author: ANU-WALD   File: update_fmc_mosaic.py    GNU Affero General Public License v3.0 6 votes vote down vote up
def update_fmc_mosaic(date, n_band, dst, tmp, comp):
            
    fmc_mean = compose_mosaic(date, n_band, "lfmc_mean", gdal.GDT_Float32)
    fmc_stdv = compose_mosaic(date, n_band, "lfmc_stdv", gdal.GDT_Float32)
    q_mask = compose_mosaic(date, n_band, "quality_mask", gdal.GDT_Byte)
    
    tmp_file = os.path.join(tmp, uuid.uuid4().hex + ".nc")
    d = datetime.utcfromtimestamp(date.astype('O')/1e9)
    pack_fmc_mosaic(d, fmc_mean, fmc_stdv, q_mask, tmp_file)

    if not os.path.isfile(dst):
        shutil.move(tmp_file, dst)
    else:
        tmp_file2 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
        os.system("cdo mergetime {} {} {}".format(tmp_file, dst, tmp_file2))
        os.remove(tmp_file)
        if comp:
            # In case the UNLIMITED time dimension is causing problems. nccopy has: 
            # [-u] convert unlimited dimensions to fixed size in output
            tmp_file3 = os.path.join(tmp, uuid.uuid4().hex + ".nc")
            os.system("nccopy -d 4 -c 'time/4,y/240,x/240' {} {}".format(tmp_file2, tmp_file3))
            os.remove(tmp_file2)
            shutil.move(tmp_file3, dst)
        else: 
            shutil.move(tmp_file2, dst) 
Example 8
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def rasterize(self, griddef=None, layer=0, invert=False, asVar=False, ldebug=False):
    ''' "burn" shapefile on a 2D raster; returns a 2D boolean array '''
    if griddef.__class__.__name__ != GridDefinition.__name__: raise TypeError 
    #if not isinstance(griddef,GridDefinition): raise TypeError # this is always False. probably due to pickling
    if not isinstance(invert,(bool,np.bool)): raise TypeError
    # fill values
    if invert: inside, outside = 0,1
    else: inside, outside = 1,0
    shp_lyr = self.getLayer(layer) # get shape layer
    # create raster to burn shape onto
    if ldebug: print(' - creating raster')
    msk_ds = ramdrv.Create(self.name, griddef.size[0], griddef.size[1], 1, gdal.GDT_Byte)
    # N.B.: this is a special case: only one band (1) and always boolean (gdal.GDT_Byte)
    # set projection parameters
    msk_ds.SetGeoTransform(griddef.geotransform)  # does the order matter?
    msk_ds.SetProjection(griddef.projection.ExportToWkt())  # is .ExportToWkt() necessary?
    # initialize raster band        
    msk_rst = msk_ds.GetRasterBand(1) # only one anyway...
    msk_rst.Fill(outside); msk_rst.SetNoDataValue(outside) # fill with zeros
    # burn shape layer onto raster band
    if ldebug: print(' - burning layer to raster')
    err = gdal.RasterizeLayer(msk_ds, [1], shp_lyr, burn_values = [inside]) # None, None, [1] # burn_value = 1
    # use argument ['ALL_TOUCHED=TRUE'] like so: None, None, [1], ['ALL_TOUCHED=TRUE']
    if err != 0: raise GDALError('ERROR CODE %i'%err)
    #msk_ds.FlushCash()  
    # retrieve mask array from raster band
    if ldebug: print(' - retrieving mask')
    mask = msk_ds.GetRasterBand(1).ReadAsArray()
    # convert to Variable object, is desired
    if asVar: 
      mask = Variable(name=self.name, units='mask', axes=(griddef.ylat,griddef.xlon), data=mask, 
                      dtype=np.bool, mask=None, fillValue=outside, atts=None, plot=None)
      mask = addGDALtoVar(mask, griddef=griddef,) # add GDAL to mask       
    # return mask array
    return mask  

# a container class that also acts as a shape for the outline 
Example 9
Project: PostgreSQL-GeoPackage   Author: EOX-A   File: gpkg-pg_dump.py    MIT License 5 votes vote down vote up
def create_gpkg(
    gpkg_name, proj_string, size=(1, 1), geotransform=[0, 1, 0, 0, 0, -1],
    creation_options=None
):
    if os.path.exists("%s.gpkg" % gpkg_name):
        sys.stderr.write(
            "ERROR: SQLite GeoPackage '%s.gpkg' already exists.\n" % gpkg_name
        )
        sys.exit(1)

    gdal.AllRegister()
    drv = gdal.GetDriverByName("GPKG")
    try:
        gpkg = drv.Create(
            "%s.gpkg" % gpkg_name, size[0], size[1], 1, gdal.GDT_Byte,
            creation_options
        )

        proj = osr.SpatialReference()
        res = proj.SetWellKnownGeogCS(proj_string)
        if res != 0:
            if proj_string[0:4] == 'EPSG':
                proj.ImportFromEPSG(int(proj_string[5:]))
        gpkg.SetProjection(proj.ExportToWkt())
        gpkg.SetGeoTransform(geotransform)
        gpkg = None
    except Exception as e:
        sys.stderr.write(
            "ERROR: Cannot create SQLite GeoPackage '%s.gpkg'. "
            "Error message was: '%s'.\n" % (gpkg_name, e.message)
        )
        sys.exit(1) 
Example 10
Project: buzzard   Author: airware   File: _gdal_gdt_conv.py    Apache License 2.0 5 votes vote down vote up
def _str_of_gdt(gdt):
    return {
        gdal.GDT_Byte: 'GDT_Byte',
        gdal.GDT_Int16: 'GDT_Int16',
        gdal.GDT_Int32: 'GDT_Int32',
        gdal.GDT_UInt16: 'GDT_UInt16',
        gdal.GDT_UInt32: 'GDT_UInt32',
        gdal.GDT_Float32: 'GDT_Float32',
        gdal.GDT_Float64: 'GDT_Float64',
        gdal.GDT_CFloat32: 'GDT_CFloat32',
        gdal.GDT_CFloat64: 'GDT_CFloat64',
        gdal.GDT_CInt16: 'GDT_CInt16',
        gdal.GDT_CInt32: 'GDT_CInt32',
    }[gdt] 
Example 11
Project: lidar   Author: giswqs   File: filling.py    MIT License 5 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 
Example 12
Project: coded   Author: bullocke   File: classify.py    MIT License 5 votes vote down vote up
def write_geotiff(fname, data, geo_transform, projection, data_type=gdal.GDT_Byte):

    """
    Create a GeoTIFF file with the given data.
    :param fname: Path to a directory with shapefiles
    :param data: Number of rows of the result
    :param geo_transform: Returned value of 
	gdal.Dataset.GetGeoTransform (coefficients for transforming between 
	pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.
    :param projection: Projection definition string (Returned by 
	gdal.Dataset.GetProjectionRef)
    """

    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, data_type)
    dataset.SetGeoTransform(geo_transform)
    dataset.SetProjection(projection)
    band = dataset.GetRasterBand(1)
    band.WriteArray(data)

    metadata = {
        'TIFFTAG_COPYRIGHT': 'CC BY 4.0',
        'TIFFTAG_DOCUMENTNAME': 'classification',
        'TIFFTAG_IMAGEDESCRIPTION': 'Supervised classification.',
        'TIFFTAG_SOFTWARE': 'Python, GDAL, scikit-learn'
    }
    dataset.SetMetadata(metadata)

    dataset = None  # Close the file
    return 
Example 13
Project: ALCD   Author: CNES   File: masks_preprocessing.py    GNU General Public License v3.0 5 votes vote down vote up
def rasterize_shp(input_shp, out_tif, reference_tif):
    '''
    From a shapefile, rasterize it
    '''
    gdalformat = 'GTiff'
    datatype = gdal.GDT_Byte
    burnVal = 0  # value for the output image pixels

    # Get projection info from reference image
    image = gdal.Open(reference_tif, gdal.GA_ReadOnly)
    print(image)

    # Open Shapefile
    shapefile = ogr.Open(input_shp)
    shapefile_layer = shapefile.GetLayer()

    # Rasterise
    output = gdal.GetDriverByName(gdalformat).Create(
        out_tif, image.RasterXSize, image.RasterYSize, 1, datatype, options=['COMPRESS=DEFLATE'])
    output.SetProjection(image.GetProjectionRef())
    output.SetGeoTransform(image.GetGeoTransform())

    # Write data to band 1
    band = output.GetRasterBand(1)
    band.SetNoDataValue(1)
    gdal.RasterizeLayer(output, [1], shapefile_layer, burn_values=[burnVal])

    # Close datasets
    band = None
    output = None
    image = None
    shapefile = None 
Example 14
Project: qgis-rasterarray   Author: rrwen   File: RasterArray.py    MIT License 5 votes vote down vote up
def Array2Raster(inArray,
                 outRaster,
                 rasterOrigin,
                 pixelWidth,
                 pixelHeight,
                 EPSG):
    
    # (H1.1) Obtain Array Information
    cols = inArray.shape[1]
    rows = inArray.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]
    
    # (H1.2) Write Array to Raster       
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(outRaster,
                              cols,
                              rows,
                              1,
                              gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX,
                               pixelWidth,
                               0,
                               originY,
                               0,
                               pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(inArray,0,0)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(EPSG)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
    outband.SetNoDataValue(-99)
    
    # (H1.3) Reset
    driver = None
    outRaster = None
    outBand = None
    inArray = None 
Example 15
Project: Semantic_Segmentation_Keras   Author: liuph0119   File: predicting.py    Apache License 2.0 5 votes vote down vote up
def save_prediction(arr, dst_path, to_prob=False, with_geo=False, spatial_reference=None):
    """ save prediction.
    :param arr: 2D/3D array.
    :param dst_path: string.
    :param to_prob: bool, whether to save probability. If True, a *.npy or a float-formatted tiff file will be created.
    :param with_geo: bool, whether to save spatial reference.
    :param spatial_reference: string, or dict. If the data type is string, this must be a path of a existing tiff file,
        else, must be a dict contains "geotransform" and "projection"
    :return:
    """
    if with_geo:
        if type(spatial_reference) is str and os.path.exists(spatial_reference):
            img_info = get_image_info(spatial_reference, get_geotransform=True, get_projection=True)
        elif type(spatial_reference) is dict \
                and "geotransform" in spatial_reference and "projection" in spatial_reference:
            img_info = spatial_reference
        else:
            raise ValueError("Invalid 'spatial_reference': None. Expected to be a dict or a existing tiff image path.")
        if to_prob:
            save_to_image_gdal(arr, dst_path, datatype=gdal.GDT_Float32,
                               geoTransform=img_info["geotransform"], proj=img_info["projection"])
        else:
            save_to_image_gdal(arr, dst_path, datatype=gdal.GDT_Byte,
                               geoTransform=img_info["geotransform"], proj=img_info["projection"])
    else:
        if to_prob:
            if ".npy" not in dst_path:
                dst_path = dst_path + ".npy"
            np.save(dst_path, arr)
        else:
            save_to_image(arr, dst_path) 
Example 16
Project: Semantic_Segmentation_Keras   Author: liuph0119   File: image_io_utils.py    Apache License 2.0 5 votes vote down vote up
def save_to_image_gdal(arr, image_path, datatype=gdal.GDT_Byte, geoTransform = (0,1,0,0,0,-1), proj=None, nodata=None):
    """ save to geo-tiff image
    :param arr: array of shape (height, width, 3) or (height, width)
    :param image_path: string
    :param datatype: data type of the geo-tiff image, default gdal.GDT_Byte.
    :param geoTransform: tuple, default (0,1,0,0,0,-1)
        the geo-transform parameters.
    :param proj: string, default None
        the spatial projection wkt string. If None, no spatial projection will added to the output image
    :param nodata: float, default None
        no-data value. If None, no no-data value will be set

    :return: None
    """
    datatype = gdal.GDT_Byte if datatype is None else datatype

    if arr.ndim == 2:
        arr = np.expand_dims(arr, -1)
    if arr.ndim != 3:
        raise ValueError("[save_to_image_gdal: the input array must have 2 or 3 dimensions!!!]")

    nBands = arr.shape[-1]
    nRows, nCols = arr.shape[:2]

    driver = gdal.GetDriverByName("GTiff")
    if os.path.exists(image_path):
        os.remove(image_path)

    outRaster = driver.Create(image_path, nCols, nRows, nBands, datatype)
    outRaster.SetGeoTransform(geoTransform)
    if proj is not None:
        outRasterSRS = osr.SpatialReference()
        outRasterSRS.ImportFromWkt(proj)
        outRaster.SetProjection(outRasterSRS.ExportToWkt())

    for _bandNum in range(nBands):
        outBand = outRaster.GetRasterBand(_bandNum+1)
        if nodata is not None:
            outBand.SetNoDataValue(nodata)
        outBand.WriteArray(arr[:,:, _bandNum])
        outBand.FlushCache() 
Example 17
Project: cv4ag   Author: worldbank   File: gdal_rasterize.py    MIT License 5 votes vote down vote up
def create_raster(xsize, ysize, driver='MEM',tmpfile='', gt=None, srs_wkt=None, nodata=None,  init=None, datatype=gdal.GDT_Byte):
    # Create a memory raster to rasterize into.
    out_ds = gdal.GetDriverByName(driver).Create(tmpfile, xsize, ysize, 1 ,datatype)
    if init is not None:out_ds.GetRasterBand(1).Fill(init)
    if nodata is not None:out_ds.GetRasterBand(1).SetNoDataValue(nodata)
    if gt:out_ds.SetGeoTransform(gt)
    if srs_wkt:out_ds.SetProjection(srs_wkt)
    return out_ds 
Example 18
Project: PyGeoC   Author: lreis2415   File: raster.py    MIT License 5 votes vote down vote up
def raster_reclassify(srcfile, v_dict, dstfile, gdaltype=GDT_Float32):
        """Reclassify raster by given classifier dict.

        Args:
            srcfile: source raster file.
            v_dict: classifier dict.
            dstfile: destination file path.
            gdaltype (:obj:`pygeoc.raster.GDALDataType`): GDT_Float32 as default.
        """
        src_r = RasterUtilClass.read_raster(srcfile)
        src_data = src_r.data
        dst_data = numpy.copy(src_data)
        if gdaltype == GDT_Float32 and src_r.dataType != GDT_Float32:
            gdaltype = src_r.dataType
        no_data = src_r.noDataValue
        new_no_data = DEFAULT_NODATA
        if gdaltype in [GDT_Unknown, GDT_Byte, GDT_UInt16, GDT_UInt32]:
            new_no_data = 0
        if not MathClass.floatequal(new_no_data, src_r.noDataValue):
            if src_r.noDataValue not in v_dict:
                v_dict[src_r.noDataValue] = new_no_data
                no_data = new_no_data

        for (k, v) in iteritems(v_dict):
            dst_data[src_data == k] = v
        RasterUtilClass.write_gtiff_file(dstfile, src_r.nRows, src_r.nCols, dst_data,
                                         src_r.geotrans, src_r.srs, no_data, gdaltype) 
Example 19
Project: computing-pipeline   Author: terraref   File: Get_FLIR.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def create_geotiff(np_arr, gps_bounds, out_file_path, base_dir):
    try:
        nrows,ncols = np.shape(np_arr)
        # gps_bounds: (lat_min, lat_max, lng_min, lng_max)
        xres = (gps_bounds[3] - gps_bounds[2])/float(ncols)
        yres = (gps_bounds[1] - gps_bounds[0])/float(nrows)
        geotransform = (gps_bounds[2],xres,0,gps_bounds[1],0,-yres) #(top left x, w-e pixel resolution, rotation (0 if North is up), top left y, rotation (0 if North is up), n-s pixel resolution)

        output_path = out_file_path

        output_raster = gdal.GetDriverByName('GTiff').Create(output_path, ncols, nrows, 3, gdal.GDT_Byte)

        output_raster.SetGeoTransform(geotransform) # specify coordinates
        srs = osr.SpatialReference() # establish coordinate encoding
        srs.ImportFromEPSG(4326) # specifically, google mercator
        output_raster.SetProjection( srs.ExportToWkt() ) # export coordinate system to file

        # TODO: Something wonky w/ uint8s --> ending up w/ lots of gaps in data (white pixels)
        output_raster.GetRasterBand(1).WriteArray(np_arr.astype('uint8')) # write red channel to raster file
        output_raster.GetRasterBand(1).FlushCache()
        output_raster.GetRasterBand(1).SetNoDataValue(-99)
        
        output_raster.GetRasterBand(2).WriteArray(np_arr.astype('uint8')) # write green channel to raster file
        output_raster.GetRasterBand(2).FlushCache()
        output_raster.GetRasterBand(2).SetNoDataValue(-99)

        output_raster.GetRasterBand(3).WriteArray(np_arr.astype('uint8')) # write blue channel to raster file
        output_raster.GetRasterBand(3).FlushCache()
        output_raster.GetRasterBand(3).SetNoDataValue(-99)


        # for test: once we've saved the image, make sure to append this path to our list of TIFs
        tif_file = os.path.join(base_dir, tif_list_file)
        f = open(tif_file,'a+')
        f.write(output_path + '\n')
    except Exception as ex:
        fail('Error creating GeoTIFF: ' + str(ex)) 
Example 20
Project: geobricks_qgis_plugin_trmm   Author: geobricks   File: conversions.py    GNU General Public License v2.0 5 votes vote down vote up
def Byte(dataset_or_band):
    return ConvertedDataset(dataset_or_band, gdal.GDT_Byte) 
Example 21
Project: surfclass   Author: Kortforsyningen   File: rasterio.py    MIT License 4 votes vote down vote up
def read_2d(self, geom):
        """Reads part of the raster into a 2D MaskedArray with a mask marking cells outside the polygon.

        This is suitable for doing analysis on cells inside a given geometry object.

        Note: The mask marks cells outside the geometry. This means that cells inside the geometry are NOT masked
            even if they are equal to the raster nodata value.

        Args:
            geom (osgeo.ogr.Geometry): OGR Geometry object

        Raises:
            TypeError: If geometry is not an `osgeo.ogr.Geometry`
            ValueError: If bbox of the geometry is entirely or partly outside raster coverage.

        Returns:
            [numpy.ma.maskedArray]: A masked array where cells outside the geometry are masked.

        """
        if not isinstance(geom, ogr.Geometry):
            raise TypeError("Must be OGR geometry")
        mem_type = geom.GetGeometryType()
        ogr_env = geom.GetEnvelope()
        geom_bbox = Bbox(ogr_env[0], ogr_env[2], ogr_env[1], ogr_env[3])
        window = self.bbox_to_pixel_window(geom_bbox)
        if window[2] <= 0 or window[3] <= 0:
            return np.ma.empty(shape=(0, 0))
        src_array = self.read_raster(window=window, masked=False)
        # calculate new geotransform of the feature subset
        new_gt = self.window_geotransform(window)

        # Create a temporary vector layer in memory
        mem_ds = self._ogr_mem_drv.CreateDataSource("out")
        mem_layer = mem_ds.CreateLayer("mem_lyr", geom.GetSpatialReference(), mem_type)
        mem_feature = ogr.Feature(mem_layer.GetLayerDefn())
        mem_feature.SetGeometry(geom)
        mem_layer.CreateFeature(mem_feature)

        # Rasterize the feature
        mem_raster_ds = self._gdal_mem_drv.Create(
            "", window[2], window[3], 1, gdal.GDT_Byte
        )
        mem_raster_ds.SetProjection(self._ds.GetProjection())
        mem_raster_ds.SetGeoTransform(new_gt)
        # Burn 1 inside our feature
        gdal.RasterizeLayer(mem_raster_ds, [1], mem_layer, burn_values=[1])
        rasterized_array = mem_raster_ds.ReadAsArray()

        # Mask the source data array with our current feature mask
        masked = np.ma.MaskedArray(src_array, mask=np.logical_not(rasterized_array))
        return masked 
Example 22
Project: Land_surface_temperature_Landsat   Author: dimejimudele   File: split_window_algorithm_Land_Surface_Temperture_Landsat8.py    MIT License 4 votes vote down vote up
def array2raster(newRasterfn, dataset, array, dtype):
    """
    save GTiff file from numpy.array
    input:
        newRasterfn: name to save file with
        dataset : original tif file to obtain geo information. You should use the Level 1 quantized and calibrated scaled Digital Numbers (DN) TIR band data (e.g Band 10 landsat 8 data)
        array : The Land surface temperature array
        dtype: Byte or Float32.
    """
    cols = array.shape[1]
    rows = array.shape[0]
    originX, pixelWidth, b, originY, d, pixelHeight = dataset.GetGeoTransform()

    driver = gdal.GetDriverByName('GTiff')

    # set data type to save.
    GDT_dtype = gdal.GDT_Unknown
    if dtype == "Byte":
        GDT_dtype = gdal.GDT_Byte
    elif dtype == "Float32":
        GDT_dtype = gdal.GDT_Float32
    else:
        print("Not supported data type.")

    # set number of band.
    if array.ndim == 2:
        band_num = 1
    else:
        band_num = array.shape[2]

    outRaster = driver.Create(newRasterfn, cols, rows, band_num, GDT_dtype)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

    # Loop over all bands.
    for b in range(band_num):
        outband = outRaster.GetRasterBand(b + 1)
        # Read in the band's data into the third dimension of our array
        if band_num == 1:
            outband.WriteArray(array)
        else:
            outband.WriteArray(array[:,:,b])

    # setteing srs from input tif file.
    prj=dataset.GetProjection()
    outRasterSRS = osr.SpatialReference(wkt=prj)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache() 
Example 23
Project: Land_surface_temperature_Landsat   Author: dimejimudele   File: mono_window_LST_computation.py    MIT License 4 votes vote down vote up
def array2raster(newRasterfn, dataset, array, dtype):
    """
    save GTiff file from numpy.array
    input:
        newRasterfn: name to save file with
        dataset : original tif file to obtain geo information. You should use the Level 1 quantized and calibrated scaled Digital Numbers (DN) TIR band data (e.g Band 10 landsat 8 data)
        array : The Land surface temperature array
        dtype: Byte or Float32.
    """
    cols = array.shape[1]
    rows = array.shape[0]
    originX, pixelWidth, b, originY, d, pixelHeight = dataset.GetGeoTransform()

    driver = gdal.GetDriverByName('GTiff')

    # set data type to save.
    GDT_dtype = gdal.GDT_Unknown
    if dtype == "Byte":
        GDT_dtype = gdal.GDT_Byte
    elif dtype == "Float32":
        GDT_dtype = gdal.GDT_Float32
    else:
        print("Not supported data type.")

    # set number of band.
    if array.ndim == 2:
        band_num = 1
    else:
        band_num = array.shape[2]

    outRaster = driver.Create(newRasterfn, cols, rows, band_num, GDT_dtype)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

    # Loop over all bands.
    for b in range(band_num):
        outband = outRaster.GetRasterBand(b + 1)
        # Read in the band's data into the third dimension of our array
        if band_num == 1:
            outband.WriteArray(array)
        else:
            outband.WriteArray(array[:,:,b])

    # setteing srs from input tif file.
    prj=dataset.GetProjection()
    outRasterSRS = osr.SpatialReference(wkt=prj)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache() 
Example 24
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 4 votes vote down vote up
def shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None):
    """Rasterize input shapefile to match existing raster Dataset (or specified res/extent/t_srs)
    """
    shp_ds = ogr.Open(shp_fn)
    lyr = shp_ds.GetLayer()
    #This returns xmin, ymin, xmax, ymax
    shp_extent = lyr_extent(lyr)
    shp_srs = lyr.GetSpatialRef()
    # dst_dt = gdal.GDT_Byte
    ndv = 0
    if r_ds is not None:
        r_extent = ds_extent(r_ds)
        res = get_res(r_ds, square=True)[0] 
        if extent is None:
            extent = r_extent
        r_srs = get_ds_srs(r_ds)
        r_geom = ds_geom(r_ds)
        # dst_ns = r_ds.RasterXSize
        # dst_nl = r_ds.RasterYSize
        #Convert raster extent to shp_srs
        cT = osr.CoordinateTransformation(r_srs, shp_srs)
        r_geom_reproj = geom_dup(r_geom)
        r_geom_reproj.Transform(cT)
        r_geom_reproj.AssignSpatialReference(t_srs)
        lyr.SetSpatialFilter(r_geom_reproj)
        #lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
    else:
        #TODO: clean this up
        if res is None:
            sys.exit("Must specify input res")
        if extent is None:
            print("Using input shp extent")
            extent = shp_extent
    if t_srs is None:
        t_srs = r_srs
    if not shp_srs.IsSame(t_srs):
        print("Input shp srs: %s" % shp_srs.ExportToProj4())
        print("Specified output srs: %s" % t_srs.ExportToProj4())
        out_ds = lyr_proj(lyr, t_srs)
        outlyr = out_ds.GetLayer()
    else:
        outlyr = lyr
    #outlyr.SetSpatialFilter(r_geom)
    m_ds = mem_ds(res, extent, srs=t_srs, dtype=gdal.GDT_Byte)
    b = m_ds.GetRasterBand(1)
    b.SetNoDataValue(ndv)
    gdal.RasterizeLayer(m_ds, [1], outlyr, burn_values=[1])
    a = b.ReadAsArray()
    a = ~(a.astype('Bool'))
    return a 
Example 25
Project: HistoricalMap   Author: nkarasiak   File: function_dataraster.py    GNU General Public License v2.0 4 votes vote down vote up
def open_data_band(filename):
    """[email protected] The function open and load the image given its name. 
    The function open and load the image given its name. 
    The type of the data is checked from the file and the scipy array is initialized accordingly.
        Input:
            filename: the name of the file
        Output:
            data : the opened data with gdal.Open() method
            im : empty table with right dimension (array)
    
    """
    data = gdal.Open(filename,gdal.GA_ReadOnly)
    if data is None:
        print 'Impossible to open '+filename
        exit()
    nc = data.RasterXSize
    nl = data.RasterYSize
#    d  = data.RasterCount
    
    # Get the type of the data
    gdal_dt = data.GetRasterBand(1).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()
    
    # Initialize the array
    im = sp.empty((nl,nc),dtype=dt) 
    return data,im 
Example 26
Project: HistoricalMap   Author: nkarasiak   File: function_dataraster.py    GNU General Public License v2.0 4 votes vote down vote up
def create_empty_tiff(outname,im,d,GeoTransform,Projection):
    '''[email protected] Write an empty image on the hard drive.
    
    Input: 
        outname: the name of the file to be written
        im: the image cube
        GeoTransform: the geotransform information 
        Projection: the projection information
    Output:
        Nothing --
    '''
    nl = im.shape[0]
    nc = im.shape[1]

    driver = gdal.GetDriverByName('GTiff')
    dt = im.dtype.name
    # Get the data type
    if dt == 'bool' or dt == 'uint8':
        gdal_dt=gdal.GDT_Byte
    elif dt == 'int8' or dt == 'int16':
        gdal_dt=gdal.GDT_Int16
    elif dt == 'uint16':
        gdal_dt=gdal.GDT_UInt16
    elif dt == 'int32':
        gdal_dt=gdal.GDT_Int32
    elif dt == 'uint32':
        gdal_dt=gdal.GDT_UInt32
    elif dt == 'int64' or dt == 'uint64' or dt == 'float16' or dt == 'float32':
        gdal_dt=gdal.GDT_Float32
    elif dt == 'float64':
        gdal_dt=gdal.GDT_Float64
    elif dt == 'complex64':
        gdal_dt=gdal.GDT_CFloat64
    else:
        print 'Data type non-suported'
        exit()
    
    dst_ds = driver.Create(outname,nc,nl, d, gdal_dt)
    dst_ds.SetGeoTransform(GeoTransform)
    dst_ds.SetProjection(Projection)
    
    return dst_ds
    
    '''
    Old function that cannot manage to write on each band outside the script
    '''
#    if d==1:
#        out = dst_ds.GetRasterBand(1)
#        out.WriteArray(im)
#        out.FlushCache()
#    else:
#        for i in range(d):
#            out = dst_ds.GetRasterBand(i+1)
#            out.WriteArray(im[:,:,i])
#            out.FlushCache()
#    dst_ds = None 
Example 27
Project: HistoricalMap   Author: nkarasiak   File: function_historical_map.py    GNU General Public License v2.0 4 votes vote down vote up
def postClassRaster(self,inRaster,sieveSize,inClassNumber,outShp):        
        """ [email protected] Sieve size with gdal.Sieve() fiunction, them reclass to delete unwanted labels """
        
        try:
            rasterTemp = tempfile.mktemp('.tif')
                    
            datasrc = gdal.Open(inRaster)
            srcband = datasrc.GetRasterBand(1)
            data,im=dataraster.open_data_band(inRaster)        
            
            drv = gdal.GetDriverByName('GTiff')
            dst_ds = drv.Create(rasterTemp,datasrc.RasterXSize,datasrc.RasterXSize,1,gdal.GDT_Byte)
            
            dst_ds.SetGeoTransform(datasrc.GetGeoTransform())
            dst_ds.SetProjection(datasrc.GetProjection())
        
            dstband=dst_ds.GetRasterBand(1)
            
            
            
            def sieve(srcband,dstband,sieveSize):
                gdal.SieveFilter(srcband,None,dstband,sieveSize,8)
            
            pixelSize = datasrc.GetGeoTransform()[1] #get pixel size
            pixelSieve = int(sieveSize/(pixelSize*pixelSize)) #get number of pixel to sieve
            
            sieve(srcband,dstband,pixelSieve)
            
            dst_ds = None # close destination band
            
            self.Progress.addStep()
            

            rasterTemp = self.reclassAndFillHole(rasterTemp,inClassNumber)
            
            
            self.Progress.addStep()
            
            
            
        except:
            QgsMessageLog.logMessage("Cannot sieve with raster function")

        
        try :
            outShp = self.polygonize(rasterTemp,outShp) # vectorize raster
        except :
            self.Progress.reset()    
        
        self.Progress.addStep()        
        
        self.Progress.reset()
        return outShp