Python osgeo.ogr.GetDriverByName() Examples

The following are code examples for showing how to use osgeo.ogr.GetDriverByName(). They are extracted from open source Python projects. You can vote up the examples you like or vote down the ones you don't like. You can also save this page to your account.

Example 1
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 8 votes vote down vote up
def exporttogeojson(geojsonfilename, buildinglist):
    #
    # geojsonname should end with .geojson
    # building list should be list of dictionaries
    # list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'polyPix': poly,
    #                       'polyGeo': poly}
    # image_id is a string,
    # BuildingId is an integer,
    # poly is a ogr.Geometry Polygon
    #
    # returns geojsonfilename

    driver = ogr.GetDriverByName('geojson')
    if os.path.exists(geojsonfilename):
        driver.DeleteDataSource(geojsonfilename)
    datasource = driver.CreateDataSource(geojsonfilename)
    layer = datasource.CreateLayer('buildings', geom_type=ogr.wkbPolygon)
    field_name = ogr.FieldDefn("ImageId", ogr.OFTString)
    field_name.SetWidth(75)
    layer.CreateField(field_name)
    layer.CreateField(ogr.FieldDefn("BuildingId", ogr.OFTInteger))

    # loop through buildings
    for building in buildinglist:
        # create feature
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField("ImageId", building['ImageId'])
        feature.SetField("BuildingId", building['BuildingId'])
        feature.SetGeometry(building['polyPix'])

        # Create the feature in the layer (geojson)
        layer.CreateFeature(feature)
        # Destroy the feature to free resources
        feature.Destroy()

    datasource.Destroy()

    return geojsonfilename 
Example 2
Project: Planet-Pipeline-GUI   Author: samapriya   File: kml_aoi.py    (license) View Source Project 6 votes vote down vote up
def parsed(args):
        kml_file=args.geo
        def kml2geojson(kml_file):
            drv = ogr.GetDriverByName('KML')
            kml_ds = drv.Open(kml_file)
            for kml_lyr in kml_ds:
                for feat in kml_lyr:
                    outfile=feat.ExportToJson()
                    geom2=str(outfile).replace(", 0.0",'')
                    with open(args.loc+'./kmlout.geojson','w') as csvfile:
                        writer=csv.writer(csvfile)
                        writer.writerow([geom2])
        kml2geojson(args.geo)
        raw= open(args.loc+'./kmlout.geojson')
        for line in raw:
            fields=line.strip().split(":")[3]
            f2=fields.strip().split("}")[0]
            filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
        with open(args.loc+'./aoi.json', 'w') as outfile:
            outfile.write(filenames)
            outfile.close() 
Example 3
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 6 votes vote down vote up
def polygonize(self,rasterTemp,outShp):
            
            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)
                
            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile
      
            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)
            
            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  
            
            outDatasource.Destroy()
            sourceRaster=None
            band=None
                    
            
            ioShpFile = ogr.Open(outShp, update = 1)
            
            
            lyr = ioShpFile.GetLayerByIndex(0)
            
            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()
                        
            return outShp 
Example 4
Project: Planet-GEE-Pipeline-GUI   Author: samapriya   File: kml_aoi.py    (license) View Source Project 6 votes vote down vote up
def parsed(args):
        kml_file=args.geo
        def kml2geojson(kml_file):
            drv = ogr.GetDriverByName('KML')
            kml_ds = drv.Open(kml_file)
            for kml_lyr in kml_ds:
                for feat in kml_lyr:
                    outfile=feat.ExportToJson()
                    geom2=str(outfile).replace(", 0.0",'')
                    with open(args.loc+'./kmlout.geojson','w') as csvfile:
                        writer=csv.writer(csvfile)
                        writer.writerow([geom2])
        kml2geojson(args.geo)
        raw= open(args.loc+'./kmlout.geojson')
        for line in raw:
            fields=line.strip().split(":")[3]
            f2=fields.strip().split("}")[0]
            filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
        with open(args.loc+'./aoi.json', 'w') as outfile:
            outfile.write(filenames)
            outfile.close() 
Example 5
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def import_summary_geojson(geojsonfilename, removeNoBuildings=True):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)

    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    buildingList = []
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            if removeNoBuildings:
                if feature.GetField('BuildingId') != -1:
                    buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                                  'poly': feature.GetGeometryRef().Clone()})
            else:

                buildingList.append({'ImageId': feature.GetField('ImageId'), 'BuildingId': feature.GetField('BuildingId'),
                              'poly': feature.GetGeometryRef().Clone()})

    return buildingList 
