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")