-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGEDI_L2A_forAge_pointQuery.py
66 lines (63 loc) · 2.91 KB
/
GEDI_L2A_forAge_pointQuery.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
# Author: Milutin Milenkovic
# Copyright:
# Licence:
# ----------------------------------------------------------------------------------------------------------------------
# -- Short Description --
"""
This script assigns the forest age class and border pixel flag to GEDI shots.
"""
# -------- Input --------
# (1) GEDI shots over secondary forest in the Rondônia state, Brazil.
# (2) a forest age raster: Silva Junior et al. 2020, https://doi.org/10.6084/m9.figshare.12622025
# (3) a black-and-withe version of the forest age raster with border pixels removed
# -------- Output -------
# (1) The input GEDI shots with forest age attribute and border pixel flag
# ----------------------------------------------------------------------------------------------------------------------
import geopandas as gpd
import pandas as pd
import numpy as np
from rasterstats import point_query
# ------------------
# Forest Age rasters:
# ------------------
forestAge_file = r'/mnt/ssd/milutin/Para_upScaling/StandAgeMapBiomas/svbr-rondonia-2018.tif'
#
out_eroded_forestAge_file = r'/mnt/ssd/milutin/Para_upScaling/StandAgeMapBiomas/svbr-rondonia-2018_bw_eroded.tif'
# -----------------------------
# gedi csv file with all data
# -----------------------------
gediL2A_file = r'/mnt/raid/milutin/upScaling/Rondonia/GEDI/Rondonia_L2A_v002/output/gedi_L2A_Rondonia_all_sens_a2.csv'
# ##############################################
# data reading and processing
# ##############################################
# read the GEDI points
gedi_df = pd.read_csv(gediL2A_file)
# ------------------------
# convert to geopandas df
# ------------------------
gedi_gdf = gpd.GeoDataFrame(gedi_df, geometry=gpd.points_from_xy(gedi_df.Longitude, gedi_df.Latitude))
gedi_gdf.set_crs(epsg=4326, inplace=True)
# ------------------------------------
# calculate zonal statistics
# ------------------------------------
forestAge_list = point_query(vectors=gedi_gdf['geometry'], raster=forestAge_file, interpolate='nearest')
gedi_gdf['forestAge'] = np.asarray(forestAge_list)
# ------------------------------------
# remove points with stand age == 0
# ------------------------------------
gedi_gdf = gedi_gdf[gedi_gdf['forestAge'] != 0]
# ------------------------------------
# remove points with stand age == None
# ------------------------------------
gedi_gdf.dropna(subset=['forestAge'], inplace=True)
# ----------------------------------------------
# get pixel value from the eroded bw forest age
# ----------------------------------------------
forestAge_list = point_query(vectors=gedi_gdf['geometry'], raster=out_eroded_forestAge_file, interpolate='nearest')
gedi_gdf['eroded_age'] = np.asarray(forestAge_list)
# ------------------------------------
# save:
# ------------------------------------
outFilePath = r'/mnt/raid/milutin/upScaling/Rondonia/GEDI/Rondonia_L2A_v002/output/directAnalysis/gedi_L2A_gdf_sens_a2.json'
with open(outFilePath, 'w') as f:
f.write(gedi_gdf.to_json())