Python osgeo.ogr.Open() Examples

The following are code examples for showing how to use osgeo.ogr.Open(). 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 testRaster(dlocation, bandnum=1):
        """
        @return: a GDAL dataset object
        """
        success = True
        try:
            try:
                dataset = gdal.Open(str(dlocation), gdalconst.GA_ReadOnly)
            except Exception, e:
                raise LMError(currargs='Unable to open dataset {} with GDAL ({})'
                                    .format(dlocation, str(e)))
            try:
                band = dataset.GetRasterBand(1)
            except Exception, e:
                raise LMError(currargs='No band {} in dataset {} ({})'
                                    .format(band, dlocation, str(e))) 
Example 2
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def isValidDataset(self, dlocation=None):
        """
        @summary: Checks to see if the dataset at self.dlocations is a valid 
                     vector readable by OGR.
        @return: True if dataset is a valid OGR dataset; False if not
        """
        valid = False
        if dlocation is None:
            dlocation = self._dlocation
        if (dlocation is not None
             and (os.path.isdir(dlocation) 
                    or os.path.isfile(dlocation)) ):
            try:
                ds = ogr.Open(dlocation)
                ds.GetLayer(0)
            except Exception, e:
                pass
            else:
                valid = True 
Example 3
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 4
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def copyData(self, sourceDataLocation, targetDataLocation=None, 
                     format=LMFormat.getDefaultOGR().driver):
        """
        Copy sourceDataLocation dataset to targetDataLocation or this layer's 
        dlocation.
        """
        if sourceDataLocation is not None and os.path.exists(sourceDataLocation):
            if targetDataLocation is None:
                if self._dlocation is not None:
                    targetDataLocation = self._dlocation
                else:
                    raise LMError('Target location is None')
        else:
            raise LMError('Source location %s is invalid' % str(sourceDataLocation))

        ogr.RegisterAll()
        drv = ogr.GetDriverByName(format)
        try:
            ds = drv.Open(sourceDataLocation)
        except Exception, e:
            raise LMError(['Invalid datasource' % sourceDataLocation, str(e)]) 
Example 5
Project: core   Author: lifemapper   File: extract_environment_values.py    GNU General Public License v3.0 6 votes vote down vote up
def get_metrics_for_layer(points, layer_filename, metricFunctions):
   """
   @summary: Get layer values for each point and then generate metrics
   """
   ds = gdal.Open(layer_filename)
   band = ds.GetRasterBand(1)
   data = np.array(band.ReadAsArray())
   gt = ds.GetGeoTransform()
   nodataVal = band.GetNoDataValue()
   
   values = []
   
   for x, y in points:
      px = int((x - gt[0]) / gt[1])
      py = int((y - gt[3]) / gt[5])
      try:
         val = data[py, px]
         if not is_close(val, nodataVal):
            values.append(data[py, px])
         else:
            print 'Could not append value at ({}, {}): {}'.format(px, py, val)
      except Exception, e:
            print 'Could not append value at ({}, {}): {}'.format(px, py, str(e)) 
Example 6
Project: core   Author: lifemapper   File: extract_environment_values.py    GNU General Public License v3.0 6 votes vote down vote up
def get_point_xys(points_filename, removeDuplicateLocations=True):
   """
   @summary: Get x,y pairs for each point in a shapefile
   """
   points = []
   
   ds = ogr.Open(points_filename)
   lyr = ds.GetLayer()
   
   for feat in lyr:
      geom = feat.GetGeometryRef()
      points.append((geom.GetX(), geom.GetY()))
   
   if removeDuplicateLocations:
      points = list(set(points))
      
   return points

# ............................................................................. 
Example 7
Project: Asistente-LADM_COL   Author: AgenciaImplementacion   File: dlg_import_from_excel.py    GNU General Public License v3.0 6 votes vote down vote up
def get_excel_info(self, path, sheetname):
        data_source = ogr.Open(path, 0)
        layer = data_source.GetLayerByName(sheetname)

        if layer is None:
            # A sheetname couldn't be found
            return None, None

        feature = layer.GetNextFeature()

        # If ogr recognizes the header, the first row will contain data, otherwise it'll contain field names
        header_in_first_row = True
        for field in self.fields[sheetname]:
            if feature.GetField(self.fields[sheetname].index(field)) != field:
                header_in_first_row = False

        num_rows = layer.GetFeatureCount()
        return header_in_first_row, num_rows - 1 if header_in_first_row else num_rows 
