Python ogr.Open方法代码示例

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def reproject_modis_to_geotiff(input_folder, dest_prj=32643):
    '''
    Function to convert all modis hdf in folder to geotiff using gdal
    required libraries: osr, gdal
    :param input_folder: where hds files are stored
    :param dest_prj: default is wgs 84 / UTM 43N (EPSG 32643)
    '''
    files_list = os.listdir(input_folder)
    for item in files_list:
        if fnmatch.fnmatch(item, '*.hdf'):
            input_file_name = input_folder + '/' + item
            output_file_name = input_folder + '/' + item[0:23]
            img = gdal.Open(input_file_name)
            subset = img.GetSubDatasets()
            in_raster = subset[0][0]
            reproject_dataset(dataset=in_raster, output_file_name=output_file_name+ '.tif')

# reproject
# reproject_modis_to_geotiff(input_folder="/media/kiruba/New Volume/MODIS/ET/scratch")

# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL. 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def rasterize(data, vectorSrc, field, outFile):
    dataSrc = gdal.Open(data)
    import ogr
    shp = ogr.Open(vectorSrc)

    lyr = shp.GetLayer()

    driver = gdal.GetDriverByName('GTiff')
    dst_ds = driver.Create(
        outFile,
        dataSrc.RasterXSize,
        dataSrc.RasterYSize,
        1,
        gdal.GDT_UInt16)
    dst_ds.SetGeoTransform(dataSrc.GetGeoTransform())
    dst_ds.SetProjection(dataSrc.GetProjection())
    if field is None:
        gdal.RasterizeLayer(dst_ds, [1], lyr, None)
    else:
        OPTIONS = ['ATTRIBUTE=' + field]
        gdal.RasterizeLayer(dst_ds, [1], lyr, None, options=OPTIONS)

    data, dst_ds, shp, lyr = None, None, None, None
    return outFile 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def pointReading(fn,pt_lyrName_r):
    ds=ogr.Open(fn,0) #0为只读模式，1为编辑模式
    if ds is None:
        sys.exit('Could not open{0}'.format(fn))
    pt_lyr=ds.GetLayer(pt_lyrName_r) #可以直接数据层(文件)名或者指定索引
    vp = VectorPlotter(True)  #显示vector数据
    vp.plot(pt_lyr,'bo')
    
    i=0
    for feat in pt_lyr: #循环feature
        pt=feat.geometry()
        pt_x=pt.GetX()
        pt_y=pt.GetY()
        name=feat.GetField('NAME')
        kind=feat.GetField('KIND')
        print(name,kind,(pt_x,pt_y,))
        i+=1
        if i==12:
            break
    del ds 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_stack_images(managed_multiple_geotiff_dir):
    test_dir = managed_multiple_geotiff_dir
    images = [os.path.join(test_dir.path, image) for image in os.listdir(test_dir.path)]
    images.sort()
    result_path = os.path.join(test_dir.path, "test_out.tif") # This isn't turning up in the expected order
    pyeo.raster_manipulation.stack_images(images, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray().shape[0] == 4  # If result has 4 bands, success
    assert result.ReadAsArray()[0,0,4] == 3
    assert result.ReadAsArray()[1,0,4] == 13 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_stack_and_trim_images(managed_noncontiguous_geotiff_dir):
    # Test data is two five band 11x12 pixel geotiffs and a 10x10 polygon
    # The geotiffs upper left corners ar at 90,90 and 100,100
    test_dir = managed_noncontiguous_geotiff_dir
    old_image_path = os.path.join(test_dir.path, "se_test")
    new_image_path = os.path.join(test_dir.path, "ne_test")
    aoi_path = os.path.join(test_dir.path, "aoi")
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.raster_manipulation.stack_and_trim_images(old_image_path, new_image_path, aoi_path, result_path)
    result = gdal.Open(result_path).ReadAsArray()
    assert result.shape == (10, 10, 10)  # Result should have 10 bands and be 10 by 10 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_pixel_bounds_from_polygon(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    layer = aoi.GetLayer(0)
    aoi_feature = layer.GetFeature(0)
    polygon = aoi_feature.GetGeometryRef()
    result = pyeo.coordinate_manipulation.pixel_bounds_from_polygon(raster, polygon)
    assert result == (1, 11, 1, 11) 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_point_to_pixel_coordinates(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    point = r"POINT (101 112)"
    raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    out = pyeo.coordinate_manipulation.point_to_pixel_coordinates(raster, point)
    # Since ne_test is an 11 by 12 image, tl corner 0,0 w 10m pixels,
    # we'd expect the br coords to be 10, 11
    assert out == (10, 11) 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_check_overlap(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    assert pyeo.coordinate_manipulation.check_overlap(test_raster, test_aoi) is True 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_get_raster_bounds(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "ne_test"))
    target = (0, 110, 0, 120)
    result = pyeo.coordinate_manipulation.get_raster_bounds(test_raster)
    assert result.GetEnvelope() == target 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_get_aoi_bounds(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    target = (10, 20, 10, 20)
    result = pyeo.coordinate_manipulation.get_aoi_bounds(test_aoi)
    assert result.GetEnvelope() == target 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_get_aoi_intersection(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_raster = gdal.Open(os.path.join(test_dir.path, "se_test"))
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result = pyeo.coordinate_manipulation.get_aoi_intersection(test_raster, test_aoi)
    assert result.GetEnvelope() == (10, 20, 10, 20) 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def test_aoi_to_mask(managed_noncontiguous_geotiff_dir):
    test_dir = managed_noncontiguous_geotiff_dir
    test_aoi = ogr.Open(os.path.join(test_dir.path, "aoi"))
    result_path = os.path.join(test_dir.path, "test_out.tif")
    pyeo.aoi_to_mask(test_aoi, result_path)
    result = gdal.Open(result_path)
    assert result.ReadAsArray()[2, 2] == 1 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def clip_raster_by_vector(input_folder, mask_shapefile, output_folder, file_extension='*.tif', t_srs='EPSG:32643', no_data=32767 ):
    files_list = os.listdir(input_folder)
    ds = ogr.Open(mask_shapefile)
    lyr = ds.GetLayer(0)
    lyr.ResetReading()
    ft = lyr.GetNextFeature()
    for item in files_list:
        if fnmatch.fnmatch(item, file_extension):
            in_raster = input_folder + '/' + item
            out_raster = output_folder + '/' +'tg_' + item
            subprocess.call(['gdalwarp', in_raster, out_raster, '-cutline', mask_shapefile, '-t_srs', t_srs, '-crop_to_cutline', '-dstnodata', "%s" %no_data]) 

# 需要导入模块: import ogr [as 别名]
# 或者: from ogr import Open [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 ogr [as 别名]
# 或者: from ogr import Open [as 别名]
def process_modis(input_folder, output_folder, scale_factor=0.1, null_value=32760, file_extension='*.tif'):
    files_list = os.listdir(input_folder)
    for item in files_list:
        if fnmatch.fnmatch(item, file_extension):
            in_raster = input_folder + '/' + item
            in_raster = gdal.Open(in_raster, GA_ReadOnly)
            out_raster = output_folder + '/' + 'proc_' + item
            band1 = in_raster.GetRasterBand(1)
            in_array = BandReadAsArray(band1)
            in_array[in_array > null_value] = np.nan
            data_out = in_array * scale_factor
            gdalnumeric.SaveArray(data_out, filename=out_raster, format='GTiff', prototype=in_raster)