ea-lidar-uncertainty-review/clean.py at main · rmarowitz/ea-lidar-uncertainty-review · GitHub

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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
import os
import pathlib

import earthpy as et
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import rasterstats as rs
import xarray as xr
import seaborn as sns

class NEONDataLoader(object):
    """Parent class to load NEON height data."""

    base_dir_tmpl = os.path.join(
        'california',
        'neon-{site_name_low}-site')
    insitu_path_tmpl = os.path.join(
        '{base_dir}', '2013',
        'insitu',
        'veg{seperator}structure',
        'D17_2013_{site_name_up}_vegStr.csv')
    plots_path_tmpl = os.path.join(
        '{base_dir}',
        'vector_data',
        '{site_name_up}{plot}_centroids.shp')
    chm_path_tmpl = os.path.join(
        '{base_dir}', '2013',
        'lidar',
        '{site_name_low}_lidarCHM.tif')

    site_name = NotImplemented
    id_col_name = NotImplemented
    formatting_dict = NotImplemented
    id_modifier = None

    def __init__(self):
        self.formatting_dict['site_name_low'] = self.site_name.lower()
        self.formatting_dict['site_name_up'] = self.site_name.upper()
        self.formatting_dict['base_dir'] = (
            self.base_dir_tmpl.format(**self.formatting_dict))

        self.insitu_path = self.insitu_path_tmpl.format(**self.formatting_dict)
        self.chm_path = self.chm_path_tmpl.format(**self.formatting_dict)
        self.plots_path = self.plots_path_tmpl.format(**self.formatting_dict)

        self._insitu_height_stats = None
        self._lidar_chm_stats = None
        self._height_stats = None


    @property
    def lidar_chm_stats(self):
        """
        Calculate max and mean tree height from LiDAR.
        """
        if self._lidar_chm_stats is None:
            plots_gdf = gpd.read_file(self.plots_path)
            plots_gdf.geometry = plots_gdf.geometry.buffer(20)

            # Calculate the zonal stats
            chm_stats = rs.zonal_stats(
                plots_gdf,
                self.chm_path,
                stats=['mean', 'max'],
                geojson_out=True,
                nodata=0,
                copy_properties=True
            )
            self._lidar_chm_stats = gpd.GeoDataFrame.from_features(chm_stats)
            lidar_mean_max = self._lidar_chm_stats.rename(
                columns={'max':'lidar_max', 'mean':'lidar_mean'},
                inplace=True)
            if not self.id_modifier is None:
                self._lidar_chm_stats[self.id_col_name] = (
                    self._lidar_chm_stats[self.id_col_name].apply(self.id_modifier))
        return self._lidar_chm_stats

    @property
    def insitu_height_stats(self):
        """
        Calculate insitu tree height data max and mean.
        """
        if self._insitu_height_stats is None:
            self._insitu_height_stats = (pd.read_csv(self.insitu_path)
                          .groupby('plotid')
                          .stemheight
                          .agg(['max','mean'])
                         .rename(columns={
                                 'max':'insitu_max', 'mean':'insitu_mean'}))
        return self._insitu_height_stats

    @property
    def height_stats(self):
        """
        Calculate lidar and insitu height stats on a merged dataframe.
        """
        if self._height_stats is None:
            self._height_stats = (
                self.lidar_chm_stats
                .merge(self.insitu_height_stats,
                       right_index=True,
                       left_on=self.id_col_name))

        return self._height_stats

    def plots(self):
        """
        Create plots with regression lines.
        """

        fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(14, 7))
        # ax1.scatter(sjer_df.lidar_mean, sjer_df.insitu_mean)
        sns.regplot(data=self.height_stats, x='lidar_max',y='insitu_max',
                   color='purple',
                   ax=ax1)
        ax1.plot((0,1), (0,1), transform=ax1.transAxes, ls='--', c='k')
        ax1.set(xlabel='Max Lidar Height(m)',ylabel='Max Insitu Height(m)',
               title=('Plot 1: Lidar Max vs Insitu Max Tree Height\n'
                     '{} Field Site').format(self.site_name),
               ylim=(0,30),
               xlim=(0,30))

        sns.regplot(data=self.height_stats, x='lidar_mean',y='insitu_mean',
                   color='purple',
                   ax=ax2)
        ax2.set(xlabel='Mean Lidar Height(m)',ylabel='Mean Insitu Height(m)',
               title=('Plot 2: Lidar Mean vs Insitu Mean Tree Height\n'
                     '{} Field Site').format(self.site_name),
                ylim=(0,30),
                xlim=(0,30))
        ax2.plot((0,1), (0,1), transform=ax2.transAxes, ls='--', c='k')