Hướng dẫn raster to shapefile python - raster để shapefile trăn

Mục tiêu học tập

• Cắt một bộ dữ liệu raster bằng Python bằng cách sử dụng một đối tượng phạm vi vector có nguồn gốc từ một Shapefile.Python using a vector extent object derived from a shapefile.
• Mở một Shapefile trong Python.Python.

Trong các bài học trước, bạn đã phân loại lại một raster trong Python; Tuy nhiên, các cạnh của bộ dữ liệu raster của bạn không đồng đều.Python; however, the edges of your raster dataset were uneven.

Trong bài học này, bạn sẽ học cách cắt ra raster - để tạo một đối tượng / tệp raster mới mà bạn có thể chia sẻ với các đồng nghiệp và / hoặc mở trong các công cụ khác như công cụ GIS trên máy tính để bàn như QGIS.

Về cây trồng không gian

Cắt xén (đôi khi còn được gọi là cắt), là khi bạn tập hợp con hoặc làm cho một bộ dữ liệu nhỏ hơn, bằng cách loại bỏ tất cả dữ liệu bên ngoài khu vực cây trồng hoặc phạm vi không gian. Trong trường hợp này, bạn có một raster lớn - nhưng hãy để giả vờ rằng bạn chỉ cần làm việc với một tập hợp con nhỏ hơn của raster.

Bạn có thể sử dụng chức năng

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

1. Làm cho dữ liệu nhỏ hơn và
2. Làm cho việc xử lý và vẽ đồ thị nhanh hơn

Nói chung khi bạn có thể, nó thường là một ý tưởng tốt để cắt dữ liệu raster của bạn!

Để bắt đầu, hãy để tải các thư viện mà bạn sẽ cần trong bài học này.

Tải thư viện

import os
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from shapely.geometry import mapping
import rioxarray as rxr
import xarray as xr
import geopandas as gpd

import earthpy as et
import earthpy.plot as ep

# Prettier plotting with seaborn
sns.set(font_scale=1.5)

# Get data and set working directory
et.data.get_data("colorado-flood")
os.chdir(os.path.join(et.io.HOME,
                      'earth-analytics',
                      'data'))

Downloading from https://ndownloader.figshare.com/files/16371473
Extracted output to /root/earth-analytics/data/colorado-flood/.

Mở các lớp raster và vector

Trong các bài học trước, bạn đã làm việc với một lớp raster trông giống như lớp bên dưới. Lưu ý rằng dữ liệu có cạnh không đều ở phía bên trái. Hãy để giả vờ cạnh này nằm ngoài khu vực nghiên cứu của bạn và bạn muốn loại bỏ nó hoặc cắt nó bằng cách sử dụng phạm vi khu vực nghiên cứu của bạn. Bạn có thể làm điều này bằng cách sử dụng hàm

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

lidar_chm_path = os.path.join("colorado-flood", 
                              "spatial"
                              "boulder-leehill-rd",
                              "outputs",
                              "lidar_chm.tif")

lidar_chm_im = rxr.open_rasterio("colorado-flood/spatial/boulder-leehill-rd/outputs/lidar_chm.tif",
                                 masked=True).squeeze()

f, ax = plt.subplots(figsize=(10, 5))
lidar_chm_im.plot.imshow()
ax.set(title="Lidar Canopy Height Model (CHM)")

ax.set_axis_off()
plt.show()

Lớp vector mở

Để bắt đầu clip của bạn, hãy mở một lớp vector chứa mức độ cây trồng mà bạn muốn sử dụng để cắt dữ liệu của mình. Để mở ShapeFile, bạn sử dụng chức năng

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

aoi = os.path.join("colorado-flood",
                   "spatial",
                   "boulder-leehill-rd",
                   "clip-extent.shp")

# Open crop extent (your study area extent boundary)
crop_extent = gpd.read_file(aoi)

Tiếp theo, xem Hệ thống tham chiếu tọa độ (CRS) của cả hai bộ dữ liệu của bạn. Hãy nhớ rằng để thực hiện bất kỳ phân tích nào với hai bộ dữ liệu này cùng nhau, chúng sẽ cần phải ở trong cùng một CRS.

