Python osgeo.ogr.Feature() Examples

The following are code examples for showing how to use osgeo.ogr.Feature(). 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: unmixing   Author: arthur-e   File: lsma.py    MIT License 6 votes vote down vote up
def get_idx_as_shp(self, path, gt, wkt):
        '''
        Exports a Shapefile containing the locations of the extracted
        endmembers. Assumes the coordinates are in decimal degrees.
        '''
        coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)

        driver = ogr.GetDriverByName('ESRI Shapefile')
        ds = driver.CreateDataSource(path)
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)

        layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
        for pair in coords:
            feature = ogr.Feature(layer.GetLayerDefn())

            # Create the point from the Well Known Text
            point = ogr.CreateGeometryFromWkt('POINT(%f %f)' %  pair)
            feature.SetGeometry(point) # Set the feature geometry
            layer.CreateFeature(feature) # Create the feature in the layer
            feature.Destroy() # Destroy the feature to free resources

        # Destroy the data source to free resources
        ds.Destroy() 
Example 2
Project: surfclass   Author: Kortforsyningen   File: vectorize.py    MIT License 6 votes vote down vote up
def process(self):
        """Start processing."""
        logger.debug("Started processing")
        self._add_fields()
        vdefn = self._outlyr.GetLayerDefn()
        for f in self._featurereader:
            all_classes = dict(self._zero_counts)
            class_counts = self._count_classes_inside(f.geometry())
            # Add zero counts for classes not seen
            all_classes.update(class_counts)
            outfeat = ogr.Feature(vdefn)
            outfeat.SetFrom(f)
            total_count = 0
            for class_id, count in all_classes.items():
                total_count += count
                try:
                    field_id = self._class_field_index[class_id]
                    outfeat.SetFieldInteger64(field_id, count)
                except KeyError:
                    pass
            outfeat.SetFieldInteger64(self._total_field_index, total_count)
            self._outlyr.CreateFeature(outfeat) 
Example 3
Project: PyGeoC   Author: lreis2415   File: vector.py    MIT License 6 votes vote down vote up
def write_line_shp(line_list, out_shp):
        """Export ESRI Shapefile -- Line feature"""
        print('Write line shapefile: %s' % out_shp)
        driver = ogr_GetDriverByName(str('ESRI Shapefile'))
        if driver is None:
            print('ESRI Shapefile driver not available.')
            sys.exit(1)
        if os.path.exists(out_shp):
            driver.DeleteDataSource(out_shp)
        ds = driver.CreateDataSource(out_shp.rpartition(os.sep)[0])
        if ds is None:
            print('ERROR Output: Creation of output file failed.')
            sys.exit(1)
        lyr = ds.CreateLayer(str(out_shp.rpartition(os.sep)[2].split('.')[0]), None, wkbLineString)
        for l in line_list:
            line = ogr_Geometry(wkbLineString)
            for i in l:
                line.AddPoint(i[0], i[1])
            templine = ogr_CreateGeometryFromJson(line.ExportToJson())
            feature = ogr_Feature(lyr.GetLayerDefn())
            feature.SetGeometry(templine)
            lyr.CreateFeature(feature)
            feature.Destroy()
        ds.Destroy() 
Example 4
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 6 votes vote down vote up
def setDamHeights(self, feat, low, mid, high):
        """
        Set dam heights to be modeled for a dam.

        :param feat: OGR Feature of dam.
        :param low: Lower interval of dam height.
        :param mid: Median dam height.
        :param high: Upper interval of dam height.

        :return: Updated OGR Feature of dam.
        """
        feat.SetField("ht_lo", low)
        feat.SetField("ht_mid", mid)
        feat.SetField("ht_hi", high)
        feat.SetField("ht_lo_mod", low)
        feat.SetField("ht_mid_mod", mid)
        feat.SetField("ht_hi_mod", high)
        #feat.SetField("ht_max", max)
        return feat 
Example 5
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def features(self):
        """
        @summary: Converts the private dictionary of features into a list of 
                         LmAttObjs
        @note: Uses list comprehensions to create a list of LmAttObjs and another
                     to create a list of (key, value) pairs for the attribute 
                     dictionary
        @return: A list of LmAttObjs
        """
        return [LmAttObj(dict([
                            (self._featureAttributes[k2][0], self._features[k1][k2]) \
                            for k2 in self._featureAttributes]), 
                              "Feature") for k1 in self._features]
    
    # .................................. 
Example 6
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def getFeatureValByFieldIndex(self, fieldIdx, featureFID):
        if self._features:
            if self._features.has_key(featureFID):
                return self._features[featureFID][fieldIdx]
            else:
                raise LMError ('Feature ID %s not found in dataset %s' % 
                                        (str(featureFID), self._dlocation))
        else:
            raise LMError('Dataset features are empty.')

# ............................................... 
Example 7
Project: core   Author: lifemapper   File: createshape.py    GNU General Public License v3.0 5 votes vote down vote up
def _createFillFeat(self, lyrDef, recDict, lyr):
        feat = ogr.Feature(lyrDef)
        try:
            self._fillFeature(feat, recDict)
        except Exception, e:
            print('Failed to _createFillFeat, e = {}'.format(fromUnicode(toUnicode(e))))
            raise e 