Example 6
Project: utilities   Author: SpaceNetChallenge   File: geoTools.py    (license) View Source Project 6 votes vote down vote up
def import_chip_geojson(geojsonfilename, ImageId=''):
    # driver = ogr.GetDriverByName('geojson')
    datasource = ogr.Open(geojsonfilename, 0)
    if ImageId=='':
        ImageId = geojsonfilename
    layer = datasource.GetLayer()
    print(layer.GetFeatureCount())

    polys = []
    BuildingId = 0
    for idx, feature in enumerate(layer):

        poly = feature.GetGeometryRef()

        if poly:
            BuildingId = BuildingId + 1
            polys.append({'ImageId': ImageId, 'BuildingId': BuildingId,
                          'poly': feature.GetGeometryRef().Clone()})

    return polys 
Example 7
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 6 votes vote down vote up
def createRasterFromGeoJson(srcGeoJson, srcRasterFileName, outRasterFileName):
    NoData_value = 0
    source_ds = ogr.Open(srcGeoJson)
    source_layer = source_ds.GetLayer()

    srcRaster = gdal.Open(srcRasterFileName)


    # Create the destination data source
    target_ds = gdal.GetDriverByName('GTiff').Create(outRasterFileName, srcRaster.RasterXSize, srcRaster.RasterYSize, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform(srcRaster.GetGeoTransform())
    target_ds.SetProjection(srcRaster.GetProjection())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(NoData_value)

    # Rasterize
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
    band.FlushCache() 
Example 8
Project: Planet-GEE-Pipeline-CLI   Author: samapriya   File: kml_aoi.py    (license) View Source Project 6 votes vote down vote up
def parsed(args):
        kml_file=args.geo
        def kml2geojson(kml_file):
            drv = ogr.GetDriverByName('KML')
            kml_ds = drv.Open(kml_file)
            for kml_lyr in kml_ds:
                for feat in kml_lyr:
                    outfile=feat.ExportToJson()
                    geom2=str(outfile).replace(", 0.0",'')
                    with open(args.loc+'./kmlout.geojson','w') as csvfile:
                        writer=csv.writer(csvfile)
                        writer.writerow([geom2])
        kml2geojson(args.geo)
        raw= open(args.loc+'./kmlout.geojson')
        for line in raw:
            fields=line.strip().split(":")[3]
            f2=fields.strip().split("}")[0]
            filenames = p1+f2+p2+str(args.start)+p3+str(args.end)+p4+p5+str(args.cloud)+p6
        with open(args.loc+'./aoi.json', 'w') as outfile:
            outfile.write(filenames)
            outfile.close() 
Example 9
Project: rastercube   Author: terrai   File: shputils.py    (license) View Source Project 6 votes vote down vote up
def rasterize_shapefile_like(shpfile, model_raster_fname, dtype, options,
                             nodata_val=0,):
    """
    Given a shapefile, rasterizes it so it has
    the exact same extent as the given model_raster

    `dtype` is a gdal type like gdal.GDT_Byte
    `options` should be a list that will be passed to GDALRasterizeLayers
        papszOptions, like ["ATTRIBUTE=vegetation","ALL_TOUCHED=TRUE"]
    """
    model_dataset = gdal.Open(model_raster_fname)
    shape_dataset = ogr.Open(shpfile)
    shape_layer = shape_dataset.GetLayer()
    mem_drv = gdal.GetDriverByName('MEM')
    mem_raster = mem_drv.Create(
        '',
        model_dataset.RasterXSize,
        model_dataset.RasterYSize,
        1,
        dtype
    )
    mem_raster.SetProjection(model_dataset.GetProjection())
    mem_raster.SetGeoTransform(model_dataset.GetGeoTransform())
    mem_band = mem_raster.GetRasterBand(1)
    mem_band.Fill(nodata_val)
    mem_band.SetNoDataValue(nodata_val)

    # http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
    # http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
    err = gdal.RasterizeLayer(
        mem_raster,
        [1],
        shape_layer,
        None,
        None,
        [1],
        options
    )
    assert err == gdal.CE_None
    return mem_raster.ReadAsArray() 
Example 10
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 6 votes vote down vote up
def mem_ds(res, extent, srs=None, dtype=gdal.GDT_Float32):
    """Create a new GDAL Dataset in memory

    Useful for various applications that require a Dataset
    """
    #These round down to int
    #dst_ns = int((extent[2] - extent[0])/res)
    #dst_nl = int((extent[3] - extent[1])/res)
    #This should pad by 1 pixel, but not if extent and res were calculated together to give whole int
    dst_ns = int((extent[2] - extent[0])/res + 0.99)
    dst_nl = int((extent[3] - extent[1])/res + 0.99)
    m_ds = gdal.GetDriverByName('MEM').Create('', dst_ns, dst_nl, 1, dtype)
    m_gt = [extent[0], res, 0, extent[3], 0, -res]
    m_ds.SetGeoTransform(m_gt)
    if srs is not None:
        m_ds.SetProjection(srs.ExportToWkt())
    return m_ds

#Modify proj/gt of dst_fn in place 
Example 11
Project: pyroSAR   Author: johntruckenbrodt   File: vector.py    (license) View Source Project 6 votes vote down vote up
def __init__(self, filename=None, driver="ESRI Shapefile"):

        if driver not in ["ESRI Shapefile", "Memory"]:
            raise IOError("driver not supported")

        if filename is None:
            driver = "Memory"
        else:
            self.filename = filename

        self.driver = ogr.GetDriverByName(driver)

        self.vector = self.driver.CreateDataSource("out") if driver == "Memory" else self.driver.Open(filename)

        nlayers = self.vector.GetLayerCount()
        if nlayers > 1:
            raise IOError("multiple layers are currently not supported")
        elif nlayers == 1:
            self.init_layer() 
Example 12
Project: pyroSAR   Author: johntruckenbrodt   File: vector.py    (license) View Source Project 6 votes vote down vote up
def write(self, outfile, format="ESRI Shapefile", overwrite=True):
        (outfilepath, outfilename) = os.path.split(outfile)
        basename = os.path.splitext(outfilename)[0]

        driver = ogr.GetDriverByName(format)

        if os.path.exists(outfile):
            if overwrite:
                driver.DeleteDataSource(outfile)
            else:
                raise IOError("file already exists")

        outdataset = driver.CreateDataSource(outfile)
        outlayer = outdataset.CreateLayer(self.layername, geom_type=self.geomType)
        outlayerdef = outlayer.GetLayerDefn()

        for fieldDef in self.fieldDefs:
            outlayer.CreateField(fieldDef)

        for feature in self.layer:
            outFeature = ogr.Feature(outlayerdef)
            outFeature.SetGeometry(feature.GetGeometryRef())
            for j in range(0, self.nfields):
                outFeature.SetField(self.fieldnames[j], feature.GetField(j))
            # add the feature to the shapefile
            outlayer.CreateFeature(outFeature)
            outFeature.Destroy()

        srs_out = self.srs.Clone()
        srs_out.MorphToESRI()
        with open(os.path.join(outfilepath, basename+".prj"), "w") as prj:
            prj.write(srs_out.ExportToWkt())

        outdataset.Destroy() 
Example 13
Project: unmixing   Author: arthur-e   File: lsma.py    (license) View Source Project 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 14
Project: gee-bridge   Author: francbartoli   File: wpCalc.py    (license) View Source Project 5 votes vote down vote up
def generate_tiles(self):

        """INCOMPLETE  Split the calculated WPbm in 100 tiles facilitating the export

        Returns:
            TYPE: Description
        """

        driver = ogr.GetDriverByName('ESRI Shapefile')
        dir_shps = "tiles"
        os.chdir(dir_shps)
        file_shps = glob.glob("*.shp")

        allExportWPbm = []
        file_names = []

        for file_shp in file_shps:

            dataSource = driver.Open(file_shp, 0)

            if dataSource is None:
                sys.exit(('Could not open {0}.'.format(file_shp)))
            else:
                layer = dataSource.GetLayer(0)
                extent = layer.GetExtent()
                active_file = "tile_" + str(file_shp.split('.')[0]).split("_")[3]
                file_names.append(active_file)
                low_sx = extent[0], extent[3]
                up_sx = extent[0], extent[2]
                up_dx = extent[1], extent[2]
                low_dx = extent[1], extent[3]

                cut = [list(low_sx), list(up_sx), list(up_dx), list(low_dx)]

                Export_WPbm = {
                    "crs": "EPSG:4326",
                    "scale": 250,
                    'region': cut}
                allExportWPbm.append(Export_WPbm)

        return allExportWPbm, file_names 
Example 15
Project: sanergy-public   Author: dssg   File: kml_to_geojson.py    (license) View Source Project 5 votes vote down vote up
def kml2geojson(kml_file):
    '''
    From Gist. 
    Transform kml (from Google Maps) to the geojson format.
    '''
    drv = ogr.GetDriverByName('KML')
    kml_ds = drv.Open(kml_file)
    for kml_lyr in kml_ds:
        for feat in kml_lyr:
            print feat.ExportToJson() 
Example 16
Project: pfb-network-connectivity   Author: azavea   File: detect_utm_zone.py    (license) View Source Project 5 votes vote down vote up
def get_bbox(shapefile):
    """ Gets the bounding box of a shapefile in EPSG 4326.
        If shapefile is not in WGS84, bounds are reprojected.
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 0 means read, 1 means write
    data_source = driver.Open(shapefile, 0)
    # Check to see if shapefile is found.
    if data_source is None:
        failure('Could not open %s' % (shapefile))
    else:
        layer = data_source.GetLayer()
        shape_bbox = layer.GetExtent()

        spatialRef = layer.GetSpatialRef()
        target = osr.SpatialReference()
        target.ImportFromEPSG(4326)

        # this check for non-WGS84 projections gets some false positives, but that's okay
        if target.ExportToProj4() != spatialRef.ExportToProj4():
            transform = osr.CoordinateTransformation(spatialRef, target)
            point1 = ogr.Geometry(ogr.wkbPoint)
            point1.AddPoint(shape_bbox[0], shape_bbox[2])
            point2 = ogr.Geometry(ogr.wkbPoint)
            point2.AddPoint(shape_bbox[1], shape_bbox[3])
            point1.Transform(transform)
            point2.Transform(transform)
            shape_bbox = [point1.GetX(), point2.GetX(), point1.GetY(), point2.GetY()]

        return shape_bbox 
Example 17
Project: Planet-GEE-Pipeline-GUI   Author: samapriya   File: ogr2ft.py    (license) View Source Project 5 votes vote down vote up
def _get_ft_ds():
    refresh_token = OAuth2().get_refresh_token()
    ft_driver = ogr.GetDriverByName('GFT')

    return ft_driver.Open('GFT:refresh=' + refresh_token, True) 
Example 18
Project: utilities   Author: SpaceNetChallenge   File: externalVectorProcessing.py    (license) View Source Project 5 votes vote down vote up
def buildTindex(rasterFolder, rasterExtention='.tif'):
    rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
    print(rasterList)

    print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))

    memDriver = ogr.GetDriverByName('MEMORY')
    gTindex = memDriver.CreateDataSource('gTindex')
    srcImage = gdal.Open(rasterList[0])
    spat_ref = osr.SpatialReference()
    spat_ref.SetProjection(srcImage.GetProjection())
    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("location", ogr.OFTString)
    gTindexLayer.CreateField(idField)

    # Create the feature and set values
    featureDefn = gTindexLayer.GetLayerDefn()



    for rasterFile in rasterList:
        srcImage = gdal.Open(rasterFile)

        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)

        feature = ogr.Feature(featureDefn)
        feature.SetGeometry(polyToCut)
        feature.SetField("location", rasterFile)
        gTindexLayer.CreateFeature(feature)
        feature = None


    return gTindex, gTindexLayer 
Example 19
Project: wikiwhere   Author: mkrnr   File: countries.py    (license) View Source Project 5 votes vote down vote up
def __init__(self, country_file):
        driver = ogr.GetDriverByName('ESRI Shapefile')
        self.countryFile = driver.Open(country_file)
        self.layer = self.countryFile.GetLayer() 
Example 20
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 5 votes vote down vote up
def open_vector(filename, driver=None):
    """Open vector file, return gdal.Dataset and OGR.Layer

        .. warning:: dataset and layer have to live in the same context,
            if dataset is deleted all layer references will get lost

        .. versionadded:: 0.12.0

    Parameters
    ----------
    filename : string
        vector file name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    layer : ogr.Layer
        layer
    """
    dataset = gdal.OpenEx(filename)

    if driver:
        gdal.GetDriverByName(driver)

    layer = dataset.GetLayer()

    return dataset, layer 
Example 21
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 5 votes vote down vote up
def open_shape(filename, driver=None):
    """Open shapefile, return gdal.Dataset and OGR.Layer

        .. warning:: dataset and layer have to live in the same context,
            if dataset is deleted all layer references will get lost

        .. versionadded:: 0.6.0

    Parameters
    ----------
    filename : string
        shapefile name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    layer : ogr.Layer
        layer
    """

    if driver is None:
        driver = ogr.GetDriverByName('ESRI Shapefile')
    dataset = driver.Open(filename)
    if dataset is None:
        print('Could not open file')
        raise IOError
    layer = dataset.GetLayer()
    return dataset, layer 
Example 22
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 5 votes vote down vote up
def open_raster(filename, driver=None):
    """Open raster file, return gdal.Dataset

        .. versionadded:: 0.6.0

    Parameters
    ----------
    filename : string
        raster file name
    driver : string
        gdal driver string

    Returns
    -------
    dataset : gdal.Dataset
        dataset
    """

    dataset = gdal.OpenEx(filename)

    if driver:
        gdal.GetDriverByName(driver)

    return dataset 
Example 23
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 5 votes vote down vote up
def read_safnwc(filename):
    """Read MSG SAFNWC hdf5 file into a gdal georeferenced object

    Parameters
    ----------
    filename : string
        satellite file name

    Returns
    -------
    ds : gdal.DataSet
        with satellite data
    """

    root = gdal.Open(filename)
    ds1 = gdal.Open('HDF5:' + filename + '://CT')
    ds = gdal.GetDriverByName('MEM').CreateCopy('out', ds1, 0)

    # name = os.path.basename(filename)[7:11]
    try:
        proj = osr.SpatialReference()
        proj.ImportFromProj4(ds.GetMetadata()["PROJECTION"])
    except Exception:
        raise NameError("No metadata for satellite file %s" % filename)
    geotransform = root.GetMetadata()["GEOTRANSFORM_GDAL_TABLE"].split(",")
    geotransform[0] = root.GetMetadata()["XGEO_UP_LEFT"]
    geotransform[3] = root.GetMetadata()["YGEO_UP_LEFT"]
    ds.SetProjection(proj.ExportToWkt())
    ds.SetGeoTransform([float(x) for x in geotransform])
    return ds 
Example 24
Project: dzetsaka   Author: lennepkade   File: function_vector.py    (license) View Source Project 5 votes vote down vote up
def saveToShape(self,array,srs,outShapeFile):
        # Parse a delimited text file of volcano data and create a shapefile
        # use a dictionary reader so we can access by field name
        # set up the shapefile driver
        outDriver = ogr.GetDriverByName( 'ESRI Shapefile' )
        
        # create the data source
        if os.path.exists(outShapeFile):
            outDriver.DeleteDataSource(outShapeFile)
        # Remove output shapefile if it already exists
        
        ds = outDriver.CreateDataSource(outShapeFile) #options = ['SPATIALITE=YES'])
    
        # create the spatial reference, WGS84
        
        lyrout = ds.CreateLayer('randomSubset',srs)
        fields = [array[1].GetFieldDefnRef(i).GetName() for i in range(array[1].GetFieldCount())]
        
        for f in fields:
            field_name = ogr.FieldDefn(f, ogr.OFTString)
            field_name.SetWidth(24)
            lyrout.CreateField(field_name)
            
        
        for k in array:
            lyrout.CreateFeature(k)
    
        # Save and close the data source
        ds = None 
Example 25
Project: qgis-geogiglight-plugin   Author: boundlessgeo   File: testgpkg.py    (license) View Source Project 5 votes vote down vote up
def _getOgrTestLayer(self):
        dest = self._copyTestLayer()
        driver = ogr.GetDriverByName('GPKG')
        dataSource = driver.Open(dest, 1)
        return dataSource 
Example 26
Project: CHaMP_Metrics   Author: SouthForkResearch   File: tin.py    (license) View Source Project 5 votes vote down vote up
def export_shp_nodes(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbPoint25D)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, point in self.tin_nodes.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(point.wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None

        return 
Example 27
Project: CHaMP_Metrics   Author: SouthForkResearch   File: tin.py    (license) View Source Project 5 votes vote down vote up
def export_shp_triangles(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbPolygon)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, poly in self.triangles.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(poly.wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None 
Example 28
Project: CHaMP_Metrics   Author: SouthForkResearch   File: tin.py    (license) View Source Project 5 votes vote down vote up
def export_shp_hull(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbPolygon)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, poly in self.hull_polygons.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(poly.wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None 
Example 29
Project: CHaMP_Metrics   Author: SouthForkResearch   File: tin.py    (license) View Source Project 5 votes vote down vote up
def export_shp_breaklines(self, outshp):
        from osgeo import ogr

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        ds = driver.CreateDataSource(outshp)
        layer = ds.CreateLayer('', None, ogr.wkbLineString25D)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        fieldLineType = ogr.FieldDefn('LineType', ogr.OFTString)
        fieldLineType.SetWidth(4)
        layer.CreateField(fieldLineType)

        defn = layer.GetLayerDefn()

        ## If there are multiple geometries, put the "for" loop here
        for id, line in self.breaklines.iteritems():
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', int(id))
            feat.SetField('LineType', line['linetype'])

            # Make a geometry, from Shapely object
            geom = ogr.CreateGeometryFromWkb(line['geometry'].wkb)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None 
Example 30
Project: antares   Author: CONABIO   File: split.py    (license) View Source Project 5 votes vote down vote up
def get_convex_hull(shape_name, destination_directory):
    '''
    This method will read all objects from a shape file and create a new one
    with the convex hull of all the geometry points of the first.
    '''
    driver = ogr.GetDriverByName(str('ESRI Shapefile'))
    shape = driver.Open(shape_name, 0)
    layer = shape.GetLayer()
    layer_name = layer.GetName()
    spatial_reference = layer.GetSpatialRef()
    prefix = get_basename(shape_name)
    output_name = create_filename(destination_directory, '%s-hull.shp' % prefix)
    geometries = ogr.Geometry(ogr.wkbGeometryCollection)
    for feature in layer:
        geometries.AddGeometry(feature.GetGeometryRef())
    if os.path.exists(output_name):
        driver.DeleteDataSource(output_name)
    datasource = driver.CreateDataSource(output_name)
    out_layer = datasource.CreateLayer(str('states_convexhull'), spatial_reference, geom_type=ogr.wkbPolygon)
    out_layer.CreateField(ogr.FieldDefn(str('id'), ogr.OFTInteger))
    featureDefn = out_layer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(geometries.ConvexHull())
    feature.SetField(str('id'), 1)
    out_layer.CreateFeature(feature)
    shape.Destroy()
    datasource.Destroy() 
Example 31
Project: HistoricalMap   Author: lennepkade   File: function_historical_map.py    (license) View Source Project 4 votes vote down vote up
def postClassRaster(self,inRaster,sieveSize,inClassNumber,outShp):        
        """ [email protected] Sieve size with gdal.Sieve() fiunction, them reclass to delete unwanted labels """
        
        try:
            rasterTemp = tempfile.mktemp('.tif')
                    
            datasrc = gdal.Open(inRaster)
            srcband = datasrc.GetRasterBand(1)
            data,im=dataraster.open_data_band(inRaster)        
            
            drv = gdal.GetDriverByName('GTiff')
            dst_ds = drv.Create(rasterTemp,datasrc.RasterXSize,datasrc.RasterXSize,1,gdal.GDT_Byte)
            
            dst_ds.SetGeoTransform(datasrc.GetGeoTransform())
            dst_ds.SetProjection(datasrc.GetProjection())
        
            dstband=dst_ds.GetRasterBand(1)
            
            
            
            def sieve(srcband,dstband,sieveSize):
                gdal.SieveFilter(srcband,None,dstband,sieveSize,8)
            
            pixelSize = datasrc.GetGeoTransform()[1] #get pixel size
            pixelSieve = int(sieveSize/(pixelSize*pixelSize)) #get number of pixel to sieve
            
            sieve(srcband,dstband,pixelSieve)
            
            dst_ds = None # close destination band
            
            self.Progress.addStep()
            

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

        
        try :
            outShp = self.polygonize(rasterTemp,outShp) # vectorize raster
        except :
            self.Progress.reset()    
        
        self.Progress.addStep()        
        
        self.Progress.reset()
        return outShp 
Example 32
Project: utilities   Author: SpaceNetChallenge   File: labelTools.py    (license) View Source Project 4 votes vote down vote up
def convertLabelStringToPoly(shapeFileSrc, outGeoJSon, labelType='Airplane'):

        shapeSrc = ogr.Open(shapeFileSrc)
        source_layer = shapeSrc.GetLayer()
        source_srs = source_layer.GetSpatialRef()
        # Create the output Layer
        outDriver = ogr.GetDriverByName("geojson")
        if os.path.exists(outGeoJSon):
            outDriver.DeleteDataSource(outGeoJSon)


        outDataSource = outDriver.CreateDataSource(outGeoJSon)
        outLayer = outDataSource.CreateLayer("groundTruth", source_srs, geom_type=ogr.wkbPolygon)
        # Add input Layer Fields to the output Layer
        inLayerDefn = source_layer.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            outLayer.CreateField(fieldDefn)
        outLayer.CreateField(ogr.FieldDefn("Length_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Width_m", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("Aspect(L/W)", ogr.OFTReal))
        outLayer.CreateField(ogr.FieldDefn("compassDeg", ogr.OFTReal))

        outLayerDefn = outLayer.GetLayerDefn()
        for inFeature in source_layer:

            outFeature = ogr.Feature(outLayerDefn)

            for i in range(0, inLayerDefn.GetFieldCount()):
                outFeature.SetField(inLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))

            geom = inFeature.GetGeometryRef()
            if labelType == 'Airplane':
                poly, Length, Width, Aspect, Direction = evaluateLineStringPlane(geom, label='Airplane')
            elif labelType == 'Boat':
                poly, Length, Width, Aspect, Direction = evaluateLineStringBoat(geom, label='Boat')

            outFeature.SetGeometry(poly)
            outFeature.SetField("Length_m", Length)
            outFeature.SetField("Width_m", Width)
            outFeature.SetField("Aspect(L/W)", Aspect)
            outFeature.SetField("compassDeg", Direction)

            outLayer.CreateFeature(outFeature) 
Example 33
Project: rastercube   Author: terrai   File: shputils.py    (license) View Source Project 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 34
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 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 35
Project: pygeotools   Author: dshean   File: geolib.py    (license) View Source Project 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 36
Project: HotSpotAnalysis_Plugin   Author: stanly3690   File: hotspot_analysis.py    (license) View Source Project 4 votes vote down vote up
def load_comboBox(self, layers):
        """Load the fields into combobox when layers are changed"""

        selectedLayerIndex = self.dlg.comboBox.currentIndex()

        if selectedLayerIndex < 0 or selectedLayerIndex > len(layers):
            return
        try:
            selectedLayer = layers[selectedLayerIndex]
        except:
            return

        fieldnames = [field.name() for field in selectedLayer.pendingFields()]

        self.clear_fields()
        self.dlg.comboBox_C.addItems(fieldnames)
        (path, layer_id) = selectedLayer.dataProvider().dataSourceUri().split('|')

        inDriver = ogr.GetDriverByName("ESRI Shapefile")
        inDataSource = inDriver.Open(path, 0)
        inLayer = inDataSource.GetLayer()
        global type
        type = inLayer.GetLayerDefn().GetGeomType()

        if type == 3:  # is a polygon
            self.dlg.checkBox_queen.setChecked(True)
            self.dlg.lineEditThreshold.setEnabled(False)
            self.dlg.checkBox_optimizeDistance.setChecked(False)
            self.dlg.checkBox_optimizeDistance.setEnabled(False)
            self.dlg.lineEdit_minT.setEnabled(False)
            self.dlg.lineEdit_maxT.setEnabled(False)
            self.dlg.lineEdit_dist.setEnabled(False)

        else:
            self.dlg.checkBox_queen.setChecked(False)
            self.dlg.lineEditThreshold.setEnabled(True)
            self.dlg.checkBox_optimizeDistance.setEnabled(True)
            self.dlg.lineEdit_minT.setEnabled(True)
            self.dlg.lineEdit_dist.setEnabled(True)
            self.dlg.lineEdit_maxT.setEnabled(True)
            thresh = pysal.min_threshold_dist_from_shapefile(path)
            self.dlg.lineEditThreshold.setText(str(int(thresh))) 
Example 37
Project: gml_application_schema_toolbox   Author: BRGM   File: load_gml_as_xml.py    (license) View Source Project 4 votes vote down vote up
def _create_layer(self, type, srid, attributes, title):
        """
        Creates an empty spatialite layer
        :param type: 'Point', 'LineString', 'Polygon', etc.
        :param srid: CRS ID of the layer
        :param attributes: list of (attribute_name, attribute_type, attribute_typename)
        :param title: title of the layer
        """
        driver = ogr.GetDriverByName('GPKG')
        ds = driver.CreateDataSource(self.output_local_file)
        layer = ds.CreateLayer("meta", geom_type = ogr.wkbNone)
        layer.CreateField(ogr.FieldDefn('key', ogr.OFTString))
        layer.CreateField(ogr.FieldDefn('value', ogr.OFTString))

        if srid:
            wkbType = { 'point': ogr.wkbPoint,
                        'multipoint': ogr.wkbMultiPoint,
                        'linestring': ogr.wkbLineString,
                        'multilinestring': ogr.wkbMultiLineString,
                        'polygon': ogr.wkbPolygon,
                        'multipolygon': ogr.wkbMultiPolygon
            }[type]
            srs = osr.SpatialReference()
            srs.ImportFromEPSGA(int(srid))
        else:
            wkbType = ogr.wkbNone
            srs = None
        layer = ds.CreateLayer("data", srs, wkbType, ['FID=id'])
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64))
        layer.CreateField(ogr.FieldDefn('fid', ogr.OFTString))
        layer.CreateField(ogr.FieldDefn('_xml_', ogr.OFTString))

        att_type_map = {QVariant.String : ogr.OFTString,
                        QVariant.Int : ogr.OFTInteger,
                        QVariant.Double: ogr.OFTReal,
                        QVariant.DateTime: ogr.OFTDateTime}
        for aname, atype in attributes:
            layer.CreateField(ogr.FieldDefn(aname, att_type_map[atype]))

        # update fields
        layer.ResetReading()

        qgs_layer = QgsVectorLayer("{}|layername=data".format(self.output_local_file), title, "ogr")
        return qgs_layer 
Example 38
Project: wradlib   Author: wradlib   File: gdal.py    (license) View Source Project 4 votes vote down vote up
def gdal_create_dataset(drv, name, cols=0, rows=0, bands=0,
                        gdal_type=gdal.GDT_Unknown, remove=False):
    """Creates GDAL.DataSet object.

    .. versionadded:: 0.7.0

    .. versionchanged:: 0.11.0
        - changed parameters to keyword args
        - added 'bands' as parameter

    Parameters
    ----------
    drv : string
        GDAL driver string
    name : string
        path to filename
    cols : int
        # of columns
    rows : int
        # of rows
    bands : int
        # of raster bands
    gdal_type : raster data type
        eg. gdal.GDT_Float32
    remove : bool
        if True, existing gdal.Dataset will be
        removed before creation

    Returns
    -------
    out : gdal.Dataset
        object

    """
    driver = gdal.GetDriverByName(drv)
    metadata = driver.GetMetadata()

    if not metadata.get('DCAP_CREATE', False):
        raise IOError("Driver %s doesn't support Create() method.".format(drv))

    if remove:
        if os.path.exists(name):
            driver.Delete(name)
    ds = driver.Create(name, cols, rows, bands, gdal_type)

    return ds 
Example 39
Project: dxf2gmlcatastro   Author: sigdeletras   File: dxf2gmlcatastro.py    (license) View Source Project 4 votes vote down vote up
def crea_gml(dxf_origen_file, gml_salida_file, src):
	""" Transforma la información de la geometría de un archivo DXF al estándar de Catastro en formato GML.

	:dxf_origen_file:   Dirección del archivo en formato DXF con la geometría de origen
	:gml_salida_file:   Dirección del archivo en formato GML a sobreescribir con el resultado
	:src:               Sistema de Refencia de Coordendas del DXF. Según cógigos  EPSG
	"""
	# Accede mediante GDAL al archivo DXF
	driver = ogr.GetDriverByName('DXF')
	data_source = driver.Open(dxf_origen_file, 0)
	layer = data_source.GetLayer()

	if src not in SRC_DICT:     # Comprueba que el SRC es correcto
		print('ERROR: El código SRC ({}) indicado es incorrecto.'.format(src))
		print('Los SRC permitidos son 25828, 25829, 25830 o 25831')

		sys.exit()

	print('Archivo GML de salida: {}'.format(gml_salida_file))
	print('Código EPSG seleccionado: {}\n'.format(src))

	with open(gml_salida_file, 'w') as filegml:
		filegml.writelines(PLANTILLA_1)

		print('El archivo {} contiene {} geometría.'.format(dxf_origen_file, layer.GetFeatureCount()))

		for feature in layer:
			geom = feature.GetGeometryRef()
			
			area = geom.Area()
			print('El área del polígono es {:.4f} m2.'.format(area))
			
			filegml.writelines(str(area))       # Añade el área al GML
			
			perimetro = geom.GetGeometryRef(0)

			print('\nTotal de vértices del polígono: {}'.format(perimetro.GetPointCount()))
			print('Listado de coordenadas de los vértices:\nid,x,y')

			filegml.writelines(PLANTILLA_2.format(src=src))

			for i in range(0, perimetro.GetPointCount()):
				pt = perimetro.GetPoint(i)
				coordlist = [str(pt[0]), ' ', str(pt[1]), '\n']
				
				filegml.writelines(coordlist)       # Añade listado de coordenadas X e Y al GML
				
				print('{},{:.4f},{:.4f}'.format(i, pt[0], pt[1]))

		filegml.writelines(PLANTILLA_3)     # Añade XML 
Example 40
Project: antares   Author: CONABIO   File: split.py    (license) View Source Project 4 votes vote down vote up
def split_shape_into_features(shape_name, destination_directory, column, name, sep):
    '''
    This method will take a input shape and iterate over its features, creating
    a new shape file with each one of them. It copies all the fields and the
    same spatial reference from the original file. The created files are saved
    in the destination directory using the number of the field given. 
    '''
    driver = ogr.GetDriverByName(str('ESRI Shapefile'))
    shape = driver.Open(shape_name, 0)
    layer = shape.GetLayer()
    layer_name = layer.GetName()
    spatial_reference = layer.GetSpatialRef()
    in_feature = layer.GetNextFeature()
    shape_files = []
    
    while in_feature:
        encoding = 'utf-8'
        in_feature_name = in_feature.GetField(column)
        in_name = create_filename_from_string(in_feature.GetField(name).decode(encoding))

        final_path = destination_directory + str(in_feature.GetField(0))
        create_directory_path(final_path)
        output_name = create_filename(final_path, '%s__%s.shp' % (in_feature_name, in_name))
        shape_files.append(output_name)

        if os.path.exists(output_name):
            driver.DeleteDataSource(output_name)
        datasource = driver.CreateDataSource(output_name)

        outLayer = datasource.CreateLayer(layer_name, spatial_reference, geom_type=ogr.wkbPolygon)

        inLayerDefn = layer.GetLayerDefn()
        for i in range(0, inLayerDefn.GetFieldCount()):
            fieldDefn = inLayerDefn.GetFieldDefn(i)
            #LOGGER.debug(fieldDefn.GetName())
            outLayer.CreateField(fieldDefn)
        
        outLayerDefn = outLayer.GetLayerDefn()
        geometry = in_feature.GetGeometryRef()

        out_feature = ogr.Feature(outLayerDefn)
        out_feature.SetGeometry(geometry)

        for i in range(0, outLayerDefn.GetFieldCount()):
            out_feature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))
  
        outLayer.CreateFeature(out_feature)
        out_feature.Destroy()
        in_feature.Destroy()
        in_feature = layer.GetNextFeature()
    shape.Destroy()
    datasource.Destroy()
    return shape_files