Python osgeo.gdal.GetDriverByName() Examples

The following are code examples for showing how to use osgeo.gdal.GetDriverByName(). 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: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def getShapefileRowHeaders(cls, shapefileFilename):
        """
        @summary: Get a (sorted by feature id) list of row headers for a shapefile
        @todo: Move this to a common Vector class in LmCommon or LmBackend and 
                 use with LmCompute.  This is a rough copy of what is now used for 
                 rasterIntersect.
        """
        ogr.RegisterAll()
        drv = ogr.GetDriverByName(LMFormat.getDefaultOGR().driver)
        ds = drv.Open(shapefileFilename)
        lyr = ds.GetLayer(0)
        
        rowHeaders = []
        
        for j in range(lyr.GetFeatureCount()):
            curFeat = lyr.GetFeature(j)
            siteIdx = curFeat.GetFID()
            x, y = curFeat.geometry().Centroid().GetPoint_2D()
            rowHeaders.append((siteIdx, x, y))
            
        return sorted(rowHeaders)

# ............................................... 
Example 2
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def testVector(dlocation, driver=LMFormat.getDefaultOGR().driver):
        goodData = False
        featCount = 0
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(driver)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                goodData = False
            else:
                try:
                    slyr = ds.GetLayer(0)
                except Exception, e:
                    goodData = False
                else: 
Example 3
Project: buzzard   Author: airware   File: tools.py    Apache License 2.0 6 votes vote down vote up
def make_tif(path, tloffset=(0, 0), reso=(0.25, -0.25), rsize=(20, 10),
             proj=SRS[0]['wkt'], channel_count=1, dtype=gdal.GDT_Float32):
    """Create a tiff files and return info about it"""
    tl = ROOT_TL + tloffset
    reso = np.asarray(reso)
    fp = buzz.Footprint(tl=tl, rsize=rsize, size=np.abs(reso * rsize))
    x, y = fp.meshgrid_spatial
    x = np.abs(x) - abs(ROOT_TL[0])
    y = abs(ROOT_TL[1]) - np.abs(y)
    x *= 15
    y *= 15
    a = x / 2 + y / 2
    a = np.around(a).astype('float32')
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(path, rsize[0], rsize[1], channel_count, dtype)
    dataset.SetGeoTransform(fp.gt)
    dataset.SetProjection(proj)
    for i in range(channel_count):
        dataset.GetRasterBand(i + 1).WriteArray(a)
        dataset.GetRasterBand(i + 1).SetNoDataValue(-32000.)
    dataset.FlushCache()
    return path, fp, a 
Example 4
Project: Landsat8LST_SWA   Author: eduard-kazakov   File: l8_lst_swa_common_lib.py    GNU General Public License v2.0 6 votes vote down vote up
def save_nparray_as_raster (nparray, driver_name, cell_type, base_raster_XSize, base_raster_YSize, base_raster_projection, base_raster_transform, out_path):
    """
    Save NumPy nparray as raster via GDAL library methods
    :param nparray: input nparray
    :param driver_name: driver name, e.g. GeoTIff
    :param cell_type: type of cell, e.g. Float32
    :param base_raster_XSize: X-resolution of output raster
    :param base_raster_YSize: Y-resolution of output raster
    :param base_raster_projection: Geographic projection for output raster
    :param base_raster_transform: Geographic transformation type for output raster
    :param out_path: where output raster will be saved
    """
    cols = base_raster_XSize
    rows = base_raster_YSize
    bands = 1
    dt = cell_type
    driver = gdal.GetDriverByName(driver_name)
    out_data = driver.Create(out_path,cols,rows,bands,dt)
    out_data.SetProjection (base_raster_projection)
    out_data.SetGeoTransform (base_raster_transform)

    out_data.GetRasterBand(1).WriteArray (nparray) 
