programming python Rasterio

Hướng dẫn raster data analysis in python - phân tích dữ liệu raster trong python

Wur Geoscripting

Tuần 3, Bài 10: Xử lý dữ liệu raster với Python

Jan Verbesselt, Jorge Mendes de Jesus, Aldo Bergsma, Dainius Masiliūnas, David Swinkels, Judith Verstegen, Corné Vreugdenhil-2022-01-17

Giới thiệu

Hôm nay chúng tôi sẽ làm việc với các gói Python để phân tích raster không gian. Python có một số gói chuyên dụng để xử lý Rasters:

Owslib để tải xuống dữ liệu raster không gian địa lý từ các dịch vụ bảo hiểm web
GDAL là thư viện mạnh mẽ để đọc, viết và làm cong bộ dữ liệu raster
Rasterio đọc và viết dữ liệu raster không gian địa lý
Numpy là gói cơ bản để tính toán khoa học, chẳng hạn như tính toán mảng [do đó raster]
RasterStats tóm tắt các bộ dữ liệu raster không gian địa lý dựa trên Vector Geometires

Ngày nay, các mục tiêu học tập

Có thể đọc các định dạng raster không gian từ các dịch vụ và tệp web
Có thể viết các định dạng raster không gian vào đĩa
Biết cách áp dụng các hoạt động cơ bản trên dữ liệu raster, chẳng hạn như arithmetics
Có thể vẽ dữ liệu raster không gian với matplotlib

Thiết lập môi trường Python

Tạo cấu trúc thư mục cho hướng dẫn này:

cd ~/Documents/
mkdir PythonRaster #or give the directory a name to your liking
cd ./PythonRaster
mkdir data
mkdir output

Giống như trong bài học trước, chúng tôi sẽ tạo một môi trường Conda với một tệp

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

8 và sau đó cập nhật nó với tệp

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

8 thứ hai. Sử dụng nội dung sau cho hai tệp:

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

name: raster
channels:
  - conda-forge
dependencies:
  - rasterio
  - rasterstats
  - affine

Sử dụng hai tệp này để tạo một môi trường cho hướng dẫn này. Lưu ý: cái thứ hai có thể mất một thời gian [5-10 phút]; Chỉ cần đi và lấy cà phê/trà/nước/nước trái cây. Cuối cùng, kích hoạt môi trường, mở Spyder, tạo một tập lệnh trong thư mục gốc và bắt đầu mã hóa.

Thông qua dịch vụ bảo hiểm web

Dịch vụ bảo hiểm web [WCS] tải dữ liệu raster theo cách tương tự như dữ liệu vector tải dịch vụ tính năng web [WFS]. Các dịch vụ bảo hiểm web là một tiêu chuẩn của Hiệp hội không gian địa lý mở và cho phép tải xuống dữ liệu raster không gian địa lý với nhiều loại mã hóa định dạng: GeoTiff, NetCDF, JPEG2000, v.v. Nó cho phép tải xuống hình ảnh nhưng không có giá trị dữ liệu.

Hôm nay chúng tôi sẽ làm việc với các Rasters độ cao. Hãy xem WCS của bộ dữ liệu AHN. AHN là viết tắt của nhóm Actueel Hoogtebestand Nederland và là một mô hình độ cao kỹ thuật số [DEM] bao gồm Hà Lan. Truy cập dịch vụ bảo hiểm web để xem nội dung.

from owslib.wcs import WebCoverageService

# Access the WCS by proving the url and optional arguments [here version]
wcs = WebCoverageService['//geodata.nationaalgeoregister.nl/ahn2/wcs?service=WCS', version='1.0.0']
# Print to check the contents of the WCS
print[list[wcs.contents]]

Chạy dòng mã cuối cùng cho thấy dịch vụ bảo hiểm web của AHN2 chứa bốn raster: 0,5m được nội suy, 0,5M không được nội suy, 0,5m thô và 5M. Các mét xác định kích thước tế bào. Dữ liệu raster của AHN có hệ tọa độ dự kiến RD_NEW [EPSG: 28992].

Chúng tôi cũng có thể kiểm tra loại hoạt động nào có sẵn trong WCS:

# Get all operations and print the name of each of them
print[[op.name for op in wcs.operations]]

Bạn sẽ thấy rằng dịch vụ bảo hiểm web cho phép truy cập dữ liệu [getCoverage], siêu dữ liệu [mô tả] và các khả năng [getCapabilities]. Đây là tất cả các giao thức tiêu chuẩn được xác định bởi OGC.

Một số chức năng có sẵn để truy cập siêu dữ liệu cụ thể của từng raster riêng lẻ, ví dụ:

