学习笔记 | Python处理矢栅数据（2）

1. 导入并查看矢量数据

导入相关库（geopandas,rioxarray,pandas,matplotlib）

import geopandas as gpd
import rioxarray
import pandas as pd
import matplotlib.pyplot as plt

aoi_boundary_HARV = gpd.read_file("G:/learnpy/data/NEONGeopydata/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")

查看矢量数据类型

aoi_boundary_HARV.type
# 0    Polygon
# dtype: object

查看坐标系

aoi_boundary_HARV.crs

# <Projected CRS: EPSG:32618>
# Name: WGS 84 / UTM zone 18N
# Axis Info [cartesian]:
# - E[east]: Easting (metre)
# - N[north]: Northing (metre)
# Area of Use:
# - name: Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.
# - bounds: (-78.0, 0.0, -72.0, 84.0)
# Coordinate Operation:
# - name: UTM zone 18N
# - method: Transverse Mercator
# Datum: World Geodetic System 1984 ensemble
# - Ellipsoid: WGS 84
# - Prime Meridian: Greenwich

查看边界

aoi_boundary_HARV.bounds
# minx	miny	maxx	maxy
# 0	732128.016925	4.713209e+06	732251.102892	4.713359e+06

aoi_boundary_HARV.plot(figsize=(5,5),edgecolor="red",facecolor="green", alpha=0.2)

2. CSV表格文件转矢量数据

CHM_HARV = rioxarray.open_rasterio("G:/learnpy/data/NEONGeopydata/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")
CHM_HARV

# <xarray.DataArray (band: 1, y: 1367, x: 1697)>
# [2319799 values with dtype=float64]
# Coordinates:
#   * band         (band) int32 1
#   * x            (x) float64 7.315e+05 7.315e+05 ... 7.331e+05 7.331e+05
#   * y            (y) float64 4.714e+06 4.714e+06 ... 4.712e+06 4.712e+06
#     spatial_ref  int32 0
# Attributes:
#     STATISTICS_MAXIMUM:  38.169998168945
#     STATISTICS_MEAN:     18.097804669578
#     STATISTICS_MINIMUM:  0
#     STATISTICS_STDDEV:   5.3218337145236
#     _FillValue:          -9999.0
#     scale_factor:        1.0
#     add_offset:          0.0

plot_locations_HARV = gpd.read_file("G:/learnpy/data/NEONGeopydata/NEON-DS-Site-Layout-Files/HARV/HARV_PlotLocations.csv")
plot_locations_HARV

plot_locations_HARV1 = gpd.GeoDataFrame(plot_locations_HARV, 
                                    geometry = gpd.points_from_xy(plot_locations_HARV.easting,plot_locations_HARV.northing), 
                                    crs=CHM_HARV.rio.crs)
plot_locations_HARV1 = gpd.GeoDataFrame(plot_locations_HARV, 
                                    geometry = gpd.points_from_xy(plot_locations_HARV.easting,plot_locations_HARV.northing), 
                                    crs=CHM_HARV.rio.crs)
plot_locations_HARV1.plot(figsize=(6,6),edgecolor="red",facecolor="None")
plt.title("Map of Plot Location")
plt.ticklabel_format(style="plain")

3. 矢量叠加栅格绘图

# Shapefiles 
tower = gpd.read_file("G:/learnpy/data/NEONGeopydata/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
road = gpd.read_file("G:/learnpy/data/NEONGeopydata/NEON-DS-Site-Layout-Files/HARV/HARV_roads.shp")
aoi_boundary_HARV = gpd.read_file("G:/learnpy/data/NEONGeopydata/NEON-DS-Site-Layout-Files/HARV/HarClip_UTMZ18.shp")
CHM_HARV = rioxarray.open_rasterio("G:/learnpy/data/NEONGeopydata/NEON-DS-Airborne-Remote-Sensing/HARV/CHM/HARV_chmCrop.tif")

fg,ax = plt.subplots()
CHM_HARV1 = CHM_HARV.where(CHM_HARV!=CHM_HARV.rio.nodata)
CHM_HARV1.plot(ax=ax)
plot_locations_HARV1.plot(ax = ax, edgecolor="red",facecolor="None")
aoi_boundary_HARV1.plot(ax=ax, edgecolor="blue",facecolor="None")
tower.plot(ax=ax, facecolor="red")
road.plot(ax=ax,edgecolor = "orange", facecolor="None")

4.基于矢量裁剪栅格数据并导出为tif

CHM_HARV_clip = CHM_HARV.rio.clip(aoi_boundary_HARV.geometry)
plt.style.use("ggplot")
plt.figure()
CHM_HARV_clip.plot(figsize = (6,6))
plt.title("Crope the Raster")
plt.ticklabel_format(style="plain")
CHM_HARV_clip.rio.to_raster("G:/learnpy/data/CHM_HARV_clip.tif")