美文网首页
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