# Take the 0.5m-resolution rough raster as example
cvg = wcs.contents['ahn2_05m_ruw']
# print supported reference systems, the bounding box defined in WGS 84 coordinates, and supported file formats
print[cvg.supportedCRS]
print[cvg.boundingBoxWGS84]
print[cvg.supportedFormats]

Câu 1: Hệ thống tham chiếu tọa độ của Raster thô 0,5m là gì? Có giống nhau cho ba raster khác không?: What is the coordinate reference system of the 0.5m rough raster? Is it the same for the other three rasters?

Chúng ta hãy xem dữ liệu. Chúng tôi không muốn quá tải dịch vụ web. Do đó, chúng tôi tải xuống một lần và lưu trữ dữ liệu cục bộ.

Tải xuống mô hình bề mặt kỹ thuật số [DSM], là phiên bản ‘AHN2_05M_RUW và Mô hình địa hình kỹ thuật số [DTM], là phiên bản‘ AHN2_05M_INT, vào một tệp cục bộ. Sự khác biệt giữa một DEM, DSM và DTM được giải thích trên GIS Stackexchange.

import os

# Define a bounding box in the availble crs [see before] by picking a point and drawing a 1x1 km box around it
x, y = 174100, 444100
bbox = [x-500, y-500, x+500, y+500]

# Create a data directoty within the directory where this script is run if it does not exist yet and store file
if not os.path.exists['data']: os.makedirs['data']

# Request the DSM data from the WCS
response = wcs.getCoverage[identifier='ahn2_05m_ruw', bbox=bbox, format='GEOTIFF_FLOAT32', crs='urn:ogc:def:crs:EPSG::28992', resx=0.5, resy=0.5]

# Write the data to a local file in the 'data' folder [it should exist]
with open['./data/AHN2_05m_DSM.tif', 'wb'] as file:
    file.write[response.read[]]

# Do the same for the DTM
response = wcs.getCoverage[identifier='ahn2_05m_int', bbox=bbox, format='GEOTIFF_FLOAT32', crs='urn:ogc:def:crs:EPSG::28992', resx=0.5, resy=0.5]

with open['./data/AHN2_05m_DTM.tif', 'wb'] as file:
    file.write[response.read[]]

Từ một tập tin với gdal

GDAL xử lý các định dạng dữ liệu không gian địa lý raster và vector với Python, Java, R và C API. Khi mở tệp raster trong GDAL, đối tượng có cấu trúc phân cấp bắt đầu ở cấp dữ liệu. Một bộ dữ liệu có một địa chất chuyển đổi [siêu dữ liệu] và có thể chứa một hoặc nhiều dải. Mỗi ban nhạc có một mảng dữ liệu và tổng quan có khả năng.

Cấu trúc lớp GDAL, được chuyển thể từ Garrard, 2016, xử lý địa lý với Python.

Hãy để chúng tôi mở tập tin chúng tôi vừa lưu. Trước tiên, bạn sẽ thấy bạn lấy tập dữ liệu và cần truy cập băng tần [mặc dù chỉ có một], trước khi có thể truy cập mảng dữ liệu.

from osgeo import gdal

# Open dataset, gdal automatically selects the correct driver
ds = gdal.Open["./data/AHN2_05m_DTM.tif" ]

# Get the band [band number 1]
band = ds.GetRasterBand[1]
        
# Get the data array
data = band.ReadAsArray[]
print[data]

# Delete objects to close the file
ds = None

Câu 2: Tại sao chúng tôi đặt DS thành không có ở cuối tập lệnh của bạn? Điều gì có thể xảy ra nếu bạn không làm điều đó?: Why do we set ds to None at the end of your script? What may happen if you do not do that?

API GDAL Python không phải là mô -đun Python được ghi chép tốt nhất. Do đó, Rasterio được giải thích như một mô -đun xử lý dữ liệu raster thay thế.

Từ một tệp với rasterio

Rasterio đọc và ghi nhiều định dạng raster dựa trên GDAL, cung cấp các chức năng xử lý raster dựa trên các mảng numpy và Geojson, và tích hợp matplotlib trong mô -đun rasterio.plot để trực quan hóa. Rasterio xử lý nhiều ban nhạc, mặt nạ, phản hồi và lấy mẫu lại.

Phần còn lại của bài học dưới đây là một tuyến đường hoàn chỉnh để xử lý rasterdataset. Chúng tôi sẽ sử dụng DEMS từ WCS cho khu vực nghiên cứu của chúng tôi, xử lý nó với rasterio, tính toán thông tin mới [CHM], phủ nó với dữ liệu vector đại diện cho các tòa nhà và trực quan hóa nó.