Example 8
Project: watools   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 6 votes vote down vote up
def Open_array_info(filename=''):
    """
    Opening a tiff info, for example size of array, projection and transform matrix.

    Keyword Arguments:
    filename -- 'C:/file/to/path/file.tif' or a gdal file (gdal.Open(filename))
        string that defines the input tiff file or gdal file

    """
    f = gdal.Open(r"%s" %filename)
    if f is None:
        print('%s does not exists' %filename)
    else:
        geo_out = f.GetGeoTransform()
        proj = f.GetProjection()
        size_X = f.RasterXSize
        size_Y = f.RasterYSize
        f = None
    return(geo_out, proj, size_X, size_Y) 
Example 9
Project: watools   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 6 votes vote down vote up
def Open_tiff_array(filename='', band=''):
    """
    Opening a tiff array.

    Keyword Arguments:
    filename -- 'C:/file/to/path/file.tif' or a gdal file (gdal.Open(filename))
        string that defines the input tiff file or gdal file
    band -- integer
        Defines the band of the tiff that must be opened.
    """
    f = gdal.Open(filename)
    if f is None:
        print('%s does not exists' %filename)
    else:
        if band is '':
            band = 1
        Data = f.GetRasterBand(band).ReadAsArray()
    return(Data) 
Example 10
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def strip_bands(in_raster_path, out_raster_path, bands_to_strip):
    in_raster = gdal.Open(in_raster_path)
    out_raster_band_count = in_raster.RasterCount-len(bands_to_strip)
    out_raster = create_matching_dataset(in_raster, out_raster_path, bands=out_raster_band_count)
    out_raster_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    in_raster_array = in_raster.GetVirtualMemArray()

    bands_to_copy = [band for band in range(in_raster_array.shape[0]) if band not in bands_to_strip]

    out_raster_array[...] = in_raster_array[bands_to_copy, :,:]

    out_raster_array = None
    in_raster_array = None
    out_raster = None
    in_raster = None

    return out_raster_path 
Example 11
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def flatten_probability_image(prob_image, out_path):
    """
    Takes a probability output from classify_image and flattens it into a single layer containing only the maximum
    value from each pixel.

    Parameters
    ----------
    prob_image
        The path to a probability image.
    out_path
        The place to save the flattened image.

    """
    prob_raster = gdal.Open(prob_image)
    out_raster = create_matching_dataset(prob_raster, out_path, bands=1)
    prob_array = prob_raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    out_array[:, :] = prob_array.max(axis=0)
    out_array = None
    prob_array = None
    out_raster = None
    prob_raster = None 
Example 12
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def raster_to_array(rst_pth):
    """Reads in a raster file and returns a N-dimensional array.

    Parameters
    ----------
    rst_pth
        Path to input raster.
    Returns
    -------
    As N-dimensional array.
    """
    log = logging.getLogger(__name__)
    in_ds = gdal.Open(rst_pth)
    out_array = in_ds.ReadAsArray()

    return out_array 
Example 13
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def calc_ndvi(raster_path, output_path):
    import pdb
    raster = gdal.Open(raster_path)
    out_raster = create_matching_dataset(raster, output_path, datatype=gdal.GDT_Float32)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    R = array[2, ...]
    I = array[3, ...]
    out_array[...] = (R-I)/(R+I)
    pdb.set_trace()

    out_array[...] = np.where(out_array == -2147483648, 0, out_array)

    R = None
    I = None
    array = None
    out_array = None
    raster = None
    out_raster = None 
Example 14
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_model(image_path, model_path, model_clear=0, num_chunks=10, buffer_size=0):
    """Returns a multiplicative mask (0 for cloud, shadow or haze, 1 for clear) built from the model at model_path."""
    from pyeo.classification import classify_image  # Deferred import to deal with circular reference
    with TemporaryDirectory() as td:
        log = logging.getLogger(__name__)
        log.info("Building cloud mask for {} with model {}".format(image_path, model_path))
        temp_mask_path = os.path.join(td, "cat_mask.tif")
        classify_image(image_path, model_path, temp_mask_path, num_chunks=num_chunks)
        temp_mask = gdal.Open(temp_mask_path, gdal.GA_Update)
        temp_mask_array = temp_mask.GetVirtualMemArray()
        mask_path = get_mask_path(image_path)
        mask = create_matching_dataset(temp_mask, mask_path, datatype=gdal.GDT_Byte)
        mask_array = mask.GetVirtualMemArray(eAccess=gdal.GF_Write)
        mask_array[:, :] = np.where(temp_mask_array != model_clear, 0, 1)
        temp_mask_array = None
        mask_array = None
        temp_mask = None
        mask = None
        if buffer_size:
            buffer_mask_in_place(mask_path, buffer_size)
        log.info("Cloud mask for {} saved in {}".format(image_path, mask_path))
        return mask_path 
