Andrew Mahon

label_encoder = sklearn.preprocessing.LabelEncoder()
for key, dataframe in fires_df.items():
    fires_df[key].loc[:,'STATE_CODE'] = label_encoder.fit_transform(dataframe['STATE'])

label_encoder = sklearn.preprocessing.LabelEncoder()
for key, dataframe in fires_df.items():
    fires_df[key].loc[:,'STAT_CAUSE_CODE_ZERO'] = label_encoder.fit_transform(dataframe['STAT_CAUSE_CODE'])

 for key, dataframe in fires_df.items():
    dataframe['DISCOVERY_WEEK'] = dataframe['DISCOVERY_DOY']\
        .apply(lambda x: math.floor(x/7) + 1)

climate_regions = {
    "northwest": ["WA","OR","ID"],
    "west": ["CA","NV"],
    "southwest": ["UT","CO","AZ","NM"],
    "northern_rockies": ["MT","ND","SD","WY","NE"],
    "upper_midwest": ["KS","OK","TX","AR","LA","MS"],
    "south": ["MN","WI","MI","IA"],
    "ohio_valley": ["MO","IL","IN","OH","WV","KY","TN"],
    "southeast": ["VA","NC","SC","GA","AL","FL"],
    "northeast": ["ME","NH","VT","NY","PA","MA","RI","CT","NJ","DE","MD","DC"],
    "alaska": ["AK",],
    "hawaii": ["HI"],
    "puerto_rico": ["PR"]
}

state_region_mapping = {}
for region, region_states in climate_regions.items():
    for state in region_states:
        state_region_mapping[state] = region

for key, dataframe in fires_df.items():
    dataframe['CLIMATE_REGION'] = dataframe['STATE']\
        .apply(lambda x: state_region_mapping[x])

label_encoder = sklearn.preprocessing.LabelEncoder()
for key, dataframe in fires_df.items():
    dataframe['CLIMATE_REGION_CODE'] = label_encoder.fit_transform(dataframe['CLIMATE_REGION'])

def h3_for_chunk(chunk, precision):
    lat_lon_lines = "\n".join([
        "{} {}".format(row['LATITUDE'], row['LONGITUDE']) for index,row in chunk.iterrows()
    ])
    h3 = subprocess.run(
        ['/tools/h3/bin/geoToH3', str(precision)],
        input=str.encode(lat_lon_lines),
        stdout=subprocess.PIPE).stdout.decode('utf-8').splitlines()
    return pd.DataFrame({
        "H3": h3
    }, index=chunk.index)

def chunker(seq, size):
    return (seq.iloc[pos:pos + size] for pos in range(0, len(seq), size))

for key in fires_df:
    print('Calculating bins for fires_df["{}"]'.format(key))
    h3_df = pd.DataFrame()
    with tqdm_notebook(total=fires_df[key].shape[0]) as pbar:
        for chunk in chunker(fires_df[key], 10000):
            h3_df = h3_df.append(h3_for_chunk(chunk, 5))
            pbar.update(chunk.shape[0])
    fires_df[key].drop('H3', 1, inplace=True, errors='ignore')
    fires_df[key] = fires_df[key].merge(
        h3_df,
        how='left',
        left_index=True,
        right_index=True,
        copy=False
    )
    display(fires_df[key].sample(5))

Calculating bins for fires_df["train"]

Calculating bins for fires_df["test"]

bins = gpd.GeoDataFrame({
    'H3': np.unique(np.concatenate((fires_df["train"].H3.unique(), fires_df["test"].H3.unique())))
})

print(bins.info())
bins.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 25999 entries, 0 to 25998
Data columns (total 1 columns):
H3    25999 non-null object
dtypes: object(1)
memory usage: 203.2+ KB
None

def h3_bounds_for_bin(h3_bin):
    h3_bounds = subprocess.run(
        ['/tools/h3/bin/h3ToGeoBoundary'],
        input=str.encode(h3_bin),
        stdout=subprocess.PIPE).stdout.decode('utf-8').splitlines()
    coords = []
    for coord in h3_bounds[2:-1]:
        lon_lat = coord.lstrip().split(" ")
        coords.append(((float(lon_lat[1]) - 360.), float(lon_lat[0])))
    return shapely.geometry.Polygon(coords)

