Python osr.SpatialReference方法代碼示例

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array, nodata, EPSG):
    """This function take a regular array to create a raster with it"""
    print("I am dealing with nodata values")
    array[np.isnan(array)] = nodata # Dealing with Nodata values
    print("I am writing the raster")
    cols = array.shape[1]
    rows = array.shape[0]

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('ENVI')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float64)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    #outband.SetNoDataValue(nodata)
    outband.WriteArray(array)
    #outband.SetNoDataValue(nodata)

    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(EPSG)

    outRaster.SetProjection(outRasterSRS.ExportToWkt())

    outband.FlushCache() 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def loadMetadata(self):
        # open the raster and its spatial reference
        self.src = gdal.Open(self.tif_path)

        if self.src is None:
            raise Exception('Could not load GDAL file "%s"' % self.tif_path)
        spatial_reference_raster = osr.SpatialReference(self.src.GetProjection())

        # get the WGS84 spatial reference
        spatial_reference = osr.SpatialReference()
        spatial_reference.ImportFromEPSG(4326)  # WGS84

        # coordinate transformation
        self.coordinate_transform = osr.CoordinateTransformation(spatial_reference, spatial_reference_raster)
        gt = self.geo_transform = self.src.GetGeoTransform()
        dev = (gt[1] * gt[5] - gt[2] * gt[4])
        self.geo_transform_inv = (gt[0], gt[5] / dev, -gt[2] / dev,
                                  gt[3], -gt[4] / dev, gt[1] / dev) 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def to_wkt(self, target_wkt):
        # If we're going from WGS84 -> Spherical Mercator, use PyProj, because
        #  there seems to be a bug in OGR that gives us an offset. (GDAL
        #  does fine, though.
        if target_wkt == self.wkt:
            return self
        import osr
        dstSpatialRef = osr.SpatialReference()
        dstSpatialRef.ImportFromEPSG(d_name_to_epsg[target_wkt])
        # dstSpatialRef.ImportFromWkt(d_name_to_wkt[target_wkt])
        srcSpatialRef = osr.SpatialReference()
        srcSpatialRef.ImportFromEPSG(d_name_to_epsg[self.wkt])
        # srcSpatialRef.ImportFromWkt(self.wkt_)
        coordTransform = osr.CoordinateTransformation(srcSpatialRef, dstSpatialRef)
        a, b, c = coordTransform.TransformPoint(self.lon, self.lat)
        return Point(b, a, wkt=target_wkt) 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def getproj(self, withgrid=False, projformat='pyproj'):
        from PseudoNetCDF.coordutil import getproj4_from_cf_var
        gridmapping = self.variables['crs']
        proj4str = getproj4_from_cf_var(gridmapping, withgrid=withgrid)
        preserve_units = withgrid
        if projformat == 'proj4':
            return proj4str
        elif projformat == 'pyproj':
            import pyproj
            # pyproj adds +units=m, which is not right for latlon/lonlat
            if '+proj=lonlat' in proj4str or '+proj=latlon' in proj4str:
                preserve_units = True
            return pyproj.Proj(proj4str, preserve_units=preserve_units)
        elif projformat == 'wkt':
            import osr
            srs = osr.SpatialReference()
            # Imports WKT to Spatial Reference Object
            srs.ImportFromProj4(proj4str)
            srs.ExportToWkt()  # converts the WKT to an ESRI-compatible format
            return srs.ExportToWkt()
        else:
            raise ValueError('projformat must be pyproj, proj4 or wkt') 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def get_epsg(prjfile):
    '''Get the epsg code from a projection file of a shapefile

    Args:
        prjfile: a .prj file of a shapefile

    Returns:
        str: EPSG code

    '''

    prj_file = open(prjfile, 'r')
    prj_txt = prj_file.read()
    srs = osr.SpatialReference()
    srs.ImportFromESRI([prj_txt])
    srs.AutoIdentifyEPSG()
    # return EPSG code
    return srs.GetAuthorityCode(None) 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def Save_as_MEM(data='', geo='', projection=''):
    """
    This function save the array as a memory file

    Keyword arguments:
    data -- [array], dataset of the geotiff
    geo -- [minimum lon, pixelsize, rotation, maximum lat, rotation,
            pixelsize], (geospatial dataset)
    projection -- interger, the EPSG code
    """
    # save as a geotiff
    driver = gdal.GetDriverByName("MEM")
    dst_ds = driver.Create('', int(data.shape[1]), int(data.shape[0]), 1,
                           gdal.GDT_Float32, ['COMPRESS=LZW'])
    srse = osr.SpatialReference()
    if projection == '':
        srse.SetWellKnownGeogCS("WGS84")
    else:
        srse.SetWellKnownGeogCS(projection)
    dst_ds.SetProjection(srse.ExportToWkt())
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.SetGeoTransform(geo)
    dst_ds.GetRasterBand(1).WriteArray(data)
    return(dst_ds) 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def make_geotiff(data, length, width, latn_p, lonw_p, dlat, dlon, outfile, compress_option):

    if data.dtype == np.float32:
        dtype = gdal.GDT_Float32
        nodata = np.nan  ## or 0?
    elif data.dtype == np.uint8:
        dtype = gdal.GDT_Byte
        nodata = None

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(outfile, width, length, 1, dtype, options=compress_option)
    outRaster.SetGeoTransform((lonw_p, dlon, 0, latn_p, 0, dlat))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(data)
    if nodata is not None: outband.SetNoDataValue(nodata)
    outRaster.SetMetadataItem('AREA_OR_POINT', 'Point')
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

    return