Example 15
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_class_map(class_map_path, out_path, classes_of_interest, buffer_size=0, out_resolution=None):
    """Creates a mask from a classification mask: 1 for each pixel containing one of classes_of_interest, otherwise 0"""
    # TODO: pull this out of the above function
    class_image = gdal.Open(class_map_path)
    class_array = class_image.GetVirtualMemArray()
    mask_array = np.isin(class_array, classes_of_interest)
    out_mask = create_matching_dataset(class_image, out_path, datatype=gdal.GDT_Byte)
    out_array = out_mask.GetVirtualMemArray(eAccess=gdal.GA_Update)
    np.copyto(out_array, mask_array)
    class_array = None
    class_image = None
    out_array = None
    out_mask = None
    if out_resolution:
        resample_image_in_place(out_path, out_resolution)
    if buffer_size:
        buffer_mask_in_place(out_path)
    return out_path 
Example 16
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 6 votes vote down vote up
def create_mask_from_fmask(in_l1_dir, out_path):
    log = logging.getLogger(__name__)
    log.info("Creating fmask for {}".format(in_l1_dir))
    with TemporaryDirectory() as td:
        temp_fmask_path = os.path.join(td, "fmask.tif")
        apply_fmask(in_l1_dir, temp_fmask_path)
        fmask_image = gdal.Open(temp_fmask_path)
        fmask_array = fmask_image.GetVirtualMemArray()
        out_image = create_matching_dataset(fmask_image, out_path, datatype=gdal.GDT_Byte)
        out_array = out_image.GetVirtualMemArray(eAccess=gdal.GA_Update)
        log.info("fmask created, converting to binary cloud/shadow mask")
        out_array[:,:] = np.isin(fmask_array, (2, 3, 4), invert=True)
        out_array = None
        out_image = None
        fmask_array = None
        fmask_image = None
        resample_image_in_place(out_path, 10) 
Example 17
Project: swat-tools   Author: DigitalAgricultureDiscovery   File: process.py    GNU General Public License v3.0 6 votes vote down vote up
def update_hru_shapefile(self, inshape, new_data):
        # open hru1.shp
        src = ogr.Open(inshape, 1)

        # get the hru1 layer
        lyr = src.GetLayer()

        # create new field we wish to append (Runoff or Sediment)
        lyr.CreateField(
            ogr.FieldDefn(str(self.fieldswat_output_type).title(), ogr.OFTReal))

        feature = lyr.GetNextFeature()

        for i in range(0, len(new_data)):
            feature.SetField(str(self.fieldswat_output_type), new_data[i])
            lyr.SetFeature(feature)
            feature = lyr.GetNextFeature() 
Example 18
Project: HackZurich   Author: RefugeeMatchmaking   File: test_shp.py    MIT License 6 votes vote down vote up
def test_attributeexport(self):
        def testattributes(lyr, graph):
            feature = lyr.GetNextFeature()
            while feature:
                coords = []
                ref = feature.GetGeometryRef()
                last = ref.GetPointCount() - 1
                edge_nodes = (ref.GetPoint_2D(0), ref.GetPoint_2D(last))
                name = feature.GetFieldAsString('Name')
                assert_equal(graph.get_edge_data(*edge_nodes)['Name'], name)
                feature = lyr.GetNextFeature()

        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')

        G = nx.read_shp(self.shppath)
        nx.write_shp(G, tpath)
        shpdir = ogr.Open(tpath)
        edges = shpdir.GetLayerByName("edges")
        testattributes(edges, G) 
Example 19
Project: HackZurich   Author: RefugeeMatchmaking   File: test_shp.py    MIT License 6 votes vote down vote up
def test_wkt_export(self):
        G = nx.DiGraph()
        tpath = os.path.join(tempfile.gettempdir(), 'shpdir')
        points = (
            "POINT (0.9 0.9)",
            "POINT (4 2)"
        )
        line = (
            "LINESTRING (0.9 0.9,4 2)",
        )
        G.add_node(1, Wkt=points[0])
        G.add_node(2, Wkt=points[1])
        G.add_edge(1, 2, Wkt=line[0])
        try:
            nx.write_shp(G, tpath)
        except Exception as e:
            assert False, e
        shpdir = ogr.Open(tpath)
        self.checkgeom(shpdir.GetLayerByName("nodes"), points)
        self.checkgeom(shpdir.GetLayerByName("edges"), line) 