Example 8
Project: Data-Tools-for-ArcGIS   Author: datastory   File: Service Geometry Downloader.py    MIT License 5 votes vote down vote up
def getParameterInfo(self):
		"""Define parameter definitions"""
		param0 = arcpy.Parameter(
			displayName="ArcGis Server Endpoint URL",
			name="url",
			datatype="GPString",
			parameterType="Required",
			direction="Input")

		param1 = arcpy.Parameter(
			displayName="Scratch Folder",
			name="scratch",
			datatype="DEFolder",
			parameterType="Required",
			direction="Input")

		param2 = arcpy.Parameter(
			displayName="Output Geodatabase",
			name="outDB",
			datatype="DEWorkspace",
			parameterType="Required",
			direction="Input")

		param3 = arcpy.Parameter(
			displayName="Output Feature Class",
			name="outFe",
			datatype="GPString",
			parameterType="Required",
			direction="Output")

		params = [param0, param1, param2, param3]
		return params 
Example 9
Project: qgis-geogiglight-plugin   Author: planetfederal   File: testgpkg.py    GNU General Public License v2.0 5 votes vote down vote up
def testAddingFeatureUsingOgr(self):
        dataSource =self._getOgrTestLayer()
        layer = dataSource.GetLayer()
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(123, 456)
        featureDefn = layer.GetLayerDefn()
        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(point)
        feature.SetField("fid", 5)
        feature.SetField("n", 5)
        layer.CreateFeature(feature)
        dataSource.Destroy() 
Example 10
Project: buzzard   Author: airware   File: _a_gdal_vector.py    Apache License 2.0 5 votes vote down vote up
def iter_features_driver(slicing, mask_poly, mask_rect, lyr):
        with contextlib.ExitStack() as stack:
            stack.push(lambda *args, **kwargs: lyr.ResetReading())
            if mask_poly is not None:
                lyr.SetSpatialFilter(mask_poly)
                stack.push(lambda *args, **kwargs: lyr.SetSpatialFilter(None))
            elif mask_rect is not None:
                lyr.SetSpatialFilterRect(*mask_rect)
                stack.push(lambda *args, **kwargs: lyr.SetSpatialFilter(None))

            start, stop, step = slicing.indices(len(lyr))
            indices = range(start, stop, step)
            ftr = None # Necessary to prevent the old swig bug
            if step == 1:
                lyr.SetNextByIndex(start)
                for i in indices:
                    ftr = lyr.GetNextFeature()
                    if ftr is None: # pragma: no cover
                        raise IndexError('Feature #{} not found'.format(i))
                    yield ftr
            else:
                for i in indices:
                    lyr.SetNextByIndex(i)
                    ftr = lyr.GetNextFeature()
                    if ftr is None: # pragma: no cover
                        raise IndexError('Feature #{} not found'.format(i))
                    yield ftr

        # Necessary to prevent the old swig bug
        # https://trac.osgeo.org/gdal/ticket/6749
        del slicing, mask_poly, mask_rect, ftr

    # insert_data implementation **************************************************************** ** 
Example 11
Project: buzzard   Author: airware   File: _a_gdal_vector.py    Apache License 2.0 5 votes vote down vote up
def insert_data(self, geom, geom_type, fields, index):
        geom = self._ogr_of_geom(geom, geom_type)
        with self.acquire_driver_object() as (_, lyr):
            ftr = ogr.Feature(lyr.GetLayerDefn())

            if geom is not None:
                err = ftr.SetGeometry(geom)
                if err: # pragma: no cover
                    raise ValueError('Could not set geometry (%s)' % str(gdal.GetLastErrorMsg()).strip('\n'))

                if not self.back_ds.allow_none_geometry and ftr.GetGeometryRef() is None: # pragma: no cover
                    raise ValueError(
                        'Invalid geometry inserted '
                        '(Set `allow_none_geometry=True` in Dataset constructor to silence)'
                    )

            if index >= 0:
                err = ftr.SetFID(index)
                if err: # pragma: no cover
                    raise ValueError('Could not set field id (%s)' % str(gdal.GetLastErrorMsg()).strip('\n'))
            for i, field in enumerate(fields):
                if field is not None:
                    err = ftr.SetField2(i, self._type_of_field_index[i](field))
                    if err: # pragma: no cover
                        raise ValueError('Could not set field #{} ({}) ({})'.format(
                            i, field, str(gdal.GetLastErrorMsg()).strip('\n')
                        ))
            passed = ftr.Validate(ogr.F_VAL_ALL, True)
            if not passed: # pragma: no cover
                raise ValueError('Invalid feature {} ({})'.format(
                    err, str(gdal.GetLastErrorMsg()).strip('\n')
                ))

            err = lyr.CreateFeature(ftr)
            if err: # pragma: no cover
                raise ValueError('Could not create feature {} ({})'.format(
                    err, str(gdal.GetLastErrorMsg()).strip('\n')
                )) 