Chúng ta hãy đọc trong dữ liệu raster mà chúng ta vừa lưu trữ từ WCS với Rasterio và vẽ nó bằng rasterio.plot:

import rasterio
from rasterio.plot import show

# open the two rasters 
dsm = rasterio.open["./data/AHN2_05m_DSM.tif", driver="GTiff"]
dtm = rasterio.open["./data/AHN2_05m_DTM.tif", driver="GTiff"]

# metadata functions from rasterio
print[dsm.meta]
print[dtm.meta]

# plot with rasterio.plot, which provides matplotlib functionality
show[dsm, title='Digital Surface Model', cmap='gist_ncar']

Câu 3: Có phải bản đồ của Colormap ‘Gist_NCAR, là một lựa chọn tốt cho mô hình bề mặt kỹ thuật số? Tại sao tại sao không? Nếu không, hãy cố gắng tìm một cái phù hợp hơn.: Is the colormap ‘gist_ncar’ cartographically a good choice for a digital surface model? Why/Why not? If not, try to find a more fitting one.

Siêu dữ liệu hiển thị trình điều khiển [cách biết của GDAL cách hoạt động với định dạng tệp cụ thể], kiểu dữ liệu, giá trị NODATA, chiều rộng của raster về số lượng ô, chiều cao của raster về số lượng ô, số lượng dải raster trong bộ dữ liệu, phối hợp Hệ thống tham chiếu và giá trị chuyển đổi.

Ở phần cuối, các lớp raster trong rasterio được lưu trữ dưới dạng mảng numpy, xuất hiện khi dữ liệu được đọc bằng phương pháp

name: raster
channels:
  - conda-forge
dependencies:
  - rasterio
  - rasterstats
  - affine


# rasterio object
print[type[dsm]]

# read, show object type and data
dsm_data = dsm.read[1]
print[type[dsm_data]]
print[dsm_data]

Điện toán mô hình chiều cao tán

Một mô hình chiều cao tán [CHM] cho thấy một dấu hiệu về chiều cao của cây. Nó có thể được tạo ra bằng cách trừ một mô hình địa hình kỹ thuật số từ mô hình bề mặt kỹ thuật số. Trong raster kết quả, mỗi giá trị ô biểu thị độ cao của cây trên địa hình bề mặt bên dưới. Rasterio dựa vào Numpy để thực hiện các hoạt động Raster Point.

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Câu 4: Đâu là tán cây cao nhất trong khu vực nghiên cứu? Đó có phải là những gì bạn mong đợi?: Where is the canopy the highest in the study area? Is that what you expected?

Mặc dù chúng tôi đã áp dụng các khái niệm về mô hình chiều cao tán, nhưng khu vực chúng tôi chọn không chỉ chứa thảm thực vật trên đỉnh bề mặt mà còn cả các tòa nhà. Chiều cao của thảm thực vật và các tòa nhà do đó hiện đang được kết hợp trong raster CHM của chúng tôi.

Điện toán chiều cao của các tòa nhà

Hãy để chúng tôi xác định chiều cao của các tòa nhà. Bước đầu tiên là tải xuống dữ liệu xây dựng từ dịch vụ tính năng Web Bag mà chúng tôi đã sử dụng trong hướng dẫn Vector.

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Câu 5: Điều gì xảy ra trong dòng BB BB = ‘,, Tham gia [MAP [STR, BBOX]]? Tra cứu cách lập bản đồ hoạt động nếu bạn không biết.: What happens in the line “bb = ‘,’.join[map[str, bbox]]”? Look up how mapping works if you do not know.

Bước tiếp theo là thực hiện số liệu thống kê khu vực để có được giá trị chiều cao trung bình trên mỗi đa giác xây dựng. Chúng tôi sẽ làm điều này với các rasterstats mô -đun, sử dụng các tệp geodataFrames và .tif cho tác vụ này. Nó xuất hiện một Geojson.

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Bạn có thể nhận được một tin nhắn cảnh báo từ Rasterstats; Bạn có thể bỏ qua điều này.

Trực quan hóa nhanh cho chúng ta thấy độ cao xuất phát từ dữ liệu raster trên bản đồ:

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Câu 6: Tại sao chúng ta muốn có một thang đo bằng nhau theo hướng x và y cho con số này?: Why do we want an equal scale in the x and y direction for this figure?

Viết dữ liệu raster vào một tệp

Để lưu trữ mảng Numpy dưới dạng tệp raster, Rasterio cần siêu dữ liệu đi kèm, như bạn đã thấy trước đây. Có thể sử dụng siêu dữ liệu của một raster hiện có hoặc để tạo nó từ đầu.