Example 20
Project: CensusShapefileMaker   Author: GreenInfo-Network   File: GenerateStateMapperData.py    MIT License 6 votes vote down vote up
def main(self):
        print "    Assigning ACS fields to records in censustracts.shp"

        # 1 - open the shapefile and add attributes to it (unless it exists already)
        source = ogr.Open('censusblocks.shp', 1)
        layer  = source.GetLayer()
        defn   = layer.GetLayerDefn()

        if defn.GetFieldIndex('MHHINC') == -1:
            layer.CreateField( ogr.FieldDefn('MHHINC', ogr.OFTInteger) )

        # 2 - then iterate over all records in it, and populate them from the ACS data we cached during init
        # the tract-ID is the first 11 characters of the GEOID, we can key from that
        feature = layer.GetNextFeature()
        while feature:
            tractid    = feature.GetField("GEOID")[:11]
            attributes = self.tract_attributes[tractid]
            feature.SetField("MHHINC", attributes['MHHINC'] )
            layer.SetFeature(feature)
            feature = layer.GetNextFeature()

        # 999 - done
        source.Destroy() 
Example 21
Project: qgis-safecast-plugin   Author: OpenGeoLabs   File: radiation_toolbox_dockwidget.py    GNU General Public License v2.0 6 votes vote down vote up
def _getLayerType(self, layer):
        """Check if current map layer is managable by plugin.

        :return: layer type or None (unknown)
        """
        if not layer:
            return None

        # first check whether layer has been loaded by the plugin
        if hasattr(layer, "layerType"):
            return layer.layerType

        # then check layer metadata
        try:
            filePath = layer.dataProvider().dataSourceUri().split('|')[0]
            ds = ogr.Open(filePath)
            layer = ds.GetLayerByName('safecast_metadata') is not None
            ds = None
            if layer:
                return LayerType.Safecast
        except:
            return None

        return None 
Example 22
Project: qgis-safecast-plugin   Author: OpenGeoLabs   File: safecast.py    GNU General Public License v2.0 6 votes vote down vote up
def _getMetadata(self):
        """Get metadata."""
        if isinstance(self._layer, SafecastLayer):
            return self._layer.metadata

        metadata = {}
        try:
            ds = ogr.Open(self._fileName)
            layer = ds.GetLayerByName('safecast_metadata')
            if layer:
                layer_defn = layer.GetLayerDefn()
                layer.ResetReading()
                feat = layer.GetNextFeature()
                for i in range(layer_defn.GetFieldCount()):
                    name = layer_defn.GetFieldDefn(i).GetName()
                    value = feat.GetField(i)
                    metadata[name] = value
            ds = None
        except:
            raise LoadError(
                self.tr("Unable to retrive Safecast metadata for selected layer")
            )

        return metadata 
Example 23
Project: qgis-wps4server   Author: 3liz   File: UMN.py    GNU General Public License v3.0 6 votes vote down vote up
def getDataset(self, output):
        """
        :param output: :class:`pywps.Process.InAndOutputs.ComplexOutput`
        :returns: "raster" or "vector"
        """

        logging.debug("Importing given output [%s] using gdal" % output.value)
        # If dataset is XML it will make an error like ERROR 4:
        # `/var/www/html/wpsoutputs/vectorout-26317EUFxeb' not recognised as a
        # supported file format.
        self.dataset = gdal.Open(output.value)

        if self.dataset:
            return "raster"

        if not self.dataset:
            logging.debug(
                "Importing given output [%s] using ogr" % output.value)
            self.dataset = ogr.Open(output.value)

        if self.dataset:
            return "vector"
        else:
            return None 