bins['GEOMETRY'] = bins['H3'].apply(h3_bounds_for_bin)

print(bins.info())
bins.sample(10)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 2 columns):
H3          27099 non-null object
GEOMETRY    27099 non-null object
dtypes: object(2)
memory usage: 423.5+ KB
None

bins_land_cover = bins.copy()

land_cover_path = "/data/188-million-us-wildfires/src/nlcd_2011_landcover_2011_edition_2014_10_10/nlcd_2011_landcover_2011_edition_2014_10_10.img"

project = functools.partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj('+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'))

land_cover_data = rasterio.open(land_cover_path)

def render_masked_bin(row):
    return rasterio.mask.mask(
        land_cover_data,
        [shapely.geometry.mapping(transform(project, row['GEOMETRY']))],
        crop=True
    )

masked_bin_images = bins_land_cover.sample(6).apply(render_masked_bin, 1).tolist()

fig, plots = plt.subplots(3, 2)
for index, val in enumerate(masked_bin_images):
    my_plot = plots[math.floor(index/2)][index % 2]
    rasterio.plot.show(val[0], ax=my_plot, transform=val[1], cmap='gist_earth')
    my_plot.axis('on')

land_cover_path = "/data/188-million-us-wildfires/src/nlcd_2011_landcover_2011_edition_2014_10_10/nlcd_2011_landcover_2011_edition_2014_10_10.img"


def land_cover_for_geometry(geometry):
    stats = zonal_stats(
        [shapely.geometry.mapping(transform(project, geometry))],
        land_cover_path,
        nodata=-999,
        categorical=True
    )
    return stats[0]

bins_land_cover = bins.copy()
bins_land_cover["LAND_COVER"] = bins_land_cover["GEOMETRY"].apply(land_cover_for_geometry)

print(bins_land_cover.info())
bins_land_cover.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 3 columns):
H3            27099 non-null object
LAND_COVER    27099 non-null object
GEOMETRY      27099 non-null object
dtypes: object(3)
memory usage: 635.2+ KB
None

def primary_land_cover(land_cover):
    return max(land_cover, key=land_cover.get)

bins_land_cover["PRIMARY_LAND_COVER"] = bins_land_cover["LAND_COVER"].apply(primary_land_cover)

print(bins_land_cover.info())
bins_land_cover.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 4 columns):
H3                    27099 non-null object
LAND_COVER            27099 non-null object
GEOMETRY              27099 non-null object
PRIMARY_LAND_COVER    27099 non-null object
dtypes: object(4)
memory usage: 846.9+ KB
None

land_cover_categories = {
    "water": [11, 12],
    "developed": [21, 22, 23, 24],
    "barren": [31],
    "forest": [41, 42, 43],
    "shrubland": [51, 52],
    "herbaceous": [71, 72, 73, 74],
    "planted_cultivated": [81, 82],
    "wetlands": [90, 95]
}

land_cover_categories_map = {}
for category, values in land_cover_categories.items():
    for value in values:
        land_cover_categories_map[value] = category

def land_cover_category_for_class(land_cover_class):
    return land_cover_categories_map[int(land_cover_class)] if int(land_cover_class) in land_cover_categories_map else 'unknown'

bins_land_cover["PRIMARY_LAND_COVER_CATEGORY"] = bins_land_cover["PRIMARY_LAND_COVER"].apply(land_cover_category_for_class)

print(bins_land_cover.info())
bins_land_cover.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 5 columns):
H3                             27099 non-null object
LAND_COVER                     27099 non-null object
GEOMETRY                       27099 non-null object
PRIMARY_LAND_COVER             27099 non-null object
PRIMARY_LAND_COVER_CATEGORY    27099 non-null object
dtypes: object(5)
memory usage: 1.0+ MB
None

for key, dataframe in fires_df.items():
    label_encoder = sklearn.preprocessing.LabelEncoder()
    bins_land_cover.loc[:,'PRIMARY_LAND_COVER_CODE'] = label_encoder.fit_transform(bins_land_cover['PRIMARY_LAND_COVER'])
    label_encoder = sklearn.preprocessing.LabelEncoder()
    bins_land_cover.loc[:,'PRIMARY_LAND_COVER_CATEGORY_CODE'] = label_encoder.fit_transform(bins_land_cover['PRIMARY_LAND_COVER_CATEGORY'])

