Python ogr.GetDriverByName方法代碼示例

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def create_temp_tiff(self, name, content=np.ones([3, 3, 3]), geotransform=(10, 10, 0, 10, 0, -10)):
        """Creates a temporary geotiff in self.path
        """
        if len(content.shape) != 3:
            raise IndexError
        path = os.path.join(self.path, name)
        driver = gdal.GetDriverByName('GTiff')
        new_image = driver.Create(
            path,
            xsize=content.shape[1],
            ysize=content.shape[2],
            bands=content.shape[0],
            eType=gdal.GDT_Byte
        )
        new_image.SetGeoTransform(geotransform)
        for band in range(content.shape[0]):
            raster_band = new_image.GetRasterBand(band+1)
            raster_band.WriteArray(content[band, ...].T)
        new_image.SetProjection(self.srs.ExportToWkt())
        new_image.FlushCache()
        self.images.append(new_image)
        self.image_paths.append(path) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def create_temp_shape(self, name, point_list):
        vector_file = os.path.join(self.temp_dir.name, name)
        shape_driver = ogr.GetDriverByName("ESRI Shapefile")  # Depreciated; replace at some point
        vector_data_source = shape_driver.CreateDataSource(vector_file)
        vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for point in point_list:
            ring.AddPoint(point[0], point[1])
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        vector_feature_definition = vector_layer.GetLayerDefn()
        vector_feature = ogr.Feature(vector_feature_definition)
        vector_feature.SetGeometry(poly)
        vector_layer.CreateFeature(vector_feature)
        vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
        feature = ogr.Feature(vector_layer.GetLayerDefn())
        feature.SetField("class", 3)

        vector_data_source.FlushCache()
        self.vectors.append(vector_data_source)  # Check this is the right thing to be saving here
        self.vector_paths.append(vector_file) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def geotiff_dir():
    """

    Returns
    -------
    A pointer to a temporary folder that contains a 3-band geotiff
    of 3x3, with all values being 1.

    """
    tempDir = tempfile.TemporaryDirectory()
    fileformat = "GTiff"
    driver = gdal.GetDriverByName(fileformat)
    metadata = driver.GetMetadata()
    tempPath = os.path.join(tempDir.name)
    testDataset = driver.Create(os.path.join(tempDir.name, "tempTiff.tif"),
                                xsize=3, ysize=3, bands=3, eType=gdal.GDT_CFloat32)
    for i in range(3):
        testDataset.GetRasterBand(i+1).WriteArray(np.ones([3, 3]))
    testDataset = None
    yield tempPath
    tempDir.cleanup() 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Get_Extent(input_lyr):
    """
    Obtain the input layer extent (xmin, ymin, xmax, ymax)
    """
    # Input
    filename, ext = os.path.splitext(input_lyr)
    if ext.lower() == '.shp':
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr)
        inp_lyr = inp_source.GetLayer()
        x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
        inp_lyr = None
        inp_source = None
    elif ext.lower() == '.tif':
        inp_lyr = gdal.Open(input_lyr)
        inp_transform = inp_lyr.GetGeoTransform()
        x_min = inp_transform[0]
        x_max = x_min + inp_transform[1] * inp_lyr.RasterXSize
        y_max = inp_transform[3]
        y_min = y_max + inp_transform[5] * inp_lyr.RasterYSize
        inp_lyr = None
    else:
        raise Exception('The input data type is not recognized')
    return (x_min, y_min, x_max, y_max) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def create_100x100_shp(self, name):
        """Cretes  a shapefile with a vector layer named "geometry" containing a 100mx100m square , top left corner
        being at wgs coords 10,10.
        This polygon has a field, 'class' with a value of 3. Left in for back-compatability"""
        # TODO Generalise this
        vector_file = os.path.join(self.temp_dir.name, name)
        shape_driver = ogr.GetDriverByName("ESRI Shapefile")  # Depreciated; replace at some point
        vector_data_source = shape_driver.CreateDataSource(vector_file)
        vector_layer = vector_data_source.CreateLayer("geometry", self.srs, geom_type=ogr.wkbPolygon)
        ring = ogr.Geometry(ogr.wkbLinearRing)
        ring.AddPoint(10.0, 10.0)
        ring.AddPoint(10.0, 110.0)
        ring.AddPoint(110.0, 110.0)
        ring.AddPoint(110.0, 10.0)
        ring.AddPoint(10.0, 10.0)
        poly = ogr.Geometry(ogr.wkbPolygon)
        poly.AddGeometry(ring)
        vector_feature_definition = vector_layer.GetLayerDefn()
        vector_feature = ogr.Feature(vector_feature_definition)
        vector_feature.SetGeometry(poly)
        vector_layer.CreateFeature(vector_feature)
        vector_layer.CreateField(ogr.FieldDefn("class", ogr.OFTInteger))
        feature = ogr.Feature(vector_layer.GetLayerDefn())
        feature.SetField("class", 3)

        vector_data_source.FlushCache()
        self.vectors.append(vector_data_source)  # Check this is the right thing to be saving here
        self.vector_paths.append(vector_file) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
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

        # options = ['SPATIALITE=YES'])
        ds = outDriver.CreateDataSource(outShapeFile)

        # 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 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def List_Fields(input_lyr):
    """
    Lists the field names of input layer
    """
    # Input
    if isinstance(input_lyr, str):
        inp_driver = ogr.GetDriverByName('ESRI Shapefile')
        inp_source = inp_driver.Open(input_lyr, 0)
        inp_lyr = inp_source.GetLayer()
        inp_lyr_defn = inp_lyr.GetLayerDefn()
    elif isinstance(input_lyr, ogr.Layer):
        inp_lyr_defn = input_lyr.GetLayerDefn()

    # List
    names_ls = []

    # Loop
    for j in range(0, inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(j)
        names_ls.append(field_defn.GetName())

    # Save and/or close the data sources
    inp_source = None

    # Return
    return names_ls 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Array_to_Raster(input_array, output_tiff, ll_corner, cellsize,
                    srs_wkt):
    """
    Saves an array into a raster file
    """
    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    y_ncells, x_ncells = input_array.shape
    gdal_datatype = gdaltype_from_dtype(input_array.dtype)

    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal_datatype)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(-9999)

    out_top_left_x = ll_corner[0]
    out_top_left_y = ll_corner[1] + cellsize*y_ncells

    out_source.SetGeoTransform((out_top_left_x, cellsize, 0,
                                out_top_left_y, 0, -cellsize))
    out_source.SetProjection(str(srs_wkt))
    out_band.WriteArray(input_array)

    # Save and/or close the data sources
    out_source = None

    # Return
    return output_tiff 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Extract_Band(input_tiff, output_tiff, band_number=1):
    """
    Extract and save a raster band into a new raster
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(band_number)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform(inp_transform)
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(inp_array)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def get_layer_names(geo_fpath, selected_layer_name=None):
    """
    Since layers in the same dataset can use
    different projections, we'll need to create a
    DbEntity for each layer. Use the osr library from
    gdal to loop through layers and grab their names.
    """

    _, geo_ext = os.path.splitext(geo_fpath)
    driver_name = OGR_DRIVER_MAP[geo_ext]

    driver = ogr.GetDriverByName(driver_name)
    dataset = driver.Open(geo_fpath, 0)

    layer_names = []
    for idx in xrange(dataset.GetLayerCount()):
        layer = dataset.GetLayerByIndex(idx)
        layer_name = layer.GetName()

        if selected_layer_name:
            if selected_layer_name == layer_name:
                layer_names.append(layer_name)
                break
        else:
            layer_names.append(layer_name)

    return layer_names 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def save_point_list_to_shapefile(class_sample_point_dict, out_path, geotransform, projection_wkt, produce_csv=False):
    """Saves a list of points to a shapefile at out_path. Need the gt and projection of the raster.
    GT is needed to move each point to the centre of the pixel. Can also produce a .csv file for CoolEarth"""
    log = logging.getLogger(__name__)
    log.info("Saving point list to shapefile")
    log.debug("GT: {}\nProjection: {}".format(geotransform, projection_wkt))
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection_wkt)
    layer = data_source.CreateLayer("validation_points", srs, ogr.wkbPoint)
    class_field = ogr.FieldDefn("class", ogr.OFTString)
    class_field.SetWidth(24)
    layer.CreateField(class_field)

    for map_class, point_list in class_sample_point_dict.items():
        for point in point_list:
            feature = ogr.Feature(layer.GetLayerDefn())
            coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
            offset = geotransform[1]/2   # Adds half a pixel offset so points end up in the center of pixels
            wkt = "POINT({} {})".format(coord[0]+offset, coord[1]-offset) # Never forget about negative y values in gts.
            new_point = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(new_point)
            feature.SetField("class", map_class)
            layer.CreateFeature(feature)
            feature = None

    layer = None
    data_source = None

    if produce_csv:
        csv_out_path = out_path.rsplit('.')[0] + ".csv"
        with open(csv_out_path, "w") as csv_file:
            writer = csv.writer(csv_file)
            writer.writerow(["id", "yCoordinate", "xCoordinate"])

            # Join all points create single dimesional list of points (and revise the '*' operator)
            for id,  point in enumerate(itertools.chain(*class_sample_point_dict.values())):
                coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
                offset = geotransform[1] / 2  # Adds half a pixel offset so points end up in the center of pixels
                lat = coord[0] + offset
                lon = coord[1] - offset
                writer.writerow([id, lon, lat])
        log.info("CSV out at: {}".format(csv_out_path)) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Buffer(input_shp, output_shp, distance):
    """
    Creates a buffer of the input shapefile by a given distance
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_lyr_defn = inp_lyr.GetLayerDefn()
    inp_srs = inp_lyr.GetSpatialRef()

    # Output
    out_name = os.path.splitext(os.path.basename(output_shp))[0]
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(output_shp):
        out_driver.DeleteDataSource(output_shp)
    out_source = out_driver.CreateDataSource(output_shp)

    out_lyr = out_source.CreateLayer(out_name, inp_srs, ogr.wkbPolygon)
    out_lyr_defn = out_lyr.GetLayerDefn()

    # Add fields
    for i in range(inp_lyr_defn.GetFieldCount()):
        field_defn = inp_lyr_defn.GetFieldDefn(i)
        out_lyr.CreateField(field_defn)

    # Add features
    for i in range(inp_lyr.GetFeatureCount()):
        feature_inp = inp_lyr.GetNextFeature()
        geometry = feature_inp.geometry()
        feature_out = ogr.Feature(out_lyr_defn)

        for j in range(0, out_lyr_defn.GetFieldCount()):
            feature_out.SetField(out_lyr_defn.GetFieldDefn(j).GetNameRef(),
                                 feature_inp.GetField(j))

        feature_out.SetGeometry(geometry.Buffer(distance))
        out_lyr.CreateFeature(feature_out)
        feature_out = None

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_shp 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Feature_to_Raster(input_shp, output_tiff,
                      cellsize, field_name=False, NoData_value=-9999):
    """
    Converts a shapefile into a raster
    """

    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer()
    inp_srs = inp_lyr.GetSpatialRef()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, gdal.GDT_Int16)

    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr,
                            options=["ATTRIBUTE={0}".format(field_name)])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])

    # Save and/or close the data sources
    inp_source = None
    out_source = None

    # Return
    return output_tiff 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Raster_to_Array(input_tiff, ll_corner, x_ncells, y_ncells,
                    values_type='float32'):
    """
    Loads a raster into a numpy array
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_data_type = inp_band.DataType

    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    ll_x = ll_corner[0]
    ll_y = ll_corner[1]

    top_left_x = ll_x
    top_left_y = ll_y - cellsize_y*y_ncells

    # Change start point
    temp_path = tempfile.mkdtemp()
    temp_driver = gdal.GetDriverByName('GTiff')
    temp_tiff = os.path.join(temp_path, os.path.basename(input_tiff))
    temp_source = temp_driver.Create(temp_tiff, x_ncells, y_ncells,
                                     1, inp_data_type)
    temp_source.GetRasterBand(1).SetNoDataValue(NoData_value)
    temp_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                 top_left_y, rot_2, cellsize_y))
    temp_source.SetProjection(inp_srs)

    # Snap
    gdal.ReprojectImage(inp_lyr, temp_source, inp_srs, inp_srs,
                        gdal.GRA_Bilinear)
    temp_source = None

    # Read array
    d_type = pd.np.dtype(values_type)
    out_lyr = gdal.Open(temp_tiff)
    array = out_lyr.ReadAsArray(0, 0, out_lyr.RasterXSize,
                                out_lyr.RasterYSize).astype(d_type)
    array[pd.np.isclose(array, NoData_value)] = pd.np.nan
    out_lyr = None

    return array 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import GetDriverByName [as 別名]
def Apply_Filter(input_tiff, output_tiff, number_of_passes):
    """
    Smooth a raster by replacing cell value by the average value of the
    surrounding cells
    """
    # Input
    inp_lyr = gdal.Open(input_tiff)
    inp_srs = inp_lyr.GetProjection()
    inp_transform = inp_lyr.GetGeoTransform()
    inp_band = inp_lyr.GetRasterBand(1)
    inp_array = inp_band.ReadAsArray()
    inp_data_type = inp_band.DataType

    top_left_x = inp_transform[0]
    cellsize_x = inp_transform[1]
    rot_1 = inp_transform[2]
    top_left_y = inp_transform[3]
    rot_2 = inp_transform[4]
    cellsize_y = inp_transform[5]
    NoData_value = inp_band.GetNoDataValue()

    x_ncells = inp_lyr.RasterXSize
    y_ncells = inp_lyr.RasterYSize

    # Filter
    inp_array[inp_array == NoData_value] = pd.np.nan
    out_array = array_filter(inp_array, number_of_passes)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,
                                   1, inp_data_type)
    out_band = out_source.GetRasterBand(1)
    out_band.SetNoDataValue(NoData_value)
    out_source.SetGeoTransform((top_left_x, cellsize_x, rot_1,
                                top_left_y, rot_2, cellsize_y))
    out_source.SetProjection(inp_srs)
    out_band.WriteArray(out_array)

    # Save and/or close the data sources
    inp_lyr = None
    out_source = None

    # Return
    return output_tiff