Example 24
Project: core   Author: lifemapper   File: create_mask.py    GNU General Public License v3.0 5 votes vote down vote up
def get_layer_dimensions(template_layer_filename):
    """
    @summary: Get layer information
    @param template_layer_filename: The template layer to get information for
    """
    ds = gdal.Open(template_layer_filename)

    num_cols = ds.RasterXSize
    num_rows = ds.RasterYSize
    
    minx, cell_size, _, maxy, _, _ = ds.GetGeoTransform()
    
    maxx = minx + (cell_size * num_cols)
    miny = maxy - (cell_size * num_rows)
    
    bbox = (minx, miny, maxx, maxy)

    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjectionRef())
    # TODO: Try to find a better way to fail
    try:
        epsg = int(srs.GetAttrValue('AUTHORITY', 1))
    except:
        epsg = DEFAULT_EPSG
    ds = None
    
    return bbox, cell_size, epsg

# ............................................................................. 
Example 25
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def _openWithGDAL(self, dlocation=None, bandnum=1):
        """
        @return: a GDAL dataset object
        """
        if dlocation is None:
            dlocation = self._dlocation
        try:
            dataset = gdal.Open(str(dlocation), gdalconst.GA_ReadOnly)
            band = dataset.GetRasterBand(1)
        except Exception, e:
            raise LMError(['Unable to open dataset or band {} with GDAL ({})'
                                .format(dlocation, str(e))]) 
Example 26
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 27
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)):
            ds = gdal.Open(str(self._dlocation), gdalconst.GA_ReadOnly)
            wktSRS = ds.GetProjection()
            if wktSRS is not '':
                srs = osr.SpatialReference()
                srs.ImportFromWkt(wktSRS)
            else:
                srs = self.createSRSFromEPSG()
            ds = None
            return srs
        else:
            raise LMError(currargs='Input file %s does not exist' % self._dlocation)

# ............................................... 
Example 28
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 29
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 30
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 31
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def _read_raster_layer(self, raster_filename):
        """Reads a raster layer for processing

        Args:
            raster_filename: The file path for the raster layer.

        Returns:
            A tuple containing a window function for returning a portion of the
            numpy array generated by the layer and the NODATA value to use with
            this layer.
        """
        ds = gdal.Open(raster_filename)
        band = ds.GetRasterBand(1)
        layer_array = band.ReadAsArray()
        nodata = band.GetNoDataValue()
        
        num_y, num_x = layer_array.shape
        min_x, x_res, _, max_y, _, y_res = ds.GetGeoTransform()
        max_x = min_x + (num_x * x_res)
        min_y = max_y + (y_res * num_y)
        layer_bbox = (min_x, min_y, max_x, max_y)
        window_func = self._get_window_function(
            layer_array, layer_bbox, self.shapegrid_resolution,
            num_cell_sides=self.shapegrid_sides)
        
        return (window_func, nodata)
    
    # ............................... 
Example 32
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def _read_shapegrid(self, shapegrid_filename):
        """Read the shapegrid

        Args:
            shapegrid_filename: The file location of the shapegrid
        """
        shapegrid_dataset = ogr.Open(shapegrid_filename)
        self.shapegrid_layer = shapegrid_dataset.GetLayer()
        tmp = self.shapegrid_layer.GetExtent()
        self.shapegrid_bbox = (tmp[0], tmp[2], tmp[1], tmp[3])
        
        self.num_cells = self.shapegrid_layer.GetFeatureCount()

        feature_0 = self.shapegrid_layer.GetFeature(0)
        geom = feature_0.GetGeometryRef()
        
        # Get resolution and number of sides
        geom_wkt = geom.ExportToWkt()
        boundary_points = geom_wkt.split(',')
        if len(boundary_points) == 5:
            # Square
            envelope = geom.GetEnvelope()
            self.shapegrid_resolution = (envelope[1] - envelope[0],
                                         envelope[3] - envelope[2])
        else:
            # Hexagon
            center = geom.Centroid()
            x_cent = center.GetX()
            y_cent = center.GetY()
            x1, y1 = boundary_points[1].split(' ')
            self.shapegrid_resolution = np.sqrt(
                (x_cent - x1)**2 + (y_cent - y1)**2)
        self.shapegrid_sides = len(boundary_points) - 1
        #self.shapegrid_layer.ResetReading()
        self.shapegrid_layer = None

    # ............................... 