#%% Main 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def get_pixel_latlon(raster, x, y):
    """For a given pixel in raster, gets the lat-lon value in EPSG 4326."""
    # TODO: Move to coordinate_manipulation
    native_projection = osr.SpatialReference()
    native_projection.ImportFromWkt(raster.GetProjection())
    latlon_projection = osr.SpatialReference()
    latlon_projection.ImportFromEPSG(4326)
    transformer = osr.CoordinateTransformation(native_projection, latlon_projection)

    geotransform = raster.GetGeoTransform()
    x_geo, y_geo = cm.pixel_to_point_coordinates([y,x], geotransform)  # Why did I do this reverse?
    lon, lat, _ = transformer.TransformPoint(x_geo, y_geo)
    return lat, lon 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def _generate_latlon_transformer(raster):
    native_projection = osr.SpatialReference()
    native_projection.ImportFromWkt(raster.GetProjection())
    latlon_projection = osr.SpatialReference()
    latlon_projection.ImportFromEPSG(4326)
    geotransform = raster.GetGeoTransform()
    return osr.CoordinateTransformation(native_projection, latlon_projection), geotransform 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def test_composite_across_projections():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove(r"test_outputs/composite_test.tif")
    except FileNotFoundError:
        pass
    try:
        shutil.rmtree(r"test_outputs/reprojected")
    except FileNotFoundError:
        pass
    os.mkdir(r"test_outputs/reprojected")
    epsg = 4326
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(epsg)
    projection = proj.ExportToWkt() # Refactor this terrible nonsense later
    test_data = [r"test_data/S2A_MSIL2A_20180703T073611_N0206_R092_T36MZE_20180703T094637.tif",
                 r"test_data/S2B_MSIL2A_20180728T073609_N0206_R092_T37MBV_20180728T114325.tif"]
    pyeo.raster_manipulation.reproject_image(test_data[0], r"test_outputs/reprojected/0.tif", projection)
    pyeo.raster_manipulation.reproject_image(test_data[1], r"test_outputs/reprojected/1.tif", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[0]), r"test_outputs/reprojected/0.msk", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[1]), r"test_outputs/reprojected/1.msk", projection)

    out_file = r"test_outputs/composite_test.tif"
    pyeo.raster_manipulation.composite_images_with_mask([
       r"test_outputs/reprojected/0.tif",
        r"test_outputs/reprojected/1.tif"],
        out_file)
    image = gdal.Open("test_outputs/composite_test.tif")
    assert image
    image_array = image.GetVirtualMemArray()
    assert image_array.max() > 10 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def test_composite_across_projections_meters():
    os.chdir(os.path.dirname(os.path.abspath(__file__)))
    try:
        os.remove(r"test_outputs/composite_test.tif")
    except FileNotFoundError:
        pass
    try:
        shutil.rmtree(r"test_outputs/reprojected")
    except FileNotFoundError:
        pass
    os.mkdir(r"test_outputs/reprojected")
    epsg = 32736
    proj = osr.SpatialReference()
    proj.ImportFromEPSG(epsg)
    projection = proj.ExportToWkt() # Refactor this terrible nonsense later

    test_data = [r"test_data/S2A_MSIL2A_20180703T073611_N0206_R092_T36MZE_20180703T094637.tif",
                 r"test_data/S2B_MSIL2A_20180728T073609_N0206_R092_T37MBV_20180728T114325.tif"]
    pyeo.raster_manipulation.reproject_image(test_data[0], r"test_outputs/reprojected/0.tif", projection)
    pyeo.raster_manipulation.reproject_image(test_data[1], r"test_outputs/reprojected/1.tif", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[0]), r"test_outputs/reprojected/0.msk", projection)
    pyeo.raster_manipulation.reproject_image(pyeo.filesystem_utilities.get_mask_path(test_data[1]), r"test_outputs/reprojected/1.msk", projection)

    out_file = r"test_outputs/composite_test.tif"
    pyeo.raster_manipulation.composite_images_with_mask([
       r"test_outputs/reprojected/0.tif",
        r"test_outputs/reprojected/1.tif"],
        out_file)
    image = gdal.Open("test_outputs/composite_test.tif")
    assert image
    image_array = image.GetVirtualMemArray()
    assert image_array.max() > 1 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def __init__(self, srs=4326):
        self.srs = osr.SpatialReference()
        self.srs.ImportFromEPSG(srs) 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def reproject_dataset(dataset, output_file_name, wkt_from="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs", epsg_to=32643, pixel_spacing=1000, file_format='GTiff'):
    '''
    :param dataset: Modis subset name use gdal.GetSubdatasets()
    :param output_file_name: file location along with proper extension
    :param wkt_from: Modis wkt (default)
    :param epsg_to: integer default(32643)
    :param pixel_spacing: in metres
    :param file_format: default GTiff
    :return:
    '''
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(epsg_to)
    modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected" 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def raster_to_WSG_and_UTM(raster_path, lat, lon):

    raster = gdal.Open(raster_path)
    source_projection_wkt = raster.GetProjection()
    inSRS_converter = osr.SpatialReference()
    inSRS_converter.ImportFromProj4(get_projected_coordinate_system(lat, lon))
    target_projection_wkt = inSRS_converter.ExportToWkt()
    new_raster = gdal.AutoCreateWarpedVRT(raster, source_projection_wkt, target_projection_wkt,
                                          gdal.GRA_NearestNeighbour)
    return new_raster 

# 需要導入模塊: import osr [as 別名]
# 或者: from osr import SpatialReference [as 別名]
def mk_geotiff_obj(raster, fn, bands=1, gdal_data_type=gdal.GDT_Float32,
                   lat=[46, 45], lon=[-73, -72]):
    """
    Creates a new geotiff file objects using the WGS84 coordinate system, saves
    it to disk, and returns a handle to the python file object and driver

    Parameters
    ------------
    raster : array
        Numpy array of the raster data to be added to the object
    fn : str
        Name of the geotiff file
    bands : int (optional)
        See :py:func:`gdal.GetDriverByName('Gtiff').Create
    gdal_data : gdal.GDT_<type>
        Gdal data type (see gdal.GDT_...)
    lat : list
        northern lat, southern lat
    lon : list
        [western lon, eastern lon]
    """
    NNi, NNj = raster.shape
    driver = gdal.GetDriverByName('GTiff')
    obj = driver.Create(fn, NNj, NNi, bands, gdal_data_type)
    pixel_height = -np.abs(lat[0] - lat[1]) / (NNi - 1.0)
    pixel_width = np.abs(lon[0] - lon[1]) / (NNj - 1.0)
    obj.SetGeoTransform([lon[0], pixel_width, 0, lat[0], 0, pixel_height])
    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS('WGS84')
    obj.SetProjection(srs.ExportToWkt())
    obj.GetRasterBand(1).WriteArray(raster)
    return obj, driver