Example 12
Project: pyeo   Author: clcr   File: coordinate_manipulation.py    GNU General Public License v3.0 5 votes vote down vote up
def write_geometry(geometry, out_path, srs_id=4326):
    """
    Saves the geometry in an ogr.Geometry object to a shapefile.

    Parameters
    ----------
    geometry
        An ogr.Geometry object
    out_path
        The location to save the output shapefile
    srs_id
        The projection of the output shapefile. Can be an EPSG number or a WKT string.

    Notes
    -----
    The shapefile consists of one layer named 'geometry'.


    """
    # TODO: Fix this needing an extra filepath on the end
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    if type(srs_id) is int:
        srs.ImportFromEPSG(srs_id)
    if type(srs_id) is str:
        srs.ImportFromWkt(srs_id)
    layer = data_source.CreateLayer(
        "geometry",
        srs,
        geom_type=geometry.GetGeometryType())
    feature_def = layer.GetLayerDefn()
    feature = ogr.Feature(feature_def)
    feature.SetGeometry(geometry)
    layer.CreateFeature(feature)
    data_source.FlushCache()
    data_source = None 
Example 13
Project: HackZurich   Author: RefugeeMatchmaking   File: test_shp.py    MIT License 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 14
Project: bitcoin-arbitrage-trading-bot   Author: mammuth   File: jqvmap.py    MIT License 5 votes vote down vote up
def output_ogr(self, output):
    driver = ogr.GetDriverByName( 'ESRI Shapefile' )
    if os.path.exists( output['file_name'] ):
      driver.DeleteDataSource( output['file_name'] )
    source = driver.CreateDataSource( output['file_name'] )
    layer = source.CreateLayer( self.layer_dfn.GetName(),
                                geom_type = self.layer_dfn.GetGeomType(),
                                srs = self.layer.GetSpatialRef() )

    for field in self.fields:
      fd = ogr.FieldDefn( str(field['name']), field['type'] )
      fd.SetWidth( field['width'] )
      if 'precision' in field:
        fd.SetPrecision( field['precision'] )
      layer.CreateField( fd )

    for geometry in self.geometries:
      if geometry.geom is not None:
        feature = ogr.Feature( feature_def = layer.GetLayerDefn() )
        for index, field in enumerate(self.fields):
          if field['name'] in geometry.properties:
            feature.SetField( index, geometry.properties[field['name']].encode('utf-8') )
          else:
            feature.SetField( index, '' )
        feature.SetGeometryDirectly(
          ogr.CreateGeometryFromWkb(
            shapely.wkb.dumps(
              geometry.geom
            )
          )
        )
        layer.CreateFeature( feature )
        feature.Destroy()

    source.Destroy() 
Example 15
Project: qgis-safecast-plugin   Author: OpenGeoLabs   File: __init__.py    GNU General Public License v2.0 5 votes vote down vote up
def _writeToOgrDataSource(self):
        fileExt = self.storageFormat.lower()
        filePath = os.path.splitext(self._fileName)[0] + '.{}'.format(fileExt)
        writer, msg = QgsVectorFileWriter.writeAsVectorFormat(
            self,
            filePath,
            self._provider.encoding(),
            self._provider.crs(),
            driverName=self.storageFormat
        )
        if writer != QgsVectorFileWriter.NoError:
            raise LoadError(
                self.tr("Unable to create SQLite datasource: {}").format(msg)
            )

        # set datasource to new OGR datasource
        self.setDataSource(filePath, self._layerName, 'ogr')
        self._provider = self.dataProvider()

        # create metadata layer
        if self.metadata and 'table' in self.metadata:
            ds = ogr.Open(filePath, True)
            layer_name = self.metadata['table']
            layer = ds.GetLayerByName(layer_name)
            if layer is None:
                layer = ds.CreateLayer(layer_name, None, ogr.wkbNone)
            layer_defn = layer.GetLayerDefn()
            if 'columns' in self.metadata:
                for key in self.metadata['columns']:
                    field = ogr.FieldDefn(key, ogr.OFTString)
                    layer.CreateField(field)

                feat = ogr.Feature(layer_defn)
                for key, value in list(self.metadata['columns'].items()):
                    feat.SetField(key, value)
                layer.CreateFeature(feat)
                feat = None 
Example 16
Project: graph-partition   Author: valiantljk   File: test_shp.py    GNU General Public License v2.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c']  # edgenames
        self.paths = (  [(1.0, 1.0), (2.0, 2.0)],
                        [(2.0, 2.0), (3.0, 3.0)],
                        [(0.9, 0.9), (4.0, 2.0)]
                    )
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 17
Project: cadasta-platform   Author: Cadasta   File: shape.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def write_shp_layer(self, loc_data):
        if loc_data['geometry.wkt'] == '':
            return
        geom = ogr.CreateGeometryFromWkt(loc_data['geometry.wkt'])
        layer_type = geom.GetGeometryName().lower()
        layer = self.shp_layers.get(layer_type, None)
        if layer is None:
            layer = self.create_shp_layer(layer_type)
            self.shp_layers[layer_type] = layer
        if layer:
            feature = ogr.Feature(layer.GetLayerDefn())
            feature.SetGeometry(geom)
            feature.SetField('id', loc_data['id'])
            layer.CreateFeature(feature)
            feature.Destroy() 
Example 18
Project: morpheo-plugin   Author: MorphoCity   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def create_shapedir(self):
        drv = ogr.GetDriverByName("ESRI Shapefile")
        shp = drv.CreateDataSource(self.path)
        lyr = shp.CreateLayer("nodes", None, ogr.wkbPoint)
        feature = ogr.Feature(lyr.GetLayerDefn())
        feature.SetGeometry(None)
        lyr.CreateFeature(feature)
        feature.Destroy() 