Example 33
Project: core   Author: lifemapper   File: layer_encoder.py    GNU General Public License v3.0 5 votes vote down vote up
def get_geojson(self):
        """Formats the encoded matrix as GeoJSON
        """
        ret = {
            'type' : 'FeatureCollection'
        }
        features = []
        
        column_headers = self.encoded_matrix.getColumnHeaders()
        
        column_enum = [(j, str(k)) for j, k in enumerate(column_headers)]
        
        shapegrid_dataset = ogr.Open(self.shapegrid_filename)
        shapegrid_layer = shapegrid_dataset.GetLayer()
        
        i = 0
        feat = shapegrid_layer.GetNextFeature()
        while feat is not None:
            ft_json = json.loads(feat.ExportToJson())
            # right_hand_rule(ft_json['geometry']['coordinates'])
            # TODO(CJ): Remove this if updated library adds first id correctly
            ft_json['id'] = feat.GetFID()
            ft_json['properties'] = dict(
                [(k, self.encoded_matrix.data[i, j].item()
                  ) for j, k in column_enum])
            features.append(ft_json)
            i+=1
            feat = shapegrid_layer.GetNextFeature()

        ret['features'] = features
        shapegrid_dataset = shapegrid_layer = None
        return ret 
Example 34
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def __init__(self, name=None, long_name=None, shapefile=None, folder=None, data_source=None, 
               load=False, ldebug=False, shapetype=None):
    ''' load shapefile '''
    if name is not None and not isinstance(name,str): raise TypeError
    if folder is not None and not isinstance(folder,str): raise TypeError
    if shapefile is not None and not isinstance(shapefile,str): raise TypeError
    # resolve file name and open file
    if ldebug: print(' - loading shapefile')
    if shapefile is None: shapefile = name + '.shp' # will be absolute path
    elif shapefile[-4:] != '.shp': shapefile = shapefile + '.shp' # and need extension
    if folder is not None: shapefile = folder + '/' + shapefile
    if name is None: name = os.path.splitext(os.path.basename(shapefile))[0]
    if long_name is None: long_name = name
    else: folder = os.path.dirname(shapefile)
    if not os.path.exists(shapefile): raise IOError('File \'{}\' not found!'.format(shapefile))
    if ldebug: print((' - using shapefile \'{}\''.format(shapefile)))
    # instance attributes
    self.name = name
    self.long_name = long_name
    self.folder = folder
    self.shapefile = shapefile
    self.data_source = data_source # source documentation...
    shapetype = self.__class__.__name__ if shapetype is None else shapetype
    self.shapetype = shapetype # for specific types of shapes, e.g. Basin, Lake, Prov, Natl
    # load shapefile (or not)
    self._ogr = ogr.Open(shapefile) if load else None 
Example 35
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def OGR(self):
    ''' access to OGR dataset '''
    if self._ogr is None: self._ogr = ogr.Open(self.shapefile) # load data, if not already done 
    return self._ogr 
Example 36
Project: GeoPy   Author: aerler   File: gdal.py    GNU General Public License v3.0 5 votes vote down vote up
def load(self):
    ''' load data and return self '''
    if self._ogr is None: self._ogr = ogr.Open(self.shapefile) # load data, if not already done 
    return self 
Example 37
Project: Data-Tools-for-ArcGIS   Author: datastory   File: Service Geometry Downloader.py    MIT License 5 votes vote down vote up
def _get_first_source_dataset(folder):
    """
    :return: osgeo.ogr.DataSource
    """

    first_source_file = os.listdir(folder)[0]
    first_source_ds = ogr.Open(os.path.join(folder, first_source_file))
    return first_source_ds 
Example 38
Project: Data-Tools-for-ArcGIS   Author: datastory   File: Service Geometry Downloader.py    MIT License 5 votes vote down vote up
def join_files(folder):
    """
    Merge all files in given folder into one datasource

    :param folder: name of temporary folder to be joined
    :return: resulting joined files file name
    """

    first_ds = _get_first_source_dataset(folder)
    first_layer = first_ds.GetLayer(0)

    drv = ogr.GetDriverByName('GeoJSON')
    tmpfile = tempfile.mktemp(suffix=".json")
    out_dsn = drv.CreateDataSource(tmpfile)
    out_layer = out_dsn.CopyLayer(first_layer, new_name=first_layer.GetName())

    for source in os.listdir(folder)[1:]:
        logging.info('Joining file %s to %s' % (source, tmpfile)) 
        dsn = ogr.Open(os.path.join(folder, source))
        layer = dsn.GetLayer()
        nfeatures = layer.GetFeatureCount()

        for i in range(nfeatures):
            feature = layer.GetNextFeature()
            out_layer.CreateFeature(feature.Clone())

    out_dsn.Destroy()

    return tmpfile 
