#如何使用GDAL/OGR打開矢量文件、讀取屬性表字段，并將數據范圍和每個ploygon要素的范圍
import ogr,sys,os
import numpy as np

os.chdir('D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6柵格數據操作/shiyan')

#設置driver,并打開矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open('polygon.shp', 0)
if ds is None:
  print("Could not open", 'sites.shp')
  sys.exit(1)
#獲取圖冊
layer = ds.GetLayer()

#要素數量
numFeatures = layer.GetFeatureCount()
print("Feature count: "+str(numFeatures))

#獲取范圍
extent = layer.GetExtent()
print("Extent:", extent)
print("UL:", extent[0],extent[3])
print("LR:", extent[1],extent[2])

#獲取要素
#feature = layer.GetNextFeature()
 
#循環每個要素屬性
for i in range(numFeatures):
  feature = layer.GetNextFeature(i)
  #獲取字段“id”的屬性
  id = feature.GetField('type')
  #獲取空間屬性
  print(id)
  geometry = feature.GetGeometryRef()
  #x = geometry.GetX()
  polygonextent=geometry.GetEnvelope()
  print(geometry.GetEnvelope())
  #print(y)

  #y = geometry.GetY()
  print("UL:", polygonextent[0],polygonextent[3])
  print("LR:", polygonextent[1],polygonextent[2])
  

# 先獲取polygon.shp外界矩形 的范圍，以及坐標點；  然后將坐標點轉換為 柵格的行列號，只對sample_region范圍內的柵格區域，進行柵格轉point
%matplotlib inline
import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import ogr,sys,os
from pyproj import Transformer
import numpy as np
from osgeo import gdal

os.chdir('D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/shiyan/polygon')


#設置driver,并打開矢量文件
driver = ogr.GetDriverByName('ESRI Shapefile')
dsshp = driver.Open('sample_regions.shp', 0)
if dsshp is None:
  print("Could not open", 'sites.shp')
  sys.exit(1)
#獲取圖冊
layer = dsshp.GetLayer()

#要素數量
numFeatures = layer.GetFeatureCount()
print("Feature count: "+str(numFeatures))

def world2Pixel(geotransform, x, y):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    line = int((y-originY)/pixelHeight)+1
    column = int((x-originX)/pixelWidth)+1
    return (line,column)
def Pixel2world(geotransform, line, column):
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    x = column*pixelWidth + originX - pixelWidth/2
    y = line*pixelHeight + originY - pixelHeight/2
    return(x,y)

ds = gdal.Open( "D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6柵格數據操作/landsat8_20180523/20180523.img" )

geotransform = ds.GetGeoTransform()

 
# 獲取 空間坐標 范圍 
extent = layer.GetExtent()
print("Extent:", extent)
print("UL:", extent[0],extent[3])
print("LR:", extent[1],extent[2])

# 將獲取研究區polygon.shp的空間坐標，轉換成 柵格數據的行列數 
linex1,columny1 = world2Pixel(geotransform, extent[0],extent[3])
print(f"linex1:{linex1}")
print(f"columny1:{columny1}")

linex2,columny2 = world2Pixel(geotransform,extent[1],extent[2])
print(f"linex2:{linex2}")
print(f"columny2:{columny2}")


points = []
for i in range(linex1,linex2):
    for j in range(columny1,columny2):
        x,y = Pixel2world(geotransform,i+1,j+1)
        point = Point((x,y))
        points.append(point)
gdf = gpd.GeoDataFrame()
gdf["geometry"] = points

#將柵格數據的坐標系 賦予 感興趣區的point數據
print(ds.GetProjection())
gdf.crs=ds.GetProjection()
gdf.to_file("D:/20200309/shujushouji/python+jupyternoterbook/zuoye1(1)/6柵格數據操作/test/d_points1.shp") 
gdf.plot()