Để tạo siêu dữ liệu từ đầu, CRS có thể được xác định với hàm từ rasterio và chuyển đổi với affine. Affine là một mô -đun Python tạo điều kiện cho các phép biến đổi affine, tức là & nbsp; tỷ lệ, xoay, phản chiếu hoặc sai lệch hình ảnh/raster/mảng.

Rasterio có thể viết hầu hết các định dạng raster từ GDAL. Các nhà phát triển khuyên bạn nên sử dụng trình điều khiển địa lý để viết vì đây là định dạng được thử nghiệm tốt nhất và được hỗ trợ tốt nhất.

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Trực quan hóa dữ liệu raster

Dữ liệu raster có thể được hiển thị bằng cách chuyển các mảng numpy cho matplotlib trực tiếp hoặc thứ hai thông qua một phương thức trong rasterio truy cập matplotlib cho bạn. Sử dụng matplotlib trực tiếp cho phép linh hoạt hơn, chẳng hạn như điều chỉnh huyền thoại, trục và nhãn, và phù hợp hơn cho các mục đích chuyên nghiệp. Việc trực quan hóa thông qua Rasterio yêu cầu ít mã hơn và có thể đưa ra ý tưởng nhanh về dữ liệu raster của bạn. Chúng tôi sẽ hiển thị cả hai cách tiếp cận. Hãy để chúng tôi thực hiện trực quan hóa DSM và DTM với matplotlib:

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Nếu bạn không thích colormap màu cam của matplotlib, có thể chọn một colormap khác.

Cách tiếp cận thứ hai với Rasterio chỉ yêu cầu một dòng mã để tạo ra một cốt truyện. Bằng cách tạo các ô con, các số liệu có thể được kết hợp [cũng có thể được thực hiện với matplotlib trực tiếp].

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Rasterio cũng trực quan hóa biểu đồ đơn giản bằng cách gọi các chức năng của matplotlib.

name: raster
dependencies:
  - python
  - numpy
  - matplotlib
  - spyder
  - owslib
  - gdal
  - geopandas

Câu 7: Điều gì được biểu diễn trên trục X và Y? Các nhãn trục mặc định là DN [x] và tần số [y]; Nếu bạn thay đổi chúng, bạn sẽ chọn nhãn nào để phản ánh tốt hơn nội dung của các lô?: What is represented on the x and y axis? The default axis labels are DN [x] and Frequency [y]; if you were to change them, what labels would you pick to better reflect the content of the plots?

Hướng dẫn làm việc với Rasters ở Python với rasterio
Hướng dẫn làm việc với Raster trong Python với GDAL [cho Python 2]
Hình ảnh vệ tinh Landsat
Hình ảnh vệ tinh Landsat được lấy lại
Hình ảnh vệ tinh Sentinel

Dữ liệu raster trong Python là gì?

Dữ liệu raster hoặc lưới lưới được lưu trữ dưới dạng một lưới các giá trị được hiển thị trên bản đồ dưới dạng pixel. Mỗi giá trị pixel đại diện cho một khu vực trên bề mặt trái đất. Một tệp raster bao gồm lưới các ô thông thường, tất cả đều có cùng kích thước.a grid of values which are rendered on a map as pixels. Each pixel value represents an area on the Earth's surface. A raster file is composed of regular grid of cells, all of which are the same size.

Phân tích dữ liệu raster là gì?

Phân tích dữ liệu raster dựa trên các ô và rasters.● Phân tích dữ liệu raster có thể được thực hiện ở cấp độ của từng tế bào hoặc các nhóm ô hoặc tế bào trong toàn bộ raster.● Một số hoạt động dữ liệu raster sử dụng một raster duy nhất;Những người khác sử dụng hai hoặc nhiều rasters.based on cells and rasters. ● Raster data analysis can be performed at the level of individual cells, or groups of cells, or cells within an entire raster. ● Some raster data operations use a single raster; others use two or more rasters.

Làm thế nào để Python xử lý dữ liệu raster?

Xử lý dữ liệu raster với Python và GDAL..

Trích xuất tệp.Tải xuống dữ liệu cho bài học này từ liên kết này.....

Mở một tập tin raster trong Python với GDAL.Với GDAL, bạn có thể đọc và viết một số định dạng raster khác nhau trong Python.....

Mở tệp raster dưới dạng một mảng số.....

Tiện ích dòng lệnh GDAL ..

Ba loại dữ liệu raster là gì?

Có ba loại dữ liệu raster có thể được lưu trữ trong cơ sở dữ liệu địa lý: bộ dữ liệu raster, danh mục raster và raster là thuộc tính.Các bộ dữ liệu raster là hình ảnh đơn được lưu trữ trong cơ sở dữ liệu.raster datasets, raster catalogs, and raster as attributes. Raster datasets are single images that are stored in the database.