Example 39
Project: pyq-osrm   Author: mthh   File: pyq-osrm.py    MIT License 5 votes vote down vote up
def read_shp(file_path):
    my_dict = {}
    coord_liste = []
    datasource = ogr.Open(file_path)
    layer = datasource.GetLayer(0)
    if layer.GetGeomType() != 1:
        sys.exit('Input shapefile must be a Points layer\n')
    for feature in layer:
        concat = str(feature.geometry().GetY()) \
            + ',' + str(feature.geometry().GetX())
        coord_liste.append(concat)
        my_dict[concat] = feature.GetField(0)

    return my_dict, coord_liste 
Example 40
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_VectorTools.py    MIT License 5 votes vote down vote up
def readFile(filename):
	print("Hey buddy, Reading the file: "+filename)

	filehandle = gdal.Open(filename, GA_ReadOnly )
	if filehandle == None:
		raise Exception("Unable to read the data file")

	band1 = filehandle.GetRasterBand(1)
	geotransform = filehandle.GetGeoTransform()
	geoproj = filehandle.GetProjection()
	Z = band1.ReadAsArray()
	xsize = filehandle.RasterXSize
	ysize = filehandle.RasterYSize
	return xsize,ysize,geotransform,geoproj,Z 
Example 41
Project: LSDMappingTools   Author: LSDtopotools   File: LSD_GeologyTools.py    MIT License 5 votes vote down vote up
def readFile(filename):
    print("Hey buddy, Reading the file: "+filename)

    filehandle = gdal.Open(filename, GA_ReadOnly )
    if filehandle == None:
        raise Exception("Unable to read the data file")

    band1 = filehandle.GetRasterBand(1)
    geotransform = filehandle.GetGeoTransform()
    geoproj = filehandle.GetProjection()
    Z = band1.ReadAsArray()
    xsize = filehandle.RasterXSize
    ysize = filehandle.RasterYSize
    return xsize,ysize,geotransform,geoproj,Z 
Example 42
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def ogrWkt2Shapely(input_shape):
    # this throws away the other attributes of the feature, but is 
    # sufficient in this use case
    shapely_objects=[]
    shp = ogr.Open(input_shape)
    lyr = shp.GetLayer()
    for n in range(0, lyr.GetFeatureCount()):
        feat = lyr.GetFeature(n)
        wkt_feat = loads(feat.geometry().ExportToWkt())
        shapely_objects.append(wkt_feat)
    return shapely_objects 
Example 43
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def log_file_fields(filename):
    print("File: " + filename)
    source = ogr.Open(filename)

    for i in range(source.GetLayerCount()):
        layer = source.GetLayerByIndex(i)
        layerName = layer.GetName()
        print("Layer: " + layerName)
        stringFields = []
        layerDefinition = layer.GetLayerDefn()
        for n in range(layerDefinition.GetFieldCount()):
            fieldDefinition = layerDefinition.GetFieldDefn(n)
            fieldName =  fieldDefinition.GetName()
            fieldTypeCode = fieldDefinition.GetType()
            fieldType = fieldDefinition.GetFieldTypeName(fieldTypeCode)
            if fieldType == "String":
                stringFields.append(fieldName)
    
        print("String Fields:")
        for field in stringFields:
            sql = 'SELECT %s FROM %s' % (field, layerName)
            fieldLayer = source.ExecuteSQL(sql)
            values = {}
            for i, feature in enumerate(fieldLayer):
                values[feature.GetField(0)] = values.get(feature.GetField(0), 0) + 1
        
            print("Field: " + field)
            for key in sorted(values.keys()):
                print("'%s': %d" % (key, values[key]))
            print("\n") 
Example 44
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 5 votes vote down vote up
def ogrWkt2Shapely(input_shape):
    # this throws away the other attributes of the feature, but is 
    # sufficient in this use case
    shapely_objects=[]
    shp = ogr.Open(input_shape)
    lyr = shp.GetLayer()
    for n in range(0, lyr.GetFeatureCount()):
        feat = lyr.GetFeature(n)
        wkt_feat = loads(feat.geometry().ExportToWkt())
        shapely_objects.append(wkt_feat)
    return shapely_objects 
Example 45
Project: eo_tools   Author: neel9102   File: MAGIC_clip_raster_by_poly_release.py    GNU General Public License v3.0 5 votes vote down vote up
def coord2pixelOffset(rasterfn,x,y):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    xOffset = int((x - originX)/pixelWidth)
    yOffset = int((y - originY)/pixelHeight)
    return xOffset,yOffset 
