Python osr.SpatialReference方法代码示例

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def GetGeoInfo(FileName):

    if exists(FileName) is False:
            raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'')    
    
    
    SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly)
    if SourceDS == None:
        raise Exception("Unable to read the data file")
    
    NDV = SourceDS.GetRasterBand(1).GetNoDataValue()
    xsize = SourceDS.RasterXSize
    ysize = SourceDS.RasterYSize
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
    DataType = SourceDS.GetRasterBand(1).DataType
    DataType = gdal.GetDataTypeName(DataType)
    
    return NDV, xsize, ysize, GeoT, Projection, DataType
#==============================================================================

#==============================================================================
# Function to read the original file's projection: 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def __init__(self, back_ds, wkt_stored, rect, **kwargs):
        wkt_virtual = back_ds.virtual_of_stored_given_mode(
            wkt_stored, back_ds.wkt_work, back_ds.wkt_fallback, back_ds.wkt_forced,
        )

        if wkt_virtual is not None:
            sr_virtual = osr.SpatialReference(wkt_virtual)
        else:
            sr_virtual = None

        to_work, to_virtual = back_ds.get_transforms(sr_virtual, rect)

        self.back_ds = back_ds
        self.wkt_stored = wkt_stored
        self.wkt_virtual = wkt_virtual
        self.to_work = to_work
        self.to_virtual = to_virtual
        super().__init__(**kwargs) 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def __init__(self, wkt_work, wkt_fallback, wkt_forced, analyse_transformation, **kwargs):

        if wkt_work is not None:
            sr_work = osr.SpatialReference(wkt_work)
        else:
            sr_work = None
        if wkt_fallback is not None:
            sr_fallback = osr.SpatialReference(wkt_fallback)
        else:
            sr_fallback = None
        if wkt_forced is not None:
            sr_forced = osr.SpatialReference(wkt_forced)
        else:
            sr_forced = None

        self.wkt_work = wkt_work
        self.wkt_fallback = wkt_fallback
        self.wkt_forced = wkt_forced
        self.sr_work = sr_work
        self.sr_fallback = sr_fallback
        self.sr_forced = sr_forced
        self.analyse_transformations = analyse_transformation
        super(BackDatasetConversionsMixin, self).__init__(**kwargs) 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def reproject(geom, from_epsg, to_epsg):
    """
    Reproject the given geometry from the given EPSG code to another
    """
    # Note: this is currently only accurate for the U.S.
    source = osr.SpatialReference()
    source.ImportFromEPSG(from_epsg)

    target = osr.SpatialReference()
    target.ImportFromEPSG(to_epsg)

    transform = osr.CoordinateTransformation(source, target)

    geom.Transform(transform)

    return geom 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def get_coord_transform(source_epsg, target_epsg):
    '''
    Creates an OGR-framework coordinate transformation for use in projecting
    coordinates to a new coordinate reference system (CRS). Used as, e.g.:
        transform = get_coord_transform(source_epsg, target_epsg)
        transform.TransformPoint(x, y)
    Arguments:
        source_epsg     The EPSG code for the source CRS
        target_epsg     The EPSG code for the target CRS
    '''
    # Develop a coordinate transformation, if desired
    transform = None
    source_ref = osr.SpatialReference()
    target_ref = osr.SpatialReference()
    source_ref.ImportFromEPSG(source_epsg)
    target_ref.ImportFromEPSG(target_epsg)
    return osr.CoordinateTransformation(source_ref, target_ref) 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
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() 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False):
    """
    Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path.
    Assumes both raster_to_clip and extent_raster are in the same projection.
    Parameters
    ----------
    raster_to_clip_path
        The location of the raster to be clipped.
    extent_raster_path
        The location of the raster that will provide the extent to clip to
    out_raster_path
        A location for the finished raster
    """

    with TemporaryDirectory() as td:
        temp_aoi_path = os.path.join(td, "temp_clip.shp")
        get_extent_as_shp(extent_raster_path, temp_aoi_path)
        ext_ras = gdal.Open(extent_raster_path)
        proj = osr.SpatialReference(wkt=ext_ras.GetProjection())
        srs_id = int(proj.GetAttrValue('AUTHORITY', 1))
        clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat) 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max,
                                        y_min):
    # local variables:
    temp_shapefile = locator.get_temporary_file("terrain.shp")
    cols = int((x_max - x_min) / grid_size)
    rows = int((y_max - y_min) / grid_size)
    shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]])
    geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes])
    geodataframe.to_file(temp_shapefile)
    # 1) opening the shapefile
    source_ds = ogr.Open(temp_shapefile)
    source_layer = source_ds.GetLayer()
    target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32)  ##COMMENT 2
    target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size))  ##COMMENT 3
    # 5) Adding a spatial reference ##COMMENT 4
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromProj4(crs)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(-9999)  ##COMMENT 5
    gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation])  ##COMMENT 6
    target_ds = None  # closing the file 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def setup_input_srs(input_dataset, options):
    """
    Determines and returns the Input Spatial Reference System (SRS) as an osr object and as a
    WKT representation

    Uses in priority the one passed in the command line arguments. If None, tries to extract them
    from the input dataset
    """

    input_srs = None
    input_srs_wkt = None

    if options.s_srs:
        input_srs = osr.SpatialReference()
        input_srs.SetFromUserInput(options.s_srs)
        input_srs_wkt = input_srs.ExportToWkt()
    else:
        input_srs_wkt = input_dataset.GetProjection()
        if not input_srs_wkt and input_dataset.GetGCPCount() != 0:
            input_srs_wkt = input_dataset.GetGCPProjection()
        if input_srs_wkt:
            input_srs = osr.SpatialReference()
            input_srs.ImportFromWkt(input_srs_wkt)

    return input_srs, input_srs_wkt 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def toWgs(mgrs):
    """ Converts an MGRS coordinate string to geodetic (latitude and longitude)
    coordinates

    @param mgrs - MGRS coordinate string
    @returns - tuple containning latitude and longitude values
    """
    if _checkZone(mgrs):
        zone, hemisphere, easting, northing = _mgrsToUtm(mgrs)
    else:
        zone, hemisphere, easting, northing = _mgrsToUps(mgrs)

    epsg = _epsgForUtm(zone, hemisphere)
    src = osr.SpatialReference()
    src.ImportFromEPSG(epsg)
    dst = osr.SpatialReference()
    dst.ImportFromEPSG(4326)
    ct = osr.CoordinateTransformation(src, dst)
    longitude, latitude, z = ct.TransformPoint(easting, northing)

    return latitude, longitude 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def test_dump_raster(self):
        proj = osr.SpatialReference()
        proj.ImportFromEPSG(31466)
        filename = util.get_wradlib_data_file("shapefiles/agger/" "agger_merge.shp")
        test = zonalstats.DataSource(filename, srs=proj)
        test.dump_raster(
            tempfile.NamedTemporaryFile(mode="w+b").name,
            driver="netCDF",
            pixel_size=100.0,
        )
        test.dump_raster(
            tempfile.NamedTemporaryFile(mode="w+b").name,
            driver="netCDF",
            pixel_size=100.0,
            attr="FID",
        ) 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def get_radar_projection(sitecoords):
    """Get the native radar projection which is an azimuthal equidistant projection
    centered at the site using WGS84.

    Parameters
    ----------
    sitecoords : a sequence of two floats
        the WGS84 lon / lat coordinates of the radar location

    Returns
    -------
    proj : osr.SpatialReference
        projection definition

    """
    proj = osr.SpatialReference()
    proj.SetProjCS("Unknown Azimuthal Equidistant")
    proj.SetAE(sitecoords[1], sitecoords[0], 0, 0)

    return proj 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def get_extent(coords):
    """Get the extent of 2d coordinates

    Parameters
    ----------
    coords : :class:`numpy:numpy.ndarray`
        coordinates array with shape (...,(x,y))

    Returns
    -------
    proj : osr.SpatialReference
        GDAL/OSR object defining projection
    """

    xmin = coords[..., 0].min()
    xmax = coords[..., 0].max()
    ymin = coords[..., 1].min()
    ymax = coords[..., 1].max()

    return xmin, xmax, ymin, ymax 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def read_gdal_projection(dataset):
    """Get a projection (OSR object) from a GDAL dataset.

    Parameters
    ----------
    dataset : gdal.Dataset
        raster image with georeferencing

    Returns
    -------
    srs : OSR.SpatialReference
        dataset projection object

    Examples
    --------

    See :ref:`/notebooks/classify/wradlib_clutter_cloud_example.ipynb`.

    """
    wkt = dataset.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    # src = None
    return srs 

# 需要导入模块: from osgeo import osr [as 别名]
# 或者: from osgeo.osr import SpatialReference [as 别名]
def write_geotiff(tiff_path, data, upper_left_lon, upper_left_lat, step, AREA_OR_POINT='Point'):
    """
    写入geotiff。默认pixel is point，默认WGS84
    """
    rows, columns = data.shape
    bands = 1
    pixel_width = step
    pixel_height = -step
    half_step = step / 2
    upper_left_x = upper_left_lon - half_step
    upper_left_y = upper_left_lat + half_step
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(tiff_path, columns, rows, bands, gdal.GDT_Float32)
    dataset.SetMetadataItem('AREA_OR_POINT', AREA_OR_POINT)
    dataset.SetGeoTransform([upper_left_x, pixel_width, 0,
                             upper_left_y, 0, pixel_height])
    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')
    dataset.SetProjection(sr.ExportToWkt())  # 给新建图层赋予投影信息
    dataset.GetRasterBand(1).WriteArray(data)
    dataset.FlushCache()  # 将数据写入硬盘
    dataset = None