Example 5
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 6
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 7
Project: watools   Author: wateraccounting   File: data_conversions.py    Apache License 2.0 6 votes vote down vote up
def Convert_bil_to_tiff(input_bil, output_tiff):
    """
    This function converts the bil files into tiff files

    Keyword Arguments:
    input_bil -- name, name of the bil file
    output_tiff -- Name of the output tiff file
    """
    import gdalconst
    
    gdal.GetDriverByName('EHdr').Register()
    dest = gdal.Open(input_bil, gdalconst.GA_ReadOnly)
    Array = dest.GetRasterBand(1).ReadAsArray()
    geo_out = dest.GetGeoTransform()
    Save_as_tiff(output_tiff, Array, geo_out, "WGS84")  
    
    return(output_tiff) 
Example 8
Project: watools   Author: wateraccounting   File: data_conversions.py    Apache License 2.0 6 votes vote down vote up
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)
    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) 
Example 9
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 10
Project: gnm_qgis   Author: nextgis   File: gnm_manager.py    GNU General Public License v2.0 6 votes vote down vote up
def onRemoveNetworkClicked(self):
        #gdal.SetConfigOption('CPL_DEBUG','ON')
        #gdal.SetConfigOption('CPL_LOG','D:/nextgis/gnmmanager/gdal_log.txt')
        #gdal.SetConfigOption('CPL_LOG_ERRORS','ON')
        self.dlg_remove.show() 
        result = self.dlg_remove.my_exec_()
        if result == 1:
            self.NETWORK_DS = None
            ok_to_del = self.dlg_remove.FULLY_DELETE
            if ok_to_del:
# TODO: really pass here the selected GNM format:           
                gdal.GetDriverByName('GNMFile').Delete(str(self.NETWORK_FULLPATH))
            self.clickFlagButton(self.PRESSED_TOOLB,None) # unpress any pressed flag button and unset map tool
            self.removeGnmLayersGroup()
            self.enableMenusForNetwork(False)
            self.updateLayersToSearchForFlags()                
            self.NETWORK_FULLPATH = ''
            self.NETWORK_NAME = ''
            self.GFID_STARTFLAG = -1
            self.GFID_ENDFLAG = -1
            self.GFIDS_BLOCKFLAGS = [] 
Example 11
Project: eetc   Author: george-azzari   File: gdaltools.py    MIT License 6 votes vote down vote up
def get_polygon_corners(self, affine=True, **polyargs):
        """ Returns the extent of the polygon in georeferenced or affine coordinates.
            Georeferenced coordinates are directly read from the polygon file.
            Affine coordinates are calculated using the transformation set with the current instance of the class."""
        d = ogr.GetDriverByName(polyargs['driver'])
        poly = d.Open(polyargs['fpath'])
        layer = poly.GetLayerByName(polyargs['layer_name'])
        ext = layer.GetExtent()
        ulx, lrx, lry, uly = ext[0], ext[1], ext[2], ext[3]
        if not affine:
            return ulx, lrx, lry, uly
        if affine:
            lin_ul = self.get_affinecoord(ulx, uly)
            lin_lr = self.get_affinecoord(lrx, lry)
            top, bottom, left, right = lin_ul[1], lin_lr[1], lin_ul[0], lin_lr[0]
            return left, bottom, right, top  # this is the order for warping 
Example 12
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 13
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 14
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRaster(array, path, geometry, srs, navalue=-9999):
    """
    Writes a single array to a raster with coordinate system and geometric
    information.

    path = target path
    srs = spatial reference system
    """
    xpixels = array.shape[1]    
    ypixels = array.shape[0]
    path = path.encode('utf-8')
    image = gdal.GetDriverByName("GTiff").Create(path, xpixels, ypixels,
                                1, gdal.GDT_Float32)
    image.SetGeoTransform(geometry)
    image.SetProjection(srs)
    image.GetRasterBand(1).WriteArray(array)
    image.GetRasterBand(1).SetNoDataValue(navalue) 