Example 19
Project: dask-geomodeling   Author: nens   File: test_geometry.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def create_geojson(abspath, polygons=10, bbox=None, ndim=2, projection="EPSG:4326"):
    """Create random triangle polygons inside bbox"""
    driver = ogr.GetDriverByName(str("GeoJSON"))
    driver.DeleteDataSource(abspath)
    datasource = driver.CreateDataSource(abspath)
    layer = datasource.CreateLayer(
        str("results"), get_sr(projection), geom_type=ogr.wkbPolygon
    )
    field_definition = ogr.FieldDefn(str("name"), ogr.OFTString)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn(str("id"), ogr.OFTInteger)
    layer.CreateField(field_definition)
    layer_definition = layer.GetLayerDefn()

    if np.isscalar(polygons):
        polygons = np.random.random((polygons, 3, ndim))
        bbox_min = np.asarray(bbox[:ndim])
        bbox_max = np.asarray(bbox[-ndim:])
        polygons = polygons * (bbox_max - bbox_min) + bbox_min

    for feature_id, coords in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for coord in coords:
            ring.AddPoint_2D(*coord)
        ring.AddPoint_2D(*coords[0])  # close the ring
        polygon = ogr.Geometry(ogr.wkbPolygon)
        polygon.AddGeometry(ring)
        feature = ogr.Feature(layer_definition)
        feature.SetGeometry(polygon)
        feature.SetField(str("name"), str("test"))
        feature.SetField(str("id"), feature_id + 10)
        layer.CreateFeature(feature)
    layer.SyncToDisk()
    datasource.SyncToDisk()

    return polygons 
Example 20
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 5 votes vote down vote up
def geom2mask(geom, r_ds):
    geom_srs = geom.GetSpatialReference()
    geom2 = geom_dup(geom) 
    #Create memory vector dataset and add geometry as new feature
    ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('tmp_ds')
    m_lyr = ogr_ds.CreateLayer('tmp_lyr', srs=geom_srs)
    feat = ogr.Feature(m_lyr.GetLayerDefn())
    feat.SetGeometryDirectly(geom2)
    m_lyr.CreateFeature(feat)
    mask = lyr2mask(m_lyr, r_ds)
    m_lyr = None
    ogr_ds = None
    return mask

#Old function, does not work with inner rings or complex geometries 
Example 21
Project: uv-align-distribute   Author: c30ra   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 22
Project: honours_project   Author: JFriel   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 23
Project: honours_project   Author: JFriel   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 24
Project: MEOPAR_AIS   Author: casey-h   File: 1_generate_tracks_from_AIS_DB_vectorized.py    MIT License 5 votes vote down vote up
def write_short_track_record(track_layer_val, first_pts_AIS, second_pts_AIS, track_obj, elapsed_seconds_val, track_len_m_val, speed_knots_val):
    
    # Generate a feature and and populate its fields.
    track_feature = ogr.Feature(track_layer_val.GetLayerDefn())
    track_feature.SetField("TrackID" , first_pts_AIS['segidx'])
    track_feature.SetField("mmsi" , first_pts_AIS['mmsi'])
    track_feature.SetField("elp_sec" , elapsed_seconds_val)
    # The date conversion below is madness, but apparrently the way to go: https://stackoverflow.com/questions/13703720/converting-between-datetime-timestamp-and-datetime64 
    track_feature.SetField("st_date" , pd.to_datetime(str(first_pts_AIS['time'])).strftime("%Y-%m-%d %H:%M:%S"))
    track_feature.SetField("en_date" , pd.to_datetime(str(second_pts_AIS['time'])).strftime("%Y-%m-%d %H:%M:%S"))
    track_feature.SetField("seg_len_km" , track_len_m_val / float(1000))
    track_feature.SetField("distspdkts" , speed_knots_val)
    track_feature.SetField("sogpt1",first_pts_AIS['sog'])
    track_feature.SetField("sogpt2",second_pts_AIS['sog'])
    track_feature.SetField("avg_sog",(first_pts_AIS['sog']+second_pts_AIS['sog'])/2)
    
    # Assign the geometry to the feature.
    track_feature.SetGeometry(track_obj)
    
    # Create the feature within the output layer, then reclaim assigned memory.
    track_layer_val.CreateFeature(track_feature)
    track_feature.Destroy()

# End Function write_track_record

# Function generate_short_GIS - Write out a GIS representation of the output
# of a short / point-to-point segmentation. 
Example 25
Project: DiffQue   Author: DiffQueDifficultyEsti   File: test_shp.py    MIT License 5 votes vote down vote up
def create_shapedir(self):
        drv = ogr.GetDriverByName("ESRI Shapefile")
        shp = drv.CreateDataSource(self.path)
        lyr = shp.CreateLayer("nodes", None, ogr.wkbPoint)
        feature = ogr.Feature(lyr.GetLayerDefn())
        feature.SetGeometry(None)
        lyr.CreateFeature(feature)
        feature.Destroy() 