Example 46
Project: eo_tools   Author: neel9102   File: MAGIC_clip_raster_by_poly_release.py    GNU General Public License v3.0 5 votes vote down vote up
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(fillvalue)

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds 
Example 47
Project: watools   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 5 votes vote down vote up
def Open_bil_array(bil_filename, band = 1):
    """
    Opening a bil array.

    Keyword Arguments:
    bil_filename -- 'C:/file/to/path/file.bil'
        string that defines the input tiff file or gdal file
    band -- integer
        Defines the band of the tiff that must be opened.
    """    
    gdal.GetDriverByName('EHdr').Register()
    img = gdal.Open(bil_filename)
    Data = img.GetRasterBand(band).ReadAsArray()
    
    return(Data) 
Example 48
Project: watools   Author: wateraccounting   File: raster_conversions.py    Apache License 2.0 5 votes vote down vote up
def clip_data(input_file, latlim, lonlim):
    """
    Clip the data to the defined extend of the user (latlim, lonlim) or to the
    extend of the DEM tile

    Keyword Arguments:
    input_file -- output data, output of the clipped dataset
    latlim -- [ymin, ymax]
    lonlim -- [xmin, xmax]
    """
    try:
        if input_file.split('.')[-1] == 'tif':
            dest_in = gdal.Open(input_file)
        else:
            dest_in = input_file
    except:
        dest_in = input_file

    # Open Array
    data_in = dest_in.GetRasterBand(1).ReadAsArray()

    # Define the array that must remain
    Geo_in = dest_in.GetGeoTransform()
    Geo_in = list(Geo_in)
    Start_x = np.max([int(np.floor(((lonlim[0]) - Geo_in[0])/ Geo_in[1])),0])
    End_x = np.min([int(np.ceil(((lonlim[1]) - Geo_in[0])/ Geo_in[1])),int(dest_in.RasterXSize)])

    Start_y = np.max([int(np.floor((Geo_in[3] - latlim[1])/ -Geo_in[5])),0])
    End_y = np.min([int(np.ceil(((latlim[0]) - Geo_in[3])/Geo_in[5])), int(dest_in.RasterYSize)])

    #Create new GeoTransform
    Geo_in[0] = Geo_in[0] + Start_x * Geo_in[1]
    Geo_in[3] = Geo_in[3] + Start_y * Geo_in[5]
    Geo_out = tuple(Geo_in)

    data = np.zeros([End_y - Start_y, End_x - Start_x])

    data = data_in[Start_y:End_y,Start_x:End_x]
    dest_in = None

    return(data, Geo_out) 
Example 49
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def reproject_image(in_raster, out_raster_path, new_projection,  driver = "GTiff",  memory = 2e3, do_post_resample=True):
    """
    Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function.

    Parameters
    ----------
    in_raster
        Either a gdal.Raster object or a path to a raster
    out_raster_path
        The path to the new output raster.
    new_projection
        The new projection in .wkt
    driver
        The format of the output raster.
    memory
        The amount of memory to give to the reprojection.
    do_post_resample
        If set to false, do not resample the image back to the original projection.

    Notes
    -----
    The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image
    can become a 10.002m pixel resolution image. To stop alignment issues,
    by deafult this function resamples the images back to their original resolution

    """
    log = logging.getLogger(__name__)
    log.info("Reprojecting {} to {}".format(in_raster, new_projection))
    if type(in_raster) is str:
        in_raster = gdal.Open(in_raster)
    res = in_raster.GetGeoTransform()[1]
    gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver)
    # After warping, image has irregular gt; resample back to previous pixel size
    # TODO: Make this an option
    if do_post_resample:
        resample_image_in_place(out_raster_path, res)
    return out_raster_path 
Example 50
Project: pyeo   Author: clcr   File: raster_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def apply_band_function(in_path, function, bands, out_path, out_datatype = gdal.GDT_Int32):
    """Applys an arbitrary band mathemtics function to an image at in_path and saves the result at out_map.
    Function should be a function ofblect of the form f(band_input_A, band_input_B, ...)"""
    raster = gdal.Open(in_path)
    out_raster = create_matching_dataset(raster, out_path=out_path, datatype=out_datatype)
    array = raster.GetVirtualMemArray()
    out_array = out_raster.GetVirtualMemArray(eAccess=gdal.GA_Update)
    band_views = [array[band, ...] for band in bands]
    out_array[...] = function(*band_views)
    out_array = None
    for view in band_views:
        view = None
    raster = None
    out_raster = None