Example 15
Project: Drought-Index-Portal   Author: WilliamsTravis   File: functions_v3.py    MIT License 6 votes vote down vote up
def toRasters(arraylist, path, geometry, srs):
    """
    Writes a list of 2d arrays, or a 3d array, to a series of tif files.

    Arraylist format = [[name,array],[name,array],....]
    path = target path
    geometry = gdal geometry object
    srs = spatial reference system object
    """
    if path[-2:] == "\\":
        path = path
    else:
        path = path + "\\"
    sample = arraylist[0][1]
    ypixels = sample.shape[0]
    xpixels = sample.shape[1]
    for ray in  tqdm(arraylist):
        image = gdal.GetDriverByName("GTiff").Create(os.path.join(path,
                                                              ray[0] + ".tif"),
                                    xpixels, ypixels, 1, gdal.GDT_Float32)
        image.SetGeoTransform(geometry)
        image.SetProjection(srs)
        image.GetRasterBand(1).WriteArray(ray[1]) 
Example 16
Project: scarplet   Author: stgl   File: dem.py    MIT License 6 votes vote down vote up
def save(self, filename):
        """Save grid as georeferenced TIFF
        """

        ncols = self._georef_info.nx
        nrows = self._georef_info.ny

        driver = gdal.GetDriverByName(GDAL_DRIVER_NAME)
        out_raster = driver.Create(filename, ncols, nrows, 1, self.dtype)
        out_raster.SetGeoTransform(self._georef_info.geo_transform)
        out_band = out_raster.GetRasterBand(1)
        out_band.WriteArray(self._griddata)

        proj = self._georef_info.projection.ExportToWkt()
        out_raster.SetProjection(proj)
        out_band.FlushCache() 
Example 17
Project: core   Author: lifemapper   File: create_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def write_tiff(out_filename, bbox, cell_size, data, epsg,
               nodata=DEFAULT_NODATA):
    """
    @summary: Write a GeoTiff raster layer from a numpy array
    @param out_filename: The file location to write the raster to
    @param bbox: (minx, miny, maxx, maxy) tuple of raster coordinates
    @param cell_size: The cell size, in map units, of each cell in the raster
    @param data: The data array to use for raster content
    @param epsg: The epsg code of the map projection to use for this raster
    @param nodata: A value to use for NODATA
    @todo: Probably should establish a common layer package in LmCommon
    """
    minx, miny, maxx, maxy = bbox
    (num_rows, num_cols) = data.shape
    
    readyFilename(out_filename, overwrite=True)
    
    drv = gdal.GetDriverByName(LMFormat.GTIFF.driver)
    ds = drv.Create(out_filename, num_cols, num_rows, 1, gdalconst.GDT_Byte)
    ds.SetGeoTransform([minx, cell_size, 0, maxy, 0, -cell_size])

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    ds.SetProjection(srs.ExportToWkt())
    
    outBand = ds.GetRasterBand(1)
    outBand.WriteArray(data)
    outBand.FlushCache()
    outBand.SetNoDataValue(nodata)

# .............................................................................
# Testing code 
Example 18
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getFormatLongName(self):
        name = ''
        drv = gdal.GetDriverByName(self._dataFormat)
        if drv is not None:
            name = drv.GetMetadataItem('DMD_LONGNAME')
        return name
            
# ............................................... 
Example 19
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def _copyGDALData(self, bandnum, infname, outfname, format='GTiff', kwargs={}):
        """
        @summary: Copy the dataset into a new file.  
        @param bandnum: The band number to read.
        @param outfname: Filename to write this dataset to.
        @param format: GDAL-writeable raster format to use for new dataset. 
                            http://www.gdal.org/formats_list.html
        @param doFlip: True if data begins at the southern edge of the region
        @param doShift: True if the leftmost edge of the data should be shifted 
                 to the center (and right half shifted around to the beginning) 
        @param nodata: Value used to indicate nodata in the new file.
        @param srs: Spatial reference system to use for the data. This is only 
                        necessary if the dataset does not have an SRS present.  This
                        will NOT project the dataset into a different projection.
        """
        options = []
        if format == 'AAIGrid':
            options = ['FORCE_CELLSIZE=True']
            #kwargs['FORCE_CELLSIZE'] = True
            #kwargs['DECIMAL_PRECISION'] = 4
        driver = gdal.GetDriverByName(format)
        metadata = driver.GetMetadata()
        if not (metadata.has_key(gdal.DCAP_CREATECOPY) 
                     and metadata[gdal.DCAP_CREATECOPY] == 'YES'):
            raise LMError(currargs='Driver %s does not support CreateCopy() method.' 
                              % format)
        inds = gdal.Open( infname )
        try:
            outds = driver.CreateCopy(outfname, inds, 0, options)
        except Exception, e:
            raise LMError(currargs='Creation failed for %s from band %d of %s (%s)'
                                          % (outfname, bandnum, infname, str(e))) 