print('crop extent crs: ', crop_extent.crs)
print('lidar crs: ', lidar_chm_im.rio.crs)

crop extent crs:  epsg:32613
lidar crs:  EPSG:32613

# Plot the crop boundary layer
# Note this is just an example so you can see what it looks like
# You don't need to plot this layer in your homework!
fig, ax = plt.subplots(figsize=(6, 6))

crop_extent.plot(ax=ax)

ax.set_title("Shapefile Crop Extent",
             fontsize=16)
plt.show()

Bây giờ bạn đã nhập Shapefile. Bạn có thể sử dụng chức năng

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

Clip dữ liệu raster bằng Rioxarray import rasterio import rasterio.features from shapely.geometry import shape ds = rasterio.open('20211124_092953_Data-001.jpg') ds2 = rasterio.open('20211124_092953_Data-001.vrt') band = ds2.read(1) # make a copy of original single band numpy array raster_array = band.copy() # replace multiple values with 1 to speed up the process of polgenizing raster_array = np.where(raster_array>0,1,raster_array) # numpy array binary mask loaded however needed maskShape = rasterio.features.shapes(raster_array.astype('uint8')) mypoly=[] for vec in maskShape: mypoly.append(vec[0]) 0

Nếu bạn muốn cắt dữ liệu, bạn có thể sử dụng chức năng

import rasterio
import rasterio.features
from shapely.geometry import shape

ds = rasterio.open('20211124_092953_Data-001.jpg')
ds2 = rasterio.open('20211124_092953_Data-001.vrt')
band = ds2.read(1)
# make a copy of original single band numpy array
raster_array = band.copy()
# replace multiple values with 1 to speed up the process of polgenizing
raster_array = np.where(raster_array>0,1,raster_array)


#  numpy array binary mask loaded however needed
maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
mypoly=[]
for vec in maskShape:
    mypoly.append(vec[0])

Để thực hiện clip bạn:

1. Mở bộ dữ liệu raster mà bạn muốn cắt bằng Xarray hoặc Rioxarray.
2. Mở Shapefile của bạn như một đối tượng Geopandas.
3. Cắt dữ liệu bằng hàm

   import rasterio
   import rasterio.features
   from shapely.geometry import shape
   
   ds = rasterio.open('20211124_092953_Data-001.jpg')
   ds2 = rasterio.open('20211124_092953_Data-001.vrt')
   band = ds2.read(1)
   # make a copy of original single band numpy array
   raster_array = band.copy()
   # replace multiple values with 1 to speed up the process of polgenizing
   raster_array = np.where(raster_array>0,1,raster_array)
   
   
   #  numpy array binary mask loaded however needed
   maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
   mypoly=[]
   for vec in maskShape:
       mypoly.append(vec[0])
   

   2.

import rasterio
import rasterio.features
from shapely.geometry import shape

ds = rasterio.open('20211124_092953_Data-001.jpg')
ds2 = rasterio.open('20211124_092953_Data-001.vrt')
band = ds2.read(1)
# make a copy of original single band numpy array
raster_array = band.copy()
# replace multiple values with 1 to speed up the process of polgenizing
raster_array = np.where(raster_array>0,1,raster_array)


#  numpy array binary mask loaded however needed
maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
mypoly=[]
for vec in maskShape:
    mypoly.append(vec[0])

• import rasterio
  import rasterio.features
  from shapely.geometry import shape
  
  ds = rasterio.open('20211124_092953_Data-001.jpg')
  ds2 = rasterio.open('20211124_092953_Data-001.vrt')
  band = ds2.read(1)
  # make a copy of original single band numpy array
  raster_array = band.copy()
  # replace multiple values with 1 to speed up the process of polgenizing
  raster_array = np.where(raster_array>0,1,raster_array)
  
  
  #  numpy array binary mask loaded however needed
  maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
  mypoly=[]
  for vec in maskShape:
      mypoly.append(vec[0])
  

  4: Mặc định. Đặt nó sẽ thả tất cả các pixel bên ngoài phạm vi clip