Example 26
Project: aop-helpFinder   Author: jecarvaill   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def create_shapedir(self):
        drv = ogr.GetDriverByName("ESRI Shapefile")
        shp = drv.CreateDataSource(self.path)
        lyr = shp.CreateLayer("nodes", None, ogr.wkbPoint)
        feature = ogr.Feature(lyr.GetLayerDefn())
        feature.SetGeometry(None)
        lyr.CreateFeature(feature)
        feature.Destroy() 
Example 27
Project: voteswap   Author: sbuss   File: test_shp.py    MIT License 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 28
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 5 votes vote down vote up
def setDamFieldValues(self, feat, damType):
        """
        Set attribute for dam type.

        :param feat: OGR Feature of dam.
        :param damType: String describing dam type ('primary' or 'secondary').

        :return: Updated OGR Feature of dam.
        """
        feat.SetField("damType", damType)
        return feat 
Example 29
Project: Data-Tools-for-ArcGIS   Author: datastory   File: Service Geometry Downloader.py    MIT License 4 votes vote down vote up
def transform(infile, output, insrs, format_name):
    """Transform input file to output target file, based on input coordinate
    reference system to WGS84

    :param infile: name of the input file
    :param output: name of the output file
    :param insrs: epsg code of input file
    :param format_name: ogr format name
    """

    logging.info('Transforming %s from %s to %s' % (infile, insrs, output)) 
    in_srs = osr.SpatialReference()
    in_srs.ImportFromEPSG(insrs)
    out_srs = osr.SpatialReference()
    out_srs.ImportFromEPSG(4324)
    coordTrans = osr.CoordinateTransformation(in_srs, out_srs)

    in_dsn = ogr.Open(infile)
    in_layer = in_dsn.GetLayer()
    in_feature_definition = in_layer.GetLayerDefn()

    out_driver = ogr.GetDriverByName(format_name)
    out_dsn = out_driver.CreateDataSource(output)
    out_layer = out_dsn.CreateLayer(in_layer.GetName(),
            geom_type=in_layer.GetGeomType())

    # add fields
    for i in range(0, in_feature_definition.GetFieldCount()):
        fieldDefn = in_feature_definition.GetFieldDefn(i)
        out_layer.CreateField(fieldDefn)

    # get the output layer's feature definition
    out_feature_definition = out_layer.GetLayerDefn()

    # loop through the input features
    inFeature = in_layer.GetNextFeature()
    while inFeature:
        # get the input geometry
        geom = inFeature.GetGeometryRef().Clone()
        # reproject the geometry
        geom.Transform(coordTrans)
        # create a new feature
        outFeature = ogr.Feature(out_feature_definition)
        # set the geometry and attribute
        outFeature.SetGeometry(geom)
        for i in range(0, out_feature_definition.GetFieldCount()):
            outFeature.SetField(out_feature_definition.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
        # add the feature to the shapefile
        out_layer.CreateFeature(outFeature)
        # destroy the features and get the next input feature
        outFeature.Destroy()
        inFeature.Destroy()
        inFeature = in_layer.GetNextFeature()

    # close the shapefiles
    in_dsn.Destroy()
    out_dsn.Destroy() 
Example 30
Project: swat-tools   Author: DigitalAgricultureDiscovery   File: process.py    GNU General Public License v3.0 4 votes vote down vote up
def create_new_field_shapefile(self, field_shapefile, output_data):

        # set up the shapefile driver
        driver = ogr.GetDriverByName('ESRI Shapefile')

        # create the data source
        data_source = driver.CreateDataSource(self.results_dir + '/Output')

        # create the layer
        layer = data_source.CreateLayer('Field_Response', None, ogr.wkbPolygon)

        # add fields
        layer.CreateField(ogr.FieldDefn('Shape_Area', ogr.OFTReal))
        layer.CreateField(
            ogr.FieldDefn(str(self.fieldswat_output_type), ogr.OFTReal))

        for i in range(0, len(field_shapefile.shapes())):

            # get the field's shape (coordinates)
            field_shape = field_shapefile.shapes()[i]

            # get the field's shape area (float)
            field_record = field_shapefile.records()[i][0]

            # create feature
            feature = ogr.Feature(layer.GetLayerDefn())

            # set the attributes
            feature.SetField('Shape_Area', field_record)
            feature.SetField(str(self.fieldswat_output_type), output_data[i])

            # create polygon ring
            ring = ogr.Geometry(ogr.wkbLinearRing)
            for point in field_shape.points:
                ring.AddPoint(point[0], point[1])

            # create polygon
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # set the feature geometry using the wkbpolygon
            feature.SetGeometry(poly)

            # create the feature in the layer
            layer.CreateFeature(feature)

            # free up resources
            feature.Destroy()

        # free up resources
        data_source.Destroy() 
Example 31
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 32
Project: TAGGS   Author: jensdebruijn   File: shapefiles.py    MIT License 4 votes vote down vote up
def merge_shapefile_by_id(infile, outfile, layer_name, ID):
    dirname = os.path.dirname(outfile)
    try:
        os.makedirs(dirname)
    except OSError:
        pass
    inshp = ogr.Open(infile)
    inLayer = inshp.GetLayer()

    driver = ogr.GetDriverByName("ESRI Shapefile")

    outshp = driver.CreateDataSource(outfile)
    outLayer = outshp.CreateLayer(layer_name, geom_type=ogr.wkbPolygon)

    idField = ogr.FieldDefn(ID, ogr.OFTInteger)
    outLayer.CreateField(idField)

    gn_ids = [abs(feature.GetField(ID)) for feature in inLayer]
    duplicate_ids = function.find_duplicates(gn_ids)

    duplicate_ids_features = dd(list)

    outLayerDefn = outLayer.GetLayerDefn()

    for i in range(inLayer.GetFeatureCount()):
        inFeature = inLayer.GetFeature(i)
        outFeature = ogr.Feature(outLayerDefn)

        gn_id = abs(inFeature.GetField(ID))
        if gn_id not in duplicate_ids:

            outFeature.SetField(ID, gn_id)

            geom = inFeature.GetGeometryRef()
            outFeature.SetGeometry(geom)

            outLayer.CreateFeature(outFeature)
            outFeature = None

        else:
            duplicate_ids_features[gn_id].append(i)

    for gn_id, features in duplicate_ids_features.items():
        geom = ogr.Geometry(ogr.wkbPolygon)
        for i in features:
            inFeature = inLayer.GetFeature(i)
            geom = geom.Union(inFeature.GetGeometryRef())

        outFeature = ogr.Feature(outLayerDefn)
        outFeature.SetField(ID, gn_id)
        outFeature.SetGeometry(geom)
        outLayer.CreateFeature(outFeature)
        outFeature = None

    inshp = None
    outshp = None 
Example 33
Project: Planet-GEE-Pipeline-GUI   Author: samapriya   File: ogr2ft.py    Apache License 2.0 4 votes vote down vote up
def copy_features(src_layer, dst_layer, fix_geometry, simplify_geometry, start_index, total):
    index = 0
    batch_size = 200
    index_batch = 0
    for feat in src_layer:
        if index < start_index:
            index = index + 1
            continue

        try:
            geom = shapely.wkt.loads(feat.GetGeometryRef().ExportToWkt())
        except Exception as e:
            print('Error({0}), skipping geometry.'.format(e))
            continue

        if fix_geometry and not geom.is_valid:
            geom = geom.buffer(0.0)

        if simplify_geometry:
            geom = geom.simplify(0.004)

        f = ogr.Feature(dst_layer.GetLayerDefn())

        # set field values
        for i in range(feat.GetFieldCount()):
            fd = feat.GetFieldDefnRef(i)
            f.SetField(fd.GetName(), feat.GetField(fd.GetName()))

        # set geometry    
        f.SetGeometry(ogr.CreateGeometryFromWkt(geom.to_wkt()))

        if index_batch == 0:
            dst_layer.StartTransaction()

        # create feature
        feature = dst_layer.CreateFeature(f)

        f.Destroy()

        index_batch = index_batch + 1

        if index_batch >= batch_size or index == total - 1:
            dst_layer.CommitTransaction()
            count = dst_layer.GetFeatureCount() # update number of inserted features
            print('Inserted {0} of {1} features ({2:.2f}%)'.format(count, total, 100. * float(count) / total))

            index_batch = 0

            if index == total - 1:
                break

        index = index + 1 
Example 34
Project: Bandit   Author: paknorton   File: prms_geo.py    MIT License 4 votes vote down vote up
def write_shapefile3(self, filename, attr_name, attr_values, included_fields=None):
        # NOTE: This is for messing around with the geopackage format

        # Create a shapefile for the current selected layer
        # If a filter is set then a subset of features is written
        print(attr_values)
        out_driver = ogr.GetDriverByName('GPKG')

        out_ds = out_driver.CreateDataSource('crap.gpkg')
        # out_ds = out_driver.CreateDataSource(filename)
        out_layer = out_ds.CreateLayer(self.__selected_layer.GetName(), self.__selected_layer.GetSpatialRef())

        # Copy field definitions from input to output file
        in_layer_def = self.__selected_layer.GetLayerDefn()

        for ii in range(in_layer_def.GetFieldCount()):
            fld_def = in_layer_def.GetFieldDefn(ii)
            fld_name = fld_def.GetName()

            if included_fields and fld_name not in included_fields:
                continue
            out_layer.CreateField(fld_def)

        # Get feature definitions for the output layer
        out_layer_def = out_layer.GetLayerDefn()

        # Create blank output feature
        out_feat = ogr.Feature(out_layer_def)

        # Add features to the output layer
        for in_feat in self.__selected_layer:
            if in_feat.GetField(attr_name) in attr_values:
                print(in_feat.GetField(attr_name))
                # Add field values from the input layer
                for ii in range(out_layer_def.GetFieldCount()):
                    fld_def = out_layer_def.GetFieldDefn(ii)
                    fld_name = fld_def.GetName()

                    if included_fields and fld_name not in included_fields:
                        continue
                    out_feat.SetField(out_layer_def.GetFieldDefn(ii).GetNameRef(), in_feat.GetField(ii))

                # Set geometry as centroid
                # geom = in_feat.GetGeometryRef()
                # out_feat.SetGeometry(geom.Clone())

                # Set geometry
                geom = in_feat.geometry()
                out_feat.SetGeometry(geom)

                # Add the new feature to the output layer
                out_layer.CreateFeature(out_feat)

        # Close the output datasource
        out_ds.Destroy() 
Example 35
Project: RAPIDpy   Author: erdc   File: voronoi.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def pointsToVoronoiGridShapefile(lat, lon, vor_shp_path, extent=None):
    """
    Converts points to shapefile grid via voronoi
    """
    voronoi_centroids = _get_voronoi_centroid_array(lat, lon, extent)

    # set-up output polygon shp
    log("Creating output polygon shp {0}"
        .format(os.path.basename(vor_shp_path)))
    if os.path.exists(vor_shp_path):
        os.remove(vor_shp_path)
    drv = ogr.GetDriverByName('ESRI Shapefile')
    outShp = drv.CreateDataSource(vor_shp_path)
    osr_geographic_proj = osr.SpatialReference()
    osr_geographic_proj.ImportFromEPSG(4326)
    layer = outShp.CreateLayer('', osr_geographic_proj, ogr.wkbPolygon)
    layer.CreateField(ogr.FieldDefn('GRID_LAT', ogr.OFTReal))
    layer.CreateField(ogr.FieldDefn('GRID_LON', ogr.OFTReal))
    layerDefn = layer.GetLayerDefn()

    # find nodes surrounding polygon centroid
    # sort nodes in counterclockwise order
    # create polygon perimeter through nodes
    log("Building Voronoi polygons...")
    # compute voronoi
    voronoi_manager = Voronoi(voronoi_centroids)
    voronoi_vertices = voronoi_manager.vertices
    voronoi_regions = voronoi_manager.regions
    for point_id, region_index in enumerate(voronoi_manager.point_region):
        vert_index_list = np.array(voronoi_regions[region_index])
        voronoi_centroid = voronoi_centroids[point_id]
        voronoi_poly_points = _get_voronoi_poly_points(vert_index_list,
                                                       voronoi_vertices,
                                                       voronoi_centroid)
        if len(voronoi_poly_points) == 4:
            poly = ogr.Geometry(ogr.wkbPolygon)
            ring = ogr.Geometry(ogr.wkbLinearRing)
            for node in voronoi_poly_points:
                ring.AddPoint(node[0], node[1])

            # grab first node to close ring
            ring.AddPoint(voronoi_poly_points[0][0], voronoi_poly_points[0][1])

            poly.AddGeometry(ring)
            feat = ogr.Feature(layerDefn)
            feat.SetField('GRID_LON', float(voronoi_centroid[0]))
            feat.SetField('GRID_LAT', float(voronoi_centroid[1]))
            feat.SetGeometry(poly)
            layer.CreateFeature(feat) 
Example 36
Project: rastercube   Author: terrai   File: shputils.py    MIT License 4 votes vote down vote up
def polygon_to_shapefile(polygons, poly_sr, fname, fields_defs=None,
                         poly_attrs=None):
    """
    Write a set of polygons to a shapefile
    Args:
        polygons: a list of (lat, lng) tuples
        fields_defs: The list of fields for those polygons, as a tuple
                     (name, ogr type) for each field. For example :
                     [('field1', ogr.OFTReal), ('field2', ogr.OFTInteger)]
        poly_attrs: A list of dict containing the attributes for each polygon
                        [{'field1' : 1.0, 'field2': 42},
                         {'field1' : 3.0, 'field2': 60}]
    """
    shp_driver = ogr.GetDriverByName("ESRI Shapefile")
    out_ds = shp_driver.CreateDataSource(fname)
    assert out_ds is not None, "Failed to create temporary %s" % fname
    out_layer = out_ds.CreateLayer(fname, poly_sr, geom_type=ogr.wkbPolygon)

    has_attrs = fields_defs is not None
    if has_attrs:
        attrs_name = []
        for field_name, field_type in fields_defs:
            out_layer.CreateField(ogr.FieldDefn(field_name, field_type))
            attrs_name.append(field_name)

    layer_defn = out_layer.GetLayerDefn()
    for i, poly in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        # gdal uses the (x, y) convention => (lng, lat)
        for point in poly:
            ring.AddPoint(point[1], point[0])
        ring.AddPoint(poly[0][1], poly[0][0])  # re-add the start to close
        p = ogr.Geometry(ogr.wkbPolygon)
        p.AddGeometry(ring)

        out_feature = ogr.Feature(layer_defn)
        out_feature.SetGeometry(p)

        if has_attrs:
            attrs = poly_attrs[i]
            for field_name in attrs_name:
                out_feature.SetField(field_name, attrs[field_name])

        out_layer.CreateFeature(out_feature)

    out_feature.Destroy()
    out_ds.Destroy() 
Example 37
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 4 votes vote down vote up
def lyr_proj(lyr, t_srs, preserve_fields=True):
    """Reproject an OGR layer
    """
    #Need to check t_srs
    s_srs = lyr.GetSpatialRef()
    cT = osr.CoordinateTransformation(s_srs, t_srs)

    #Do everything in memory
    drv = ogr.GetDriverByName('Memory')

    #Might want to save clipped, warped shp to disk?
    # create the output layer
    #drv = ogr.GetDriverByName('ESRI Shapefile')
    #out_fn = '/tmp/temp.shp'
    #if os.path.exists(out_fn):
    #    driver.DeleteDataSource(out_fn)
    #out_ds = driver.CreateDataSource(out_fn)
    
    out_ds = drv.CreateDataSource('out')
    outlyr = out_ds.CreateLayer('out', srs=t_srs, geom_type=lyr.GetGeomType())

    if preserve_fields:
        # add fields
        inLayerDefn = lyr.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            outlyr.CreateField(fieldDefn)
        # get the output layer's feature definition
    outLayerDefn = outlyr.GetLayerDefn()

    # loop through the input features
    inFeature = lyr.GetNextFeature()
    while inFeature:
        # get the input geometry
        geom = inFeature.GetGeometryRef()
        # reproject the geometry
        geom.Transform(cT)
        # create a new feature
        outFeature = ogr.Feature(outLayerDefn)
        # set the geometry and attribute
        outFeature.SetGeometry(geom)
        if preserve_fields:
            for i in range(0, outLayerDefn.GetFieldCount()):
                outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
        # add the feature to the shapefile
        outlyr.CreateFeature(outFeature)
        # destroy the features and get the next input feature
        inFeature = lyr.GetNextFeature()
    #NOTE: have to operate on ds here rather than lyr, otherwise segfault
    return out_ds

#See https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#convert-vector-layer-to-array
#Should check srs, as shp could be WGS84 
Example 38
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 4 votes vote down vote up
def geom2shp(geom, out_fn, fields=False):
    """Write out a new shapefile for input geometry
    """
    from pygeotools.lib import timelib
    driverName = "ESRI Shapefile"
    drv = ogr.GetDriverByName(driverName)
    if os.path.exists(out_fn):
        drv.DeleteDataSource(out_fn)
    out_ds = drv.CreateDataSource(out_fn)
    out_lyrname = os.path.splitext(os.path.split(out_fn)[1])[0]
    geom_srs = geom.GetSpatialReference()
    geom_type = geom.GetGeometryType()
    out_lyr = out_ds.CreateLayer(out_lyrname, geom_srs, geom_type)
    if fields:
        field_defn = ogr.FieldDefn("name", ogr.OFTString)
        field_defn.SetWidth(128)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("path", ogr.OFTString)
        field_defn.SetWidth(254)
        out_lyr.CreateField(field_defn)
        #field_defn = ogr.FieldDefn("date", ogr.OFTString)
        #This allows sorting by date
        field_defn = ogr.FieldDefn("date", ogr.OFTInteger)
        field_defn.SetWidth(32)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("decyear", ogr.OFTReal)
        field_defn.SetPrecision(8)
        field_defn.SetWidth(64)
        out_lyr.CreateField(field_defn)
    out_feat = ogr.Feature(out_lyr.GetLayerDefn())
    out_feat.SetGeometry(geom)
    if fields:
        #Hack to force output extesion to tif, since out_fn is shp
        out_path = os.path.splitext(out_fn)[0] + '.tif'
        out_feat.SetField("name", os.path.split(out_path)[-1])
        out_feat.SetField("path", out_path)
        #Try to extract a date from input raster fn
        out_feat_date = timelib.fn_getdatetime(out_fn)
        if out_feat_date is not None:
            datestamp = int(out_feat_date.strftime('%Y%m%d'))
            #out_feat_date = int(out_feat_date.strftime('%Y%m%d%H%M'))
            out_feat.SetField("date", datestamp)
            decyear = timelib.dt2decyear(out_feat_date)
            out_feat.SetField("decyear", decyear)
    out_lyr.CreateFeature(out_feat)
    out_ds = None
    #return status? 
Example 39
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 4 votes vote down vote up
def createDams(self):
        """
        Place primary and secondary dams at a specific location on stream reaches.

        :return: None
        """
        i = 0
        while i < self.nFeat:
            bratFt = self.bratLyr.GetFeature(self.capRank[i, 0])
            bratLine = bratFt.GetGeometryRef()
            length = bratLine.Length()
            #read number of dams and dam complexes from shapefile field
            nDamCt = bratFt.GetFieldAsInteger("totdams")
            nCompCt = bratFt.GetFieldAsInteger("totcomp")
            spacing = 0.0
            if nDamCt > 0:
                spacing = length / (nDamCt * 1.0)

            #create a point for each dam to be modeled
            nDamRm = nDamCt
            nCompRm = nCompCt
            for j in range(0, nDamCt):
                #determine if dam is primary or secondary and create height distribution
                damFeat = ogr.Feature(self.outLyr.GetLayerDefn())
                rnum = np.random.random()
                if rnum < ((nCompRm*1.0)/(nDamRm*1.0)) or nDamCt == 1:
                    htDist = np.random.lognormal(0.22, 0.36, 30)
                    damType = "primary"
                    nCompRm -= 1
                else:
                    htDist = np.random.lognormal(-0.21, 0.39, 30)
                    damType = "secondary"

                nDamRm -= 1
                #location of dam on stream segment
                pointDist = length - (spacing * (j * 1.0))
                damPoint = bratLine.Value(pointDist)
                #bratLine.Value will return None if called on multipart line. This will cause a crash later in the code
                if damPoint is not None:
                    #set any field values here
                    damFeat = self.setDamFieldValues(damFeat, damType)
                    #set dam heights
                    damFeat = self.setDamHeights(damFeat, np.percentile(htDist, 2.5), np.median(htDist), np.percentile(htDist, 97.5))
                    damFeat.SetGeometry(damPoint)
                    self.outLyr.CreateFeature(damFeat)
                damFeat = None

            i += 1
            bratFt = None