print(bins_land_cover.info())
bins_land_cover.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 7 columns):
H3                                  27099 non-null object
LAND_COVER                          27099 non-null object
GEOMETRY                            27099 non-null object
PRIMARY_LAND_COVER                  27099 non-null object
PRIMARY_LAND_COVER_CATEGORY         27099 non-null object
PRIMARY_LAND_COVER_CODE             27099 non-null int64
PRIMARY_LAND_COVER_CATEGORY_CODE    27099 non-null int64
dtypes: int64(2), object(5)
memory usage: 1.4+ MB
None

drop_cols = [c for c in bins_land_cover.columns if c[:15] == "LAND_COVER_PCT_"]
bins_land_cover = bins_land_cover.drop(drop_cols, 1)

def land_cover_percentages(row):
    total = 0
    for lc_class, val in row['LAND_COVER'].items():
        total = total + int(val)
    output = {}
    for lc_class, val in row['LAND_COVER'].items():
        pct = val / total
        row['LAND_COVER_PCT_{}'.format(lc_class)] = pct
    return row

bins_land_cover = bins_land_cover.apply(land_cover_percentages, 1)

bins_land_cover.loc[
    :,[c for c in bins_land_cover.columns if c[:15] == "LAND_COVER_PCT_"]
] = bins_land_cover[[c for c in bins_land_cover.columns if c[:15] == "LAND_COVER_PCT_"]].fillna(0)

print(bins_land_cover.info())
bins_land_cover.sample(5)

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 7 columns):
H3                                  27099 non-null object
LAND_COVER                          27099 non-null object
PRIMARY_LAND_COVER                  27099 non-null object
PRIMARY_LAND_COVER_CATEGORY         27099 non-null object
PRIMARY_LAND_COVER_CATEGORY_CODE    27099 non-null int64
PRIMARY_LAND_COVER_CODE             27099 non-null int64
GEOMETRY                            27099 non-null object
dtypes: int64(2), object(5)
memory usage: 1.4+ MB
None

for key, dataframe in fires_df.items():
    fires_df[key] = dataframe.merge(bins_land_cover, on="H3", how="left")
    display(fires_df[key].sample(5))

census_data = pd.DataFrame(
    {"H3": bins['H3'].copy()},
    index=bins.index
)

census_files = {
    "Population Density": "/data/188-million-us-wildfires/src/usgrid_data_2010/uspop10.tif",
    "Seasonal Housing": "/data/188-million-us-wildfires/src/usgrid_data_2010/ussea10.tif",
    "Housing Units": "/data/188-million-us-wildfires/src/usgrid_data_2010/ushu10.tif",
    "Occupied Housing Units": "/data/188-million-us-wildfires/src/usgrid_data_2010/usocc10.tif"
}

def render_bin(row):
    output = []
    for index, value in census_files.items():
        census_data = rasterio.open(value)
        output.append(rasterio.mask.mask(
            census_data,
            [shapely.geometry.mapping(row['GEOMETRY'])],
            crop=True,
            nodata=0
        ))
    return output

masked_bin_images = bins_land_cover.sample(2).apply(render_bin, 1).tolist()

fig, plots = plt.subplots(4, 2)
for index, val in enumerate(masked_bin_images):
    for subindex, image in enumerate(val):
        my_plot = plots[subindex][index]
        rasterio.plot.show(image[0], ax=my_plot, transform=image[1], cmap='inferno')

census_population_path = "/data/188-million-us-wildfires/src/usgrid_data_2010/uspop10.tif"

def mean_population_density_for_geometry(geometry):
    stats = zonal_stats(
        [shapely.geometry.mapping(geometry)],
        census_population_path,
        nodata=0,
        stats=["mean"]
    )
    return stats[0]["mean"]

census_data = census_data.drop("MEAN_POPULATION_DENSITY", 1, errors='ignore')
census_data["MEAN_POPULATION_DENSITY"] = \
    bins["GEOMETRY"].apply(mean_population_density_for_geometry).fillna(value=0).replace([-np.inf], 0)