Example 20
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def deleteData(self, dlocation=None, isTemp=False):
        """
        @summary: Deletes the local data file(s) on disk
        @note: Does NOT clear the dlocation attribute
        """
        success = False
        if dlocation is None:
            dlocation = self._dlocation
        if (dlocation is not None and os.path.isfile(dlocation)):
            drv = gdal.GetDriverByName(self._dataFormat)
            result = drv.Delete(dlocation)
            if result == 0:
                success = True
        if not isTemp:
            pth, fname = os.path.split(dlocation)
            if os.path.isdir(pth) and len(os.listdir(pth)) == 0:
                try:
                    os.rmdir(pth)
                except:
                    print 'Unable to rmdir %s' % pth
        return success

# .............................................................................
# Superclass _Layer methods overridden
# .............................................................................
# ............................................... 
Example 21
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def verifyField(self, dlocation, ogrFormat, attrname):
        """
        @summary: Read OGR-accessible data and save the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        """
        if attrname is None:
            fieldOk = True
        else:
            fieldOk = False
            if dlocation is not None and os.path.exists(dlocation):
                ogr.RegisterAll()
                drv = ogr.GetDriverByName(ogrFormat)
                try:
                    ds = drv.Open(dlocation)
                except Exception, e:
                    raise LMError(['Invalid datasource' % dlocation, str(e)])
                
                lyrDef = ds.GetLayer(0).GetLayerDefn()
                # Make sure given field exists and is the correct type
                for i in range(lyrDef.GetFieldCount()):
                    fld = lyrDef.GetFieldDefn(i)
                    fldname = fld.GetNameRef()
                    if attrname == fldname:
                        fldtype = fld.GetType()
                        if fldtype in (ogr.OFTInteger, ogr.OFTReal, ogr.OFTBinary):
                            fieldOk = True
                        break 