• import rasterio
  import rasterio.features
  from shapely.geometry import shape
  
  ds = rasterio.open('20211124_092953_Data-001.jpg')
  ds2 = rasterio.open('20211124_092953_Data-001.vrt')
  band = ds2.read(1)
  # make a copy of original single band numpy array
  raster_array = band.copy()
  # replace multiple values with 1 to speed up the process of polgenizing
  raster_array = np.where(raster_array>0,1,raster_array)
  
  
  #  numpy array binary mask loaded however needed
  maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
  mypoly=[]
  for vec in maskShape:
      mypoly.append(vec[0])
  

  5: Mặc định. Nếu được đặt thành true, nó sẽ kẹp tất cả dữ liệu bên trong phạm vi clip

• import rasterio
  import rasterio.features
  from shapely.geometry import shape
  
  ds = rasterio.open('20211124_092953_Data-001.jpg')
  ds2 = rasterio.open('20211124_092953_Data-001.vrt')
  band = ds2.read(1)
  # make a copy of original single band numpy array
  raster_array = band.copy()
  # replace multiple values with 1 to speed up the process of polgenizing
  raster_array = np.where(raster_array>0,1,raster_array)
  
  
  #  numpy array binary mask loaded however needed
  maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
  mypoly=[]
  for vec in maskShape:
      mypoly.append(vec[0])
  

  6: Nếu Shapefile của bạn ở trong một CR khác so với dữ liệu raster, hãy truyền CRS để đảm bảo dữ liệu được cắt chính xác.

Dưới đây bạn cắt dữ liệu đến mức của AOI GeodataFrame được nhập ở trên. Dữ liệu sau đó được vẽ.

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

Tùy chọn - xuất khẩu raster mới cắt

Khi bạn đã cắt dữ liệu của mình, bạn có thể muốn xuất nó sang một tệp địa lý mới, giống như bạn đã làm trong các bài học trước.

Bạn cũng có thể sử dụng Rioxarray!

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    with rasterio.open('20211124_092953_Data-001.vrt') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms = list(results)

from shapely.geometry import shape
shape(geoms[0]['geometry'])

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)

Thử thách tùy chọn: Thay đổi cây trồng theo các lớp thời gian

Trong bài học trước, bạn đã tạo 2 lô:

1. Một bản đồ raster được phân loại cho thấy sự thay đổi tích cực và tiêu cực trong mô hình chiều cao tán trước và sau trận lụt. Để làm điều này, bạn sẽ cần tính toán sự khác biệt giữa hai mô hình chiều cao tán.positive and negative change in the canopy height model before and after the flood. To do this you will need to calculate the difference between two canopy height models.
2. Một bản đồ raster được phân loại cho thấy sự thay đổi tích cực và tiêu cực trong địa hình được trích xuất từ ​​các mô hình địa hình kỹ thuật số trước và sau lũ trước và sau trận lụt.positive and negative change in terrain extracted from the pre and post flood Digital Terrain Models before and after the flood.

Tạo hai lô giống nhau ngoại trừ lần này cắt mỗi người trong số các Rasters mà bạn đã vẽ bằng cách sử dụng Shapefile:

import rasterio
import rasterio.features
from shapely.geometry import shape

ds = rasterio.open('20211124_092953_Data-001.jpg')
ds2 = rasterio.open('20211124_092953_Data-001.vrt')
band = ds2.read(1)
# make a copy of original single band numpy array
raster_array = band.copy()
# replace multiple values with 1 to speed up the process of polgenizing
raster_array = np.where(raster_array>0,1,raster_array)


#  numpy array binary mask loaded however needed
maskShape = rasterio.features.shapes(raster_array.astype('uint8'))
mypoly=[]
for vec in maskShape:
    mypoly.append(vec[0])

Đối với mỗi lô, hãy chắc chắn:

• Thêm một huyền thoại hiển thị rõ ràng những gì mỗi màu trong raster phân loại của bạn đại diện.
• Sử dụng màu sắc thích hợp.
• Thêm một tiêu đề vào cốt truyện của bạn.

Bạn sẽ bao gồm những lô này trong báo cáo cuối cùng của bạn vào tuần tới.