print(census_data.info())
census_data.sample(5)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 2 columns):
H3                         27099 non-null object
MEAN_POPULATION_DENSITY    27099 non-null float64
dtypes: float64(1), object(1)
memory usage: 423.5+ KB
None

census_seasonal_path = "/data/188-million-us-wildfires/src/usgrid_data_2010/ussea10.tif"

def mean_seasonal_housing_for_geometry(geometry):
    stats = zonal_stats(
        [shapely.geometry.mapping(geometry)],
        census_seasonal_path,
        nodata=0,
        stats=["mean"]
    )
    return stats[0]["mean"]

census_data = census_data.drop("MEAN_SEASONAL_HOUSING", 1, errors='ignore')
census_data["MEAN_SEASONAL_HOUSING"] = \
    bins["GEOMETRY"].apply(mean_seasonal_housing_for_geometry).fillna(value=0).replace([-np.inf], 0)

print(census_data.info())
census_data.sample(5)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 3 columns):
H3                         27099 non-null object
MEAN_POPULATION_DENSITY    27099 non-null float64
MEAN_SEASONAL_HOUSING      27099 non-null float64
dtypes: float64(2), object(1)
memory usage: 635.2+ KB
None

census_housing_units = "/data/188-million-us-wildfires/src/usgrid_data_2010/ushu10.tif"

def mean_housing_units_for_geometry(geometry):
    stats = zonal_stats(
        [shapely.geometry.mapping(geometry)],
        census_housing_units,
        nodata=0,
        stats=["mean"]
    )
    return stats[0]["mean"]

census_data = census_data.drop("MEAN_HOUSING_UNITS", 1, errors='ignore')
census_data["MEAN_HOUSING_UNITS"] = \
    bins["GEOMETRY"].apply(mean_housing_units_for_geometry).fillna(value=0).replace([-np.inf], 0)

print(census_data.info())
census_data.sample(5)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 4 columns):
H3                         27099 non-null object
MEAN_POPULATION_DENSITY    27099 non-null float64
MEAN_SEASONAL_HOUSING      27099 non-null float64
MEAN_HOUSING_UNITS         27099 non-null float64
dtypes: float64(3), object(1)
memory usage: 846.9+ KB
None

census_occupied_housing_units = "/data/188-million-us-wildfires/src/usgrid_data_2010/usocc10.tif"

def mean_occupied_housing_units_for_geometry(geometry):
    stats = zonal_stats(
        [shapely.geometry.mapping(geometry)],
        census_occupied_housing_units,
        nodata=0,
        stats=["mean"]
    )
    return stats[0]["mean"]

census_data = census_data.drop("MEAN_OCCUPIED_HOUSING_UNITS", 1, errors='ignore')
census_data["MEAN_OCCUPIED_HOUSING_UNITS"] = \
    bins["GEOMETRY"]\
        .apply(mean_occupied_housing_units_for_geometry)\
        .fillna(value=0)\
        .replace([-np.inf], 0)

print(census_data.info())
census_data.sample(5)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 27099 entries, 0 to 27098
Data columns (total 5 columns):
H3                             27099 non-null object
MEAN_POPULATION_DENSITY        27099 non-null float64
MEAN_SEASONAL_HOUSING          27099 non-null float64
MEAN_HOUSING_UNITS             27099 non-null float64
MEAN_OCCUPIED_HOUSING_UNITS    27099 non-null float64
dtypes: float64(4), object(1)
memory usage: 1.0+ MB
None

census_data["MEAN_VACANT_HOUSING_UNITS"] = census_data["MEAN_HOUSING_UNITS"] - census_data["MEAN_OCCUPIED_HOUSING_UNITS"]

for key, dataframe in fires_df.items():
    fires_df[key] = fires_df[key]\
        .drop([
            "MEAN_POPULATION_DENSITY",
            "MEAN_SEASONAL_HOUSING",
            "MEAN_HOUSING_UNITS",
            "MEAN_OCCUPIED_HOUSING_UNITS",
            "MEAN_VACANT_HOUSING_UNITS"
        ], axis=1, errors='ignore')
    fires_df[key] = dataframe.merge(census_data, on="H3", how="left")
    display(fires_df[key].sample(5))