Example 22
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def readWithOGR(self, dlocation, ogrFormat, featureLimit=None, doReadData=False):
        """
        @summary: Read OGR-accessible data and set the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        @todo: remove featureLimit, read subsetDLocation if there is a limit 
        """
        thisBBox = localIdIdx = geomIdx = feats = featAttrs = None
        if dlocation is not None and os.path.exists(dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(ogrFormat)
            try:
                ds = drv.Open(dlocation)
            except Exception, e:
                raise LMError(['Invalid datasource' % dlocation, str(e)])
                          
            self.clearFeatures() 
            try:
                slyr = ds.GetLayer(0)
            except Exception, e:
                raise LMError(currargs='#### Failed to GetLayer from %s' % dlocation,
                                  prevargs=e.args, doTrace=True)

            # Get bounding box 
Example 23
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getSRS(self):
        if self._dlocation is not None and os.path.exists(self._dlocation):
            ogr.RegisterAll()
            drv = ogr.GetDriverByName(self._dataFormat)
            try:
                ds = drv.Open(self._dlocation)
            except Exception, e:
                raise LMError(['Invalid datasource' % self._dlocation, str(e)])

            vlyr = ds.GetLayer(0)
            srs = vlyr.GetSpatialRef()
            if srs is None:
                srs = self.createSRSFromEPSG()
            return srs 
Example 24
Project: core   Author: lifemapper   File: geotools.py    GNU General Public License v3.0 5 votes vote down vote up
def copyWithoutProjection(self, outfname, format='GTiff'):
      """
      @summary: Copy this dataset with no projection information.
      @param outfname: Filename to write this dataset to.
      @param format: GDAL-writeable raster format to use for new dataset. 
                     http://www.gdal.org/formats_list.html
      """
      self.openDataSet()
      driver = gdal.GetDriverByName(format)
      
      outds = driver.Create(outfname, self.xsize, self.ysize, 1, self.gdalBandType)
      if outds is None:
         print 'Creation failed for %s from band %d of %s' % (outfname, 1, 
                                                              self.dlocation)
         return 0 
      
      outds.SetGeoTransform(self.geoTransform)
      outband = outds.GetRasterBand(1)
      outband.SetNoDataValue(self.nodata)
      outds.SetProjection('')
      rstArray = self._dataset.ReadAsArray()
      outband.WriteArray( rstArray )

      # Close new dataset to flush to disk
      outds.FlushCache()
      outds = None

# ............................................. 
Example 25
Project: radiometric_normalization   Author: planetlabs   File: gimage.py    Apache License 2.0 5 votes vote down vote up
def create_ds(file_name, xsize, ysize, band_count, compress=True):
    options = ['PHOTOMETRIC=RGB']
    if compress:
        options.append('COMPRESS=DEFLATE')
        options.append('PREDICTOR=2')

    datatype = gdal.GDT_UInt16
    gdal_ds = gdal.GetDriverByName('GTIFF').Create(
        file_name, xsize, ysize, band_count, datatype,
        options=options)
    return gdal_ds 
Example 26
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def setUp(self):
        self.band = numpy.array([[0, 1], [2, 3]], dtype=numpy.uint16)
        self.mask = numpy.array([[0, 1], [0, 1]], dtype=numpy.bool)
        self.metadata = {'geotransform': (-1.0, 2.0, 0.0, 1.0, 0.0, -1.0)}

        self.test_photometric_alpha_image = 'test_photometric_alpha_image.tif'
        test_ds = gdal.GetDriverByName('GTiff').Create(
            self.test_photometric_alpha_image, 2, 2, 4, gdal.GDT_UInt16,
            options=['PHOTOMETRIC=RGB', 'ALPHA=YES'])
        gdal_array.BandWriteArray(test_ds.GetRasterBand(1), self.band)
        gdal_array.BandWriteArray(test_ds.GetRasterBand(2), self.band)
        gdal_array.BandWriteArray(test_ds.GetRasterBand(3), self.band)
        gdal_array.BandWriteArray(test_ds.GetRasterBand(4), self.mask)
        test_ds.SetGeoTransform(self.metadata['geotransform']) 
Example 27
Project: radiometric_normalization   Author: planetlabs   File: gimage_tests.py    Apache License 2.0 5 votes vote down vote up
def test__save_to_ds(self):
        output_file = 'test_save_to_ds.tif'

        test_band = numpy.array([[0, 1], [2, 3]], dtype=numpy.uint16)
        test_gimage = gimage.GImage([test_band], self.mask, self.metadata)
        output_ds = gdal.GetDriverByName('GTiff').Create(
            output_file, 2, 2, 2, gdal.GDT_UInt16,
            options=['ALPHA=YES'])
        gimage._save_to_ds(test_gimage, output_ds, nodata=3)

        # Required for gdal to write to file
        output_ds = None

        test_ds = gdal.Open(output_file)

        saved_number_of_bands = test_ds.RasterCount
        self.assertEquals(saved_number_of_bands, 2)

        saved_band = test_ds.GetRasterBand(1).ReadAsArray()
        numpy.testing.assert_array_equal(saved_band, self.band)

        saved_nodata = test_ds.GetRasterBand(1).GetNoDataValue()
        self.assertEqual(saved_nodata, 3)

        saved_alpha = test_ds.GetRasterBand(2).ReadAsArray()
        numpy.testing.assert_array_equal(saved_alpha, self.mask * 255)

        os.unlink(output_file) 
Example 28
Project: geophys_utils   Author: GeoscienceAustralia   File: _array2file.py    Apache License 2.0 5 votes vote down vote up
def array2file(data_arrays, 
               projection, 
               geotransform, 
               file_path, 
               file_format='GTiff', 
               dtype=gdal.GDT_Float32):
    '''
    Function to output array to any(?) GDAL supported raster format and return GDAL dataset
    Formats defined in http://www.gdal.org/formats_list.html
    '''
    data_array_shape = data_arrays[0].shape
    assert [data_array.shape for data_array in data_arrays].count(data_array_shape) == len(data_arrays), 'Data arrays are of different shape'

    driver = gdal.GetDriverByName(file_format)
    gdal_dataset = driver.Create(file_path, 
                                 data_array_shape[1], data_array_shape[0], # Array must be ordered yx
                                 len(data_arrays), 
                                 dtype)
    gdal_dataset.SetGeoTransform(geotransform)
    gdal_dataset.SetProjection(projection)
    
    # Write arrays to file
    for band_index in range(len(data_arrays)):
        raster_band = gdal_dataset.GetRasterBand(band_index+1)
        raster_band.WriteArray(data_arrays[band_index])
        raster_band.FlushCache()
        
    gdal_dataset.FlushCache()
    return gdal_dataset 
Example 29
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 30
Project: sentinelloader   Author: flaviostutz   File: utils.py    MIT License 5 votes vote down vote up
def saveGeoTiff(imageDataFloat, outputFile, geoTransform, projection):
    driver = gdal.GetDriverByName('GTiff')
    image_data = driver.Create(outputFile, imageDataFloat.shape[1], imageDataFloat.shape[0], 1, gdal.GDT_Float32)
    image_data.GetRasterBand(1).WriteArray(imageDataFloat)
    image_data.SetGeoTransform(geoTransform) 
    image_data.SetProjection(projection)
    image_data.FlushCache()
    image_data=None 
Example 31
Project: buzzard   Author: airware   File: _gdal_file_raster.py    Apache License 2.0 5 votes vote down vote up
def delete(self):
        super(BackGDALFileRaster, self).delete()

        success, payload = GDALErrorCatcher(gdal.GetDriverByName, none_is_error=True)(self.driver)
        if not success: # pragma: no cover
            raise ValueError('Could not find a driver named `{}` (gdal error: `{}`)'.format(
                self.driver, payload[1]
            ))
        dr = payload

        success, payload = GDALErrorCatcher(dr.Delete, nonzero_int_is_error=True)(self.path)
        if not success: # pragma: no cover
            raise RuntimeError('Could not delete `{}` using driver `{}` (gdal error: `{}`)'.format(
                self.path, dr.ShortName, payload[1]
            )) 
Example 32
Project: buzzard   Author: airware   File: test_dataset_registering.py    Apache License 2.0 5 votes vote down vote up
def random_path_tif():
    """Create a temporary path, and take care of cleanup afterward"""
    path = '{}/{}.tif'.format(tempfile.gettempdir(), uuid.uuid4())
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName('GTiff').Delete(path)
        except:
            os.remove(path) 
Example 33
Project: buzzard   Author: airware   File: test_dataset_registering.py    Apache License 2.0 5 votes vote down vote up
def random_path_shp():
    """Create a temporary path, and take care of cleanup afterward"""
    path = '{}/{}.shp'.format(tempfile.gettempdir(), uuid.uuid4())
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName('ESRI Shapefile').Delete(path)
        except:
            os.remove(path) 
Example 34
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def shp1_path(fps):
    """Create a shapefile in SR1 containing all single letter polygons from `fps` fixture"""
    path = '{}/{}.shp'.format(tempfile.gettempdir(), uuid.uuid4())

    ds = buzz.Dataset()
    ds.create_vector('poly', path, 'polygon', sr=SR1['wkt'])
    for letter in string.ascii_uppercase[:9]:
        ds.poly.insert_data(fps[letter].poly)
    ds.poly.close()
    del ds
    yield path
    gdal.GetDriverByName('ESRI Shapefile').Delete(path) 
Example 35
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def tif1_path(fps):
    """Create a tif in SR1 with all single letter footprints from `fps` fixture burnt in it"""
    path = '{}/{}.tif'.format(tempfile.gettempdir(), uuid.uuid4())

    ds = buzz.Dataset()
    ds.create_raster('rast', path, fps.AI, 'int32', 1, sr=SR1['wkt'])
    for letter in string.ascii_uppercase[:9]:
        fp = fps[letter]
        arr = np.full(fp.shape, ord(letter), dtype=int)
        ds.rast.set_data(arr, fp=fp)
    ds.rast.close()
    del ds
    yield path
    gdal.GetDriverByName('GTiff').Delete(path) 
Example 36
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def shp2_path(fps):
    """Create a shapefile in SR2 containing all single letter polygons from `fps` fixture"""
    path = '{}/{}.shp'.format(tempfile.gettempdir(), uuid.uuid4())

    ds = buzz.Dataset()
    ds.create_vector('poly', path, 'polygon', sr=SR2['wkt'])
    for letter in string.ascii_uppercase[:9]:
        ds.poly.insert_data(fps[letter].poly)
    ds.poly.close()
    del ds
    yield path
    gdal.GetDriverByName('ESRI Shapefile').Delete(path) 
Example 37
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def tif2_path(fps):
    """Create a tif in SR2 with all single letter footprints from `fps` fixture burnt in it"""
    path = '{}/{}.tif'.format(tempfile.gettempdir(), uuid.uuid4())

    ds = buzz.Dataset()
    with ds.acreate_raster(path, fps.AI, 'int32', 1, sr=SR2['wkt']).close as r:
        for letter in string.ascii_uppercase[:9]:
            fp = fps[letter]
            arr = np.full(fp.shape, ord(letter), dtype=int)
            r.set_data(arr, fp=fp)
    yield path
    gdal.GetDriverByName('GTiff').Delete(path) 
Example 38
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def shp3_path(fps):
    """Create a shapefile without SR containing all single letter polygons from `fps` fixture"""
    path = '{}/{}.shp'.format(tempfile.gettempdir(), uuid.uuid4())

    ds = buzz.Dataset()
    ds.create_vector('poly', path, 'polygon', sr=None)
    for letter in string.ascii_uppercase[:9]:
        ds.poly.insert_data(fps[letter].poly)
    ds.poly.close()
    del ds
    yield path
    gdal.GetDriverByName('ESRI Shapefile').Delete(path) 
Example 39
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def random_path_shp():
    """Create a temporary path, and take care of cleanup afterward"""
    path = '{}/{}.shp'.format(tempfile.gettempdir(), uuid.uuid4())
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName('ESRI Shapefile').Delete(path)
        except:
            os.remove(path) 
Example 40
Project: buzzard   Author: airware   File: test_dataset_modes.py    Apache License 2.0 5 votes vote down vote up
def random_path_tif():
    """Create a temporary path, and take care of cleanup afterward"""
    path = '{}/{}.tif'.format(tempfile.gettempdir(), uuid.uuid4())
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName('GTiff').Delete(path)
        except:
            os.remove(path)

# TESTS ***************************************************************************************** ** 
Example 41
Project: buzzard   Author: airware   File: test_rastersource_opencreate.py    Apache License 2.0 5 votes vote down vote up
def path(meta_file, ext):
    path = '{}/{}{}'.format(tempfile.gettempdir(), uuid.uuid4(), ext)
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName(meta['driver']).Delete(path)
        except:
            os.remove(path) 
Example 42
Project: buzzard   Author: airware   File: test_vectorsource_getsetdata_general.py    Apache License 2.0 5 votes vote down vote up
def path(suffix, driver):
    """Create a temporary path, and take care of cleanup afterward"""
    path = '{}/{}{}'.format(tempfile.gettempdir(), uuid.uuid4(), suffix)
    yield path
    if os.path.isfile(path):
        try:
            dr = gdal.GetDriverByName(driver)
            dr.Delete(path)
        except:
            os.remove(path) 
Example 43
Project: buzzard   Author: airware   File: test_vectorsource_opencreate.py    Apache License 2.0 5 votes vote down vote up
def path(driver_file):
    ext = EXTENSION_OF_DRIVER[driver_file]
    path = '{}/{}{}'.format(tempfile.gettempdir(), uuid.uuid4(), ext)
    yield path
    if os.path.isfile(path):
        try:
            gdal.GetDriverByName(driver_file).Delete(path)
        except:
            os.remove(path) 
Example 44
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_GDALIO.py    MIT License 5 votes vote down vote up
def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999):
    """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions.

    Args:
        FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written.
        newRasterfn (str): The filename (with path and extension) of the new raster.
        array (np.array): The array to be written
        driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format
        noDataValue (float): The no data value

    Return:
        np.array: A numpy array with the data from the raster.

    Author: SMM
    """

    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName(driver_name)
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outRaster.GetRasterBand(1).SetNoDataValue( noDataValue )
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
#============================================================================== 
Example 45
Project: eo_tools   Author: neel9102   File: reflectance_landsat_8_single_release.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array,geoTransform,proj,ncol,nrow):
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    outDataset = driver.Create(output_name,ncol,nrow,1,GDT_Float32)
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    outBand.FlushCache()
    outBand.SetNoDataValue(fillval)
    outDataset.SetGeoTransform(geoTransform )
    outDataset.SetProjection(proj) 
Example 46
Project: eo_tools   Author: neel9102   File: reflectance_landsat_8_multiple_release.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array,geoTransform,proj,ncol,nrow):
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    outDataset = driver.Create(output_name,ncol,nrow,1,GDT_Float32)
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    outBand.FlushCache()
    outBand.SetNoDataValue(fillval)
    outDataset.SetGeoTransform(geoTransform )
    outDataset.SetProjection(proj) 
Example 47
Project: eo_tools   Author: neel9102   File: reflectance_RE_multiple_release.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array,geoTransform,proj,ncol,nrow):
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    outDataset = driver.Create(output_name,ncol,nrow,1,GDT_Float32)
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    outBand.FlushCache()
    outBand.SetNoDataValue(fillval)
    outDataset.SetGeoTransform(geoTransform )
    outDataset.SetProjection(proj) 
Example 48
Project: eo_tools   Author: neel9102   File: mod09q1_lta_ndvi_computation_release.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array,geoTransform,proj,ncol,nrow):
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    outDataset = driver.Create(output_name,ncol,nrow,1,GDT_Float32)
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    outBand.FlushCache()
    outBand.SetNoDataValue(fillval)
    outDataset.SetGeoTransform(geoTransform )
    outDataset.SetProjection(proj) 
Example 49
Project: eo_tools   Author: neel9102   File: mod09q1_veg_index_proc_release.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array,geoTransform,proj):
    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    outDataset = driver.Create(output_name,ncol,nrow,1,GDT_Float32)
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    outBand.FlushCache()
    outBand.SetNoDataValue(fillval)
    outDataset.SetGeoTransform(geoTransform )
    outDataset.SetProjection(proj)



# finding all the hdf files in the mentioned root directories 
Example 50
Project: eo_tools   Author: neel9102   File: L2_LAC_SST_preprocessing.py    GNU General Public License v3.0 5 votes vote down vote up
def output_file(output_name,output_array):
  
    driver = gdal.GetDriverByName("GTiff")
    outDataset = driver.Create(output_name, ncols, nrows, 1, gdal.GDT_Float32)
    outDataset.SetGeoTransform(geo_transform_1)
    #srs = osr.SpatialReference()
    #srs.ImportFromEPSG(4326) #WGS84 lat long.
    outDataset.SetProjection(proj_1)
    
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(output_array,0,0)
    
    outBand.SetNoDataValue(np.nan)
    outBand.FlushCache()