#!/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()
网友评论