美文网首页
RasterizeLayer

RasterizeLayer

作者: hehehehe | 来源:发表于2021-01-26 14:35 被阅读0次
    #!/usr/bin/env python
    # -*- coding: utf-8 -*-
    
    from osgeo import ogr
    from osgeo import gdal
    
    # set pixel size
    pixel_size = 0.002
    no_data_value = -9999
    image_type = 'GTiff'
    
    # Shapefile input name
    # input projection must be in cartesian system in meters
    # input wgs 84 or EPSG: 4326 will NOT work!!!
    input_shp = "F:\\builddata\\mkt_wuh\\hzw\\12345\\12_AB_DATA\\1209_AB_SEM_DES\\AD_Arrow2.geojson"
    # TIF Raster file to be created
    output_raster = "F:\\builddata\\mkt_wuh\\hzw\\12345\\12_AB_DATA\\1209_AB_SEM_DES\\AD_Arrow.tif"
    
    # Open the data source get the layer object
    # assign extent coordinates
    #shp_driver = ogr.GetDriverByName('ESRI Shapefile')
    open_shp = ogr.Open(input_shp)
    shp_layer = open_shp.GetLayer()
    x_min, x_max, y_min, y_max = shp_layer.GetExtent()
    
    print(x_min)
    print(x_max)
    
    # calculate raster resolution
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)
    
    print(x_res)
    print(y_res)
    # set the image type for expor
    
    driver = gdal.GetDriverByName(image_type)
    
    new_raster = driver.Create(output_raster, x_res, y_res, 1, gdal.GDT_Byte)
    new_raster.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
    
    # get the raster band we want to export too
    raster_band = new_raster.GetRasterBand(1)
    
    # assign the no data value to empty cells
    raster_band.SetNoDataValue(no_data_value)
    
    # run vector to raster on new raster with input Shapefile
    gdal.RasterizeLayer(new_raster, [1], shp_layer, burn_values=[255])
    
    new_raster.FlushCache()
    
    del new_raster
    
    
    import numpy as np 
    from osgeo import ogr,osr 
    import matplotlib 
    import matplotlib.path as mpath 
    import matplotlib.patches as mpatches 
    import matplotlib.pyplot as plt 
    from matplotlib.patches import Polygon 
    from matplotlib.collections import PatchCollection
    import json
    import geojson
    from shapely.geometry import asShape
    from shapely import wkt
    from shapely.geometry import GeometryCollection,Point,LineString
    import ogr
    
    
    with open("F:\\builddata\\mkt_wuh\\hzw\\12345\\12_AB_DATA\\1209_AB_SEM_DES\\AD_Arrow.geojson") as f:
        data = f.read()
    
    jsonStr = geojson.loads(data)
    features = jsonStr['features']
    
    
    patches = []
    polygon = asShape(features[0]['geometry'])
    coord_array  = np.array(polygon.exterior)
    coords = []
    for coord in coord_array:
        tmp = coord[:2]
        coords.append(tmp)
    
    data = np.array(coords)
    max_val = np.max(data,axis=0)
    min_val = np.min(data,axis=0)
    
    patches.append(Polygon(np.array(coords)))
    
    fig = plt.figure(figsize=(20,10), dpi=250) 
    ax = fig.add_subplot(111)
    
    collection = PatchCollection(patches)  
    
    ax.add_collection(collection)  
    
    ax.set_xlim(min_val[0],max_val[0]) 
    ax.set_ylim(min_val[1],max_val[1])  
    plt.show()
    

    相关文章

      网友评论

          本文标题:RasterizeLayer

          本文链接:https://www.haomeiwen.com/subject/znpgzktx.html