# ADS 2.0: new architecture, API and services

Roman Chyla, Alberto Accomazzi, Alexandra Holachek, Carolyn S. Grant, Jonathan
Elliott, Edwin A. Henneken, Donna M. Thompson, Michael J. Kurtz, Stephen S.
Murray, Vladimir Sudilovsky  Harvard-Smithsonian Center for Astrophysics, 60
Garden Street, Cambridge, MA 02138, USA

###### Abstract

The ADS platform is undergoing the biggest rewrite of its 20-year history.
While several components have been added to its architecture over the past
couple of years, this talk will concentrate on the underpinnings of ADS’s
search layer and its API. To illustrate the design of the components in the
new system, we will show how the new ADS user interface is built exclusively
on top of the API using RESTful web services. Taking one step further, we will
discuss how we plan to expose the treasure trove of information hosted by ADS
(10 million records and fulltext for much of the Astronomy and Physics
refereed literature) to partners interested in using this API. This will
provide you (and your intelligent applications) with access to ADS’s
underlying data to enable the extraction of new knowledge and the ingestion of
these results back into the ADS. Using this framework, researchers could run
controlled experiments with content extraction, machine learning, natural
language processing, etc. In this talk, we will discuss what is already
implemented, what will be available soon, and where we are going next.

##  1 Introduction

The SAO/NASA Astrophysics Data System (Kurtz et al. 2000) has gone through
significant changes over the course of its 20 year career, but its purpose has
always been to provide relevant information and services to the scientific
community. The ADS has achieved this goal through constant innovation and
adoption of new technologies. In this paper we will describe the changes that
are happening to the ADS in its recent incarnation. They constitute the
biggest rewrite in the history of the system and mark the opening of many new
possibilities, particularly in discovery capabilities and support for
bibliometric studies.

##  2 Environment in which the ADS operates

The ADS is an _information system_ \- built by and for scientists to provide
access to scientific literature. Its user base is considerable, currently at
any given time there are $\sim 200-300$ active users on the website (which
amounts to $\sim 1$ million unique sessions per month), although its core user
base is much smaller (Kurtz et al. 2009). For the majority of these
(anonymous) users worldwide the ADS is the destination they reached after
following a link from Wikipedia or Google to one of its records. However for
about 50,000 scientists worldwide, the ADS represents a gateway to
information, they use the system every day (or every other day). And for the
ADS they are the users who really care and need to be cared for.

The ADS provides its services within a greater information space environment,
which includes several service providers. One such providers is the arXiv,
arguably one of the most indispensable academic services, at least in the
fields of astronomy and physics. It has become a de facto publishing standard
to post pre-prints on the arXiv before they are available as published and
reviewed papers through more traditional channels. However, as many things
that we take for granted and expect them to be there, the arXiv operates with
a very limited budget and resources.

The arXiv (and its volunteers, the scientists) performs the hard work of pre-
filtering submissions and categorizing them. It does serve one purpose and
serves it well, yet it needs other projects to complement it. One such project
is the ADS, which closely collaborates with arXiv, downloads and indexes every
single arXiv submission - it is perhaps a little known fact, but you can use
the powerful search capabilities of the ADS to sift through all arXiv
fulltexts.

Throughout its existence, the ADS has established such symbiotic relationships
in order to provide the best service to the scientific community. It
collaborates not only with arXiv, but also with data repositories, publishers,
and people behind production of journals and conferences, it has linked pre-
prints to the published papers and other grey-literature, so that they not
disappear into obscurity. Tracking the important channels of scientific
information logically led to demands for indicators of scientific contribution
(the citation indicators) - which in turn get (mis)used as metrics of
scientific success. Despite how we may look at it, these indicators remain
important tools and the ADS recognizes their importance. In addition to
scientists, who come in their roles of researchers and readers, there are also
administrators and evaluators, librarians, as well as technocrats and grant
agencies who use the ADS database for other, non-scientific purposes. (Kurtz
et al. 2005; Accomazzi et al. 2007)

The reason we mention these diffeent user populations is to explain the
context in which the ADS operates. It has to choose wisely which services are
feasible to implement and maintain. Given limited resources, this becomes not
only a technological, but also a sociological challenge.

The discussion becomes simpler if we consider the entities (objects) processed
by the ADS. As an information system, the ADS primarily processes _scientific
literature_ , so the most frequent data types that it deals with are:

  * •

_Preprints_ and _journal articles_ : The ADS will store them in the database
(their _fulltext_ content as well as attached objects such as _extracted
images_).

  * •

_Metadata_ : Metadata can be very expensive to create and curate but its use
may be surprisingly powerful. For example, if author affiliations are known,
one can compile cumulative reports about the scientific output of a certain
institute, department and/or university (even country or continent).

  * •

_Citation information_ : Citing and cited by are very powerful pieces of
information, they are the basis of the citation index and the metrics
available throughout the ADS. Its use in the ADS is so important that in the
new version of the search engine it has been incorporated into the index, so
you can search for citations in the same way as you search for papers.

  * •

_Usage information_ : This is information about _readers_ , what they are
reading, what papers are popular, downloaded, consulted and how often.111This
information is properly anonymized before it gets used by different algorithms
and subsystems of the ADS.

Finally, there is _raw data_ , but in case of the ADS rather _Links to science
data_. These consists of a network of links between bibliographic records
hosted by ADS and data products, astronomical object, and electronic articles
hosted elsewhere. Other projects (such as NED and the CDS) specialize in
collecting this information and do a much better job with it than what the ADS
could do. The role of the ADS here is to serve as a gateway and make these
resources seen and used.

And what we mean by saying that is not just people who sit behind computer
screens looking at the ADS web interface - to the ADS developers those people
are morphing with the machines (robots) that start to access the new ADS
infrastructure. This has nothing to do with de-humanization or other negative
connotations of information processing. As far as information systems are
concerned, it has become impractical to treat users differently from machines.
The differences between the humans and the clients (robots, programs,
machines) that access the ADS on the users’ behalf are disappearing, and this
brings us to the new era and new architecture of the ADS. Since a few years,
the system has been preparing itself for a change and it is not exactly
_evolution_ but rather _revolution_.

##  3 New architecture and the API

While the old infrastructure of the ADS is similar to a huge building which
hosts every department and which has only one entrance (and only one power
supply), the new infrastructure is much more distributed and lightweight, more
similar to a town with roads connecting different buildings. The central
component here is the _Application Programming Interface (API) layer_ 222The
API layer is itself published on GitHub, at http://github.com/adsabs/adsws \-
we have split the internal subsystems into _microservices_. They are
independent, standalone web-services that communicate with each other through
the API using the REST and OAuth protocols (and soon also exclusively through
HTTPS). Every service, even the ones that we consider the most critical, are
connected to the API and available through it; can be accessed from anywhere
on the internet provided that the client has the proper authorization token.

The micro-services are virtualized, currently running off the Amazon cloud.
They can be scaled horizontally: if there is need for more power, we add
computational nodes to the cloud. Using this new organization, the ADS can
more easily add new functionality (new services), and it is relatively easy to
alter existing services. Basically, implementations can be swapped at will,
provided the semantics of the API remains unchanged. It is still challenging
to find the absolute minimal set of the API, but the system has already proven
itself. It simplifies maintenance and allows for rapid development. 333As will
be discussed later, this has also a negative consequences; at the moment, the
ADS is changing so rapidly, that is hard to keep the documentation up to date.
It requires certain mental discipline and robust practices and we are still
struggling to find the optimum balance.

###  3.1 The search

The core of the new ADS (similarly to the previous generation of the ADS) is
centered around search. The search engine is built on a heavily modified and
extended version of Apache SOLR (http://lucene.apache.org/solr/). As one would
probably expect, it has many advantages over the in-house custom search
service developed for ADS Classic 20 years ago. Unsurprisingly, it is able to
address much bigger memory spaces which translates into an ability to index
bigger volumes of data. For the first time in ADS history, it is possible to
search through fulltext content together with the metadata. Additionally, we
have integrated the citation and co-readership networks into the index, so you
can now simultaneously explore all three different search spaces: metadata,
fulltext content, and citations.

The new search engine provides many features, some of which will be totally
new to the ADS users, for example:

  * •

proximity search – ability to find a group of tokens based on their distance,
e.g. body:(weak NEAR5 lensing), which finds phrases where “weak” appear within
5 words of “lensing.” This query would find phrase ”weak lensing”, as well as
”weak gravitational lensing” in the fulltext

  * •

similarity search – (e.g. author:eisenstein$\sim 0.3$) will find different
spellings of the author name based on different similarity metric (this can be
used to correct for misspellings)

  * •

regular expressions – most often used by curators to find terms that follow a
particular pattern e.g. grant numbers or astronomical entities

The features described above are provided by the Lucene search library upon
which SOLR is built. In addition, ADS has extended Lucene capabilities by
implementing its own _query language_ , modifying the traditional search
syntax and enhancing it with new search operators and modifiers. The main
purpose of this “dialect” is to facilitate communication between clients and
servers. We do not expect normal ADS users to know the specialized syntax, but
we want to give them tools to execute complex queries if they ever have a need
to do so. More likely though, this query language will be used by clients
(programs) that search the ADS on behalf of users. To illustrate the point,
here is a short list of a few specialized modifiers/operators:

  * •

=searchterm – the equals sign modifier is perhaps known to ADS users, it
disables the synonym expansion so that the search engine will not look for
synonyms or alternate spellings.

  * •

references(author:einstein) and citations(author:einstein) – these are the
operators to retrieve papers that cite a collection of papers, or those that
are cited by that collection. Operators can be freely nested and combined with
other terms, even if that means that you are going to analyze the entire ADS
corpus, e.g. citations(references(author:einstein))

  * •

useful(), instructive(), trending() and other operators – the ADS search
engine supports so-called second order queries which use the citation and
readership network to provide advanced search capabilities. Currently these
operators use multiple algorithms to retrieve papers that are most cited by
papers on a topic (useful), citing the most relevant papers on a topic
(instructive), or most read by people interested on a topic (trending). The
ADS will be experimenting with different implementations of these operators
and will be adding more of them as needed.

The citation operators together with custom analytical functions are
potentially the most fun part of the new search engine. It is possible to
construct custom searches just for the purpose of supporting certain
visualizations (e.g. a view of a co-readership network) or use cases (e.g.
selecting the _N_ th token in an list of author names). In the future, we
expect users of the API to come up with novel ways of using them, something we
did not anticipate when they were created. 444Of course we are aware that in
order for this to happen we need to provide better documentation of the search
functionality.

###  3.2 Other services

The best example of the API in action is the new user interface for the ADS,
code-named Bumblebee. It is a rich JavaScript client that runs in your web
browser. It is the future replacement of the ADS Abstract service, but unlike
that one, virtually every interaction happens through the API. So without
knowing it, when you visit http://ui.adslabs.org, you are automatically
becoming an API user and accessing multiple web services. It will not feel any
different, but the difference in implementation is significant: to have access
to the API means that you are a first-class citizen of the ADS world,
accessing the same interface that ADS developers use every day.

Figure 1.: Screenshot of the new UI for the ADS, using normal search.

At the time of this writing, there exist multiple services besides search:

  * •

_metrics_ : the familiar citation and usage reports (see 2), if you access the
API directly, you can download the data that are used to generate the metrics
views (and build your own version of the visualization)

  * •

_visualizations_ : this service exposes different visualizations based on the
relations between authors and papers, e.g. word-clouds, paper networks,
citation networks, and bubble-charts of search results (see 3).

  * •

_ORCID_ : service that allows users to claim authorship through the ADS
interface into the ORCID database555ORCID stands for Open Researcher and
Contributor ID: http://orcid.org

Information on how to access and use the Api for these services, and any
future ones, is available through: http://adsabs.github.io/help.

Figure 2.: Citation metrics report for the most cited 100 papers of A.
Einstein. Figure 3.: Paper networks inside the first 100 most cited papers of
A. Einstein.

##  4 Future directions

We are aware that the new API has a number of shortcomings especially with
respect to documentation. It is still a very young and relatively new
component and it may not be a smooth ride for external developers, but the ADS
has already gone a long way in making this new ”mode” available and will
continue improving it. In the past year, the ADS has re-engineered many of its
subsystems from the ground up (after a few false starts). The current
implementation has shown that the system architecture is robust, and this
gives us confidence to keep going in this direction. To offer a level of
stability to the platform, major changes to the API will be versioned so that
we will be able to provide backward compatibility to applications built on it.
Our hope is that ultimately this open API platform will attract development
from outside the ADS group, leveraging potential hidden in wider community.

The new generations of astronomers and data scientists are technically very
savvy and ready to embrace emerging technologies to conduct their research. We
hope that they will appreciate the level of access the new API provides, and
that they will take the opportunity to build something that satisfies their
research needs (and share it with other people). We hope that somebody will
come up with ideas we did not have the time or acuity to consider, new ways to
explore the vast data indexed in ADS, or new ways to connect existing systems
and give insights that were not previously possible. We will continue building
and adding more services, in the near future improving the relevancy of the
search, and preparing the new user interface to replace the existing ADS
Abstract Service. Once the new user interface is rolled out, any new
contributions to it will become visible to at least a million users every
month, which will hopefully provide enough usability data to test even more of
the new ideas. Intrigued and want see more or perhaps even build something?
All of the code is published open-source and available through ADS’s Github
repo: http://github.com/adsabs. We welcome any and all comments, criticisms
and contributions.

#### Acknowledgments

This work has been supported by the NASA Astrophysics Data System project,
funded by NASA grant NNX12AG54G.

## References

  * Accomazzi et al. (2007) Accomazzi, A., Eichhorn, G., Kurtz, M. J., Grant, C. S., Henneken, E., Demleitner, M., Thompson, D., Bohlen, E., & Murray, S. S. 2007, in Library and Information Services in Astronomy V, edited by S. Ricketts, C. Birdie, & E. Isaksson, vol. 377 of Astronomical Society of the Pacific Conference Series, 69. cs/0610011
  * Kurtz et al. (2009) Kurtz, M., Accomazzi, A., & Murray, S. S. 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey, vol. 2010 of Astronomy, 28P. 0903.3228
  * Kurtz et al. (2000) Kurtz, M. J., Eichhorn, G., Accomazzi, A., Grant, C. S., Murray, S. S., & Watson, J. M. 2000, A&AS, 143, 41. astro-ph/0002104
  * Kurtz et al. (2005) Kurtz, M. J., Eichhorn, G., & Accomazzi, A. e. a. 2005, Journal of the American Society for Information Science and Technology, 56, 111. 0909.4789


# Wildland Fire Mid-story: A generative modeling approach for representative
fuels ††thanks: Los Alamos National Laboratory LA-UR-23-25924

Grant Hutchings  Statistical Sciences, Los Alamos National Laboratory  James
Gattiker  Statistical Sciences, Los Alamos National Laboratory  Braden
Scherting  Statistical Sciences, Los Alamos National Laboratory

(June 28, 2023)

###### Abstract

Computational models for understanding and predicting fire in wildland and
managed lands are increasing in impact. Data characterizing the fuels and
environment is needed to continue improvement in the fidelity and reliability
of fire outcomes. This paper addresses a gap in the characterization and
population of mid-story fuels, which are not easily observable either through
traditional survey, where data collection is time consuming, or with remote
sensing, where the mid-story is typically obscured by forest canopy. We
present a methodology to address populating a mid-story using a generative
model for fuel placement that captures key concepts of spatial density and
heterogeneity that varies by regional or local environmental conditions. The
advantage of using a parameterized generative model is the ability to
calibrate (or ‘tune’) the generated fuels based on comparison to limited
observation datasets or with expert guidance, and we show how this generative
model can balance information from these sources to capture the essential
characteristics of the wildland fuels environment. In this paper we emphasize
the connection of terrestrial LiDAR (TLS) as the observations used to
calibrate of the generative model, as TLS is a promising method for supporting
forest fuels assessment. Code for the methods in this paper is available.

Keywords: Wildfire Modeling, Cox Process, Bayesian Model Calibration, Gaussian
Process, Prescribed Fire, Environmental Assessment, Spatial Generative Model

##  1 Introduction

Wildland fire has historically been modeled as a physically-inspired and
empirically calibrated spread rate model applied to bulk fuel assessments
(Rothermel, 1972, Andrews, 2018). Computational models of wildland fire have
been explored, and range from fire perimeter spread approximations to
spatially resolved models of fuels, atmospheric dynamics, and physical
phenomenology (Linn et al., 2002, 2005, 2020). These computational simulation
models have been employed extensively in the study of the behavior and impacts
of wildfire, for example Linn et al. (2012). Modeling and simulation of
prescribed fire introduces additional challenges. In this report, we focus on
the need for detailed layout of fuels, which is a key aspect of achieving
reliable predictions of fire outcome (Linn et al., 2020, 2021, Hiers et al.,
2009).

Plot-level summary has been the basis in the past for characterization of mid-
story fuels, that is, the fuels under the main forest canopy. For specific
prescribed fire planning, direct survey consistent with the protocols of the
Forest Inventory and Analysis (FIA) program (Toney et al., 2007) are currently
used in operational settings. The use of reference plots to impute broad-area
coverage (Riley et al., 2016) has made estimated plot-level statistics
available with broad coverage. Manual survey is time-consuming and expensive
and summarizes plots of hundreds of square meters in a relatively small number
of metrics, sufficient for some purposes of large-scale environmental
monitoring (Tinkham et al., 2018). However, the spatial resolution of
technologies for analyzing the interactions between forest composition and
wildland fire now exceeds the spatial resolution sought by typical existing
survey protocols. Computational models such as FIRETEC accept input on the
meter scale, currently, allowing for the specification of explicit tree
positions, and spatial layout of heterogeneous mid-story fuels. There is
growing cognisance of the quantitative connection between heterogeneity in
fuels to fire outcomes Knapp and Keeley (2006), Linn et al. (2013), and the
influence of small-scale (meter to tens of meters) heterogeneity in fuel
structure on salient fire behavior (Parsons et al., 2017, Atchley et al.,
2021). Further, prescribed fire planning targets marginal conditions, where
outcomes have been shown to be sensitive (Jonko et al., 2021). Survey
protocols leveraging remote sensing technologies are emerging to satisfy the
need for spatial resolution in fuels for wildland fire.

Light detection and ranging (LiDAR) has found a strong and growing role in
environmental characterization. Areal LiDAR survey (ALS) achieves resolution
sufficient to inventory individual trees (e.g., Li et al. (2012)). This,
coupled with the large spatial extent of ALS deployed from aircraft, enables
exact representation of overstory fuels for large swaths, although currently
not for complete coverage over time. Recent work examines the ability of _in-
situ_ terrestrial LiDAR scanning (TLS) to estimate plot-level summaries
(Pokswinski et al., 2021), including tree identification and bulk
characterization of mid-story (Anderson et al., 2021, Silva et al., 2016,
Rowell et al., 2020, Loudermilk et al., 2023). This report will describe
results based on information extracted from TLS, but is agnostic to the method
of mid-story observation, which could be based on manual survey, remote or
_in-situ_ sensing imagery, multi-modal LiDAR, or other means.

Trees can be automatically extracted from broad-scale remote sensing survey,
and this abstraction of trees and canopy can be used to populate the
corresponding fuels in simulation models. Recent work (McDanold et al., 2023)
develops an approach to generating grass and litter surface fuels using a
physically-motivated model that includes dependence on overstory and canopy.
The mid-story is not typically characterized in a complete format that allows
exact layout, although it is correlated with canopy and other spatial
covariates. The perspective taken in this report is that mid-story fuels will
be characterized from limited relevant observations of an ecosystem, and
generated for the purposes of analysis and simulation models as
_representative_ , rather than exact, fuel distribution.

This report addresses the interpretation and generation of heterogeneous mid-
story fuels with a spatial statistical modeling approach. The primary goal is
the definition and demonstration of a statistical model that can generate
representative mid-story fuels. A spatial statistical model can be calibrated
to observed data; that is, the settings of the model that make it most
consistent with observations can be discovered, perhaps also informed by prior
constraints on the settings that control the characteristics of the generated
spatial patterns. There is a large family of spatial statistical models that
might be used to define generative processes. Setting reasonable values for
the parameters of any generative model is known to be a considerable
challenge, especially as the model grows in complexity. This work takes the
position of investigating a model that is complex enough to capture features
of heterogeneous fuels distributions, while simple enough to infer parameters
of the model proposed given observations, or a practical combination of
expertise and observational constraints.

The next section will describe the family of methodology that will be applied
for representing and generating heterogeneous fuels layout. The approach to
model calibration is then presented, along with a description of the metrics
that will be used in the calibration process. A review of our datasets and
goals leads to a demonstration validating the ability to infer generative
model settings using instances from the model to support an understanding and
validation of the capability in principle, and finally the corresponding
process of determining model settings from observed data and the application
of the method to fuels generation is shown. We conclude by describing how this
model can fit into a framework of fuel layouts including additional spatial
covariates, specifically the dependence of mid-story to canopy and to urban
features, i.e. roads.

##  2 Point-Process Model

The goal of this work is to demonstrate a model that can generate synthetic
layouts of mid-story fuels that are representative of an environment’s
characteristics. By _representative_ fuels, we mean that the model output is
similar to the target environment in established metrics that are chosen to
inform characteristics of interest in fuels layouts, including structure and
heterogeneity. We now describe some background of spatial models that lead to
the capability of generating representative fuels.

Spatial statistical models have a rich history in the description and analysis
of environmental data (Banerjee et al., 2003). The methodology here comes from
two strands of this literature: spatial point processes and binary models.

###  2.1 Spatial Point Processes

There is a deep literature in point processes, this section will introduce the
topic relevant to the application. Data that contains the locations of events
of interest (in space and/or time) are known as a point pattern. The standard
term _event_ corresponds here to the siting of a fuel element, which from
different points of view could represent a simple spatial sample of fuel, or
be interpreted as distinct features such as shrubs. Suppose we observe the
locations of $n$ plants $\\{s_{i}\\}_{n}$ in some plot or domain
$D\in\mathbb{R}^{2}$. The natural model for these data is a point process
$\mathcal{P}$, and we write $\\{s_{i}\\}_{n}\sim\mathcal{P}$. A commonly used
point process is the Poisson process, with spatial variability controlled by
an intensity function $\gamma(\mathbf{s})$. The intensity function determines
the density of points throughout the domain. Specifically, the expected number
of points in a region $B\subset D$ (i.e, “random counting measure”) is
$\int_{B}\gamma(\mathbf{s})d\mathbf{s}$. When the intensity is constant,
$\gamma(\mathbf{s})=\gamma$, we obtain a homogeneous Poisson process (HPP);
the more general nonhomogeneous Poisson process (NHPP) is controlled by a non-
constant intensity function; see Møller and Waagepetersen (2007) for further
explanation and examples.

###  2.2 Binary Mosaic Models

The mid-story fuel spatial pattern can be represented as a binary mosaic model
(BMM). In this approach, the observational data and corresponding model output
are a 2-dimensional binary present/absent pattern over the domain. Here, the
BMM is defined as a point process laying out the centers of disks, with the
disk radii following a distribution. The layout of centers is spatially
heterogeneous according to an underlying intensity distribution.

Approximating the map of mid-story fuels by a binary continuum constructed
from the union of overlapping compact sets, for example in Fig. 6, is known
generally as a germ-grain model: germs (centers) follow a stochastic point
process, and grains are given by compact sets which emanate from the germs and
may depend on the point process. The boolean model is special case of the
germ-grain model that arises when the germs follow a homogeneous Poisson
process, parameterized by a constant intensity $\gamma$, and the compact sets
are mutually independent (Molchanov, 1997) forming the binary mosaic. In
previous work, inference on this model will typically target three parameters:
the first two moments of the radii distribution and the scalar intensity
value, which is assumed to be constant across the domain. This model can be
extended to incorporate spatial heterogeneity through the use of an NHPP on
the germs, yielding a spatially-indexed intensity surface $\gamma(s)$. Such a
surface may be a parametric or nonparametric function of spatial locations.
Nonparametric NHPPs can be difficult to estimate, but they obviate the
problematic task of positing a parametric model for intensity a priori. See
ch. 8 of Banerjee et al. (2003) for an overview.

Pielou (1964) was perhaps the first to draw attention to the use of binary
mosaics in ecology. However, it was Diggle (1981) who first proposed an
explicit generative model for binary mosaics in vegetation coverage map
settings. Using a now widely-studied binary mosaic map of heather shrubs, he
proposed a minimum contrast estimation strategy for inferring the homogeneous
Poisson intensity and three Weibull parameters controlling radii. This
analysis succeeded in estimating the parameters and generating realizations
from the fitted models.

Inference on the boolean model for binary mosaics is expanded by Møller and
Helisová (2010). This simulation-based likelihood inference method relies on a
pre-specified Poisson process; the authors consider only homogeneous
intensities. Our approach extends both results by permitting the underlying
process intensity to vary over the domain, thereby enabling the modeling of
more complicated vegetation structures. The approach is naturally extended to
allow intensities that vary as a function of observed covariates, that is,
spatial maps that indicate higher or lower expectation for observing fuel.

###  2.3 Spatial Generative Model

We next present the details of the non-homogeneous Boolean process used to
model the binary map, which is in the Gaussian Cox process family. The model
is presented in component-wise manner that reflects the generative process;
deviations between the formal statistical model and generative/sampling
procedures are noted where relevant.

The specification of a union of disks used to approximate the binary map
requires the specification of locations of disk centers and the radii of the
respective disks. We first specify the model for the number of disks $n$ and
their respective center locations $s_{1},\dots,s_{n}$, henceforth ”points”.
The total number of points is sampled

$n\sim Pois(\lambda*A_{D}),$

where $A_{D}$ is the area of the domain and $\lambda$ represents points per
unit area. We assume the layout of points in the domain $D$ is controlled by a
smooth and continuous relative intensity function
$\omega(s):D\to\mathbb{R}^{+}$. This is achieved by modeling $\omega$ as an
appropriately-transformed Gaussian Process (GP),

$W(s)\sim\mathcal{GP(\mathbf{0},C_{\rho})}$ (1) $\gamma(s)=f(W(s))$ (2)

where $C_{\rho}$ is a covariance matrix generated by an appropriate covariance
function with hyperparameter $\rho$, $f$ is a function mapping from
$\mathbb{R}\to\mathbb{R}^{+}$. Taking $f(s)=\exp(s)$ corresponds to a Log-
Gaussian Cox process.

The proposed generative process takes $f(s)=logistic(s)$, which allows the
interpretation of the relative intensity function as a probability. We model
the relative intensity function $\omega(s)$, rather than the intensity
function $\gamma(s)$ for convenience in the generative procedure. For the
relative intensity function, $E[n_{B}]\neq\int_{B}\omega(s)ds$, rather
$\forall s\in D,\;P(s_{i}=s)=\omega(s)$.

Given a value $\rho=\hat{\rho}$, a process realization $W$ may be sampled from
the GP. A process realization $W$ is function-valued, but a discrete
approximation is obtained by computing the covariance function on discrete
grid $d\times d$, yielding covariance matrix $\Sigma_{\hat{\rho}}$, and
sampling from the multivariate normal distribution
$\mathcal{N}(\mathbf{0},\Sigma_{\hat{\rho}})$. Transforming this draw by $f$
gives the intensity function. Conditional on the relative intensity function,
the joint density of the number of points and point locations is

$\displaystyle s_{1},\dots,s_{n},n$
$\displaystyle\sim\prod_{i}\omega(s_{i})\times\frac{(\lambda
A_{D})^{n}\exp(-\omega A_{D})}{n!}.$ (3)

The points are collected in the set $\Psi=\\{s_{1},\dots,s_{n}\\}$, which is a
sample from the NHPP.

We assume the locations and sizes are independent, so disk radii may be
sampled from a distribution on $\mathbb{R}^{+}$. It is arguable that a spatial
dependence exists for disk radii, but we will show that the inference process
including a single spatial field can be challenging. Adding complexity to this
model introduces significant challenges in inference, and information at this
level of detail is not typically available for prior specification of
distribution. We choose the truncated Normal distribution bounded below by
zero:

$\displaystyle r_{1},\dots,r_{n}\sim N^{+}(\mu,\sigma^{2})$ (4)

Conditional on $r_{i}$, disk $i$ is the set
$\xi_{i}=\\{s\in\mathbb{R}^{2}:||s-s_{i}||\leq r_{i}\\}$ and the binary map is

$\displaystyle\Xi$ $\displaystyle=\bigcup_{i:s_{i}\in\Psi}\xi_{i}.$

##  3 Generative Model Calibration

This section describes operational details of calibrating the generative
model, including: generating spatial realizations from a model with specified
parameters, metrics used for summarizing and ultimately comparing a spatial
pattern, and the discovery of model parameters using the Bayesian model
calibration framework.

Generative modeling implements a stochastic map from unobservable parameters
and observable covariates to data. These models are useful for encoding an
approximation to the true data generating mechanism, systematically
characterizing uncertainty, and are able to simulate new data. They also
enable specification and estimation of parameters that control general data
features relevant to the research question. This is useful when we wish to
generate synthetic data with particular properties, but wish to integrate over
or ignore other characteristics such as specific location or orientation.
Developing appropriate feature sets capturing, or ignoring, qualities of the
data is an application driven development process.

###  3.1 Generating Model Realizations

A stochastic outcome of the model, that is, a spatial layout, is referred to a
_realization_ of the model conditional on the parameter settings. In the
generative setting, we use the Gaussian Process model intensity function in
Eq. 1, referenced on a $d\times d$ grid $\mathcal{X}_{D}$, to predict the
intensity function at randomly proposed test points. To produce a binary map,
candidate points are proposed and accepted or rejected based on the predicted
intensity function. Given a candidate location $\tilde{s}=(x,y)$ drawn
uniformly on $D$, the intensity function value
$W(\tilde{s})|\mathcal{X}_{D},\mathcal{C}_{\rho}$ is predicted from the
Gaussian Process model. The point $\tilde{s}$ is accepted with probability
$p=\omega(\tilde{s})=f(W(\tilde{s})|\mathcal{X}_{D},\mathcal{C}_{\rho})$. This
process is repeated until $n$ candidate points are accepted, where $n$ is a
Poisson random variable. Each point is the center of a disc of points in the
BMM. The disc radius will be drawn from the truncated Normal distribution of
Eq. 4. The outcome is the BMM shown, for example, in Fig. 2, where the top row
shows realizations of different parameters drawn from their ranges. In
practice, $d$ should be large enough such that the points in $\mathcal{X}_{D}$
are able to capture the heterogeneity in the field observations.

###  3.2 Descriptive metrics

Observational data and model output can be quantified and compared through
extracted metrics describing features of the data relevant to the application,
in this case describing spatial qualities such as area, perimeter, and
topological features such as counts of connected components and holes. We
summarize the geometric geospatial object by a vector of relevant summaries of
the binary mosaic, $y=\eta(\Xi)$.

The metrics used are described in Table 1, with some comments on procedural
details. This is not intended to be a final or exhaustive list of useful
features, but rather a sufficient set for demonstration. Extensions are
possible; one potentially useful example is metrics related to transect
segment summaries from manual survey.

area |  From a binary image, observed or derived from observations, area is the proportion of occupied pixels in the domain.
For a generated representation composed of disks on the domain, we use a Monte
Carlo estimate for the area of the union of the $n$ disks in the generated
data; this provides the flexibility to compute local areas and statistics of
spatial auto-correlation described below.
---|---
perimeter |  From a binary image, observed or derived from observations, perimeter is the number of occupied pixels adjacent to unoccupied pixels, with a Pythagorean correction for boundary pixels that abut occupied pixels diagonally.
For a generated representation, the following approximation is used: Arrange
points $n_{i}$ along the circumference of each disk; discard any circumference
points that lie within another disk; calculate the approximate perimeter as
the proportion of points retained times the analytic total disk perimeters. A
circumference point lies within another disk if the distance between the point
and the center of another disk, $s_{j}$, is less than $r_{j}$.
NCC |  Number of distinct connected components is easily calculated from raster data.
For a generated representation, we construct a graph object using the _igraph_
R package. Graph nodes are given by disk centers; edge $e_{ij}=1$ if
$d(s_{i},s_{j})\leq r_{i}+r_{j}$ and $e_{i}j=0$ otherwise. Given a graph, the
package provides efficient functionality for computing the number of connected
components.
holes |  From a binary image, the number of holes is calculated by performing the same NCC operation on unoccupied cells and subtracting the number of CCs on the boundary.
To compute the number of holes given only disk locations and radii, we use
tools from the _TDA_ (topological data analysis) R package.
grid areas |  Computing local areas on subsets of the domain enables us to investigate how area is distributed or correlated within the domain by computing metrics of these sub-domains.
For this work we use a regular grid of sub-domains of size $1m^{2}$.  •
Moran’s I and Geary’s C — measures of spatial auto-correlation based on
adjacency; several versions of each statistic may be obtained by changing the
size of the sub-area and by changing how adjacency is defined (e.g., Moore vs
Von Neumann neighborhood, neighborhood range, weights, etc.). • sum of sub-
domain areas — this is an approximation to total area • number of full sub-
domains • number of empty sub-domains • sample variance of sub-domain specific
areas
empirical parameter estimates |  Empirical estimates of some model parameters are available from disk-based data.
With the number of observed fuel elements $n$, and the vector of observed disk
radii $\boldsymbol{r}$:  • $\hat{\lambda}=n/(d_{X}*d_{Y})$. •
$\hat{\mu}=mean(\boldsymbol{r})$ • $\hat{\sigma}=sd(\boldsymbol{r})$.
Incorporating these estimates as metrics encourages parameters consistent with
those estimated from the data.
Table 1: Metrics used to compare observed and generated binary fuels patterns

###  3.3 Bayesian Model Calibration

The calibration process is represented graphically in Fig 1.

Figure 1: Model calibration concept. a) Simulations from a generative model
across a prior range of control parameters $\theta$ are compared to
observations, resulting in refined knowledge of $\theta$ distributions, and
then b) the posterior distributions for $\theta$ lead to predictions with
appropriate uncertainty / variability.

Given a parametric generative model or forward simulation, model calibration
refers to the process of estimating model parameter settings that give the
most consistent output, or in other terms _fitting_ the model. Considering
uncertainty, _calibration_ of the generative model to data is determining a
posterior distribution of the model parameters given the observation data. New
predictions from the calibrated model are made with the inferred parameter
distributions for $\theta$.

We cast this in a standard Bayesian calibration framework (Santner et al.,
2018). A posterior distribution with respect to the metrics $y=\eta(\Xi)$ is
defined. The metric set $y^{obs}$ is extracted from observations of the
natural environment. The corresponding set of metrics $y^{gen}$ can be
extracted from a model generated example. These can then be compared through a
likelihood function. Here, we assume a likelihood function associated with a
multivariate normal distribution on the metrics, given the values of the
parameters $\theta$:

$\L(y^{obs}|\theta)=(2\pi)^{k/2}det(\Sigma_{y|\theta})^{-1/2}e^{-\frac{1}{2}\sum_{i}\sum_{j}(y^{obs}_{i}-y^{gen}_{j})^{T}\Sigma_{y|\theta}^{-1}(y^{obs}_{i}-y^{gen}_{j})},$
(5)

where $k$ is the length of metrics vector. The covariance matrix
$\Sigma_{y|\theta}$ is not known in practice and must be estimated, which is
discussed in Section 4. The next step in Bayesian calibration is to specify
priors on the parameters $\theta$. In this implementation the default prior on
lengthscale of the GP intensity function is Uniform, with bounds chosen based
on the domain size. The prior on the disk mean is $N(1.5,0.5^{2})$, truncated
to $[0,3]$, and the prior on the variance is $\Gamma(1,.001)$ which has a peak
near zero. These priors are reasonable defaults for many of the datasets we
have worked with, but may need to be adjusted for different ecosystems.
Together we denote the prior on the parameters as $\pi(\theta)$. The Bayesian
expression of the posterior distribution for the parameters for a given set of
data-derived metrics is:

$P(\theta|y^{obs})\propto L(y^{obs}|\theta)*\pi(\theta)$ (6)

The above expressions are un-normalized, presented as such since the
normalization constants are not necessary operationally for inference. In
principle, Eq. 6 can be optimized to find the maximum _a posteriori_ value of
the parameters, or it can be calibrated with Markov chain Monte-Carlo (MCMC)
to sample the posterior distribution (we will not rehearse the MCMC algorithm
here, it is easily available in references such as Gelman et al. (2013), which
also includes substantial advice on practice).

##  4 Model Calibration Process and Results

An initial step in verifying the soundness of the procedure proposed is to
infer the model based on so-called _perfect data_ \- data generated from the
model, but treated as observations. Through that process the capability can be
assessed quantitatively as well as qualitatively, albeit for the ideal case of
data known to be consistent with model assumptions. We then proceed to
outcomes of the same model inference procedure applied to real-world
observations. There several aspects to the statistical calibration process
with some procedural complexity, some of which are unique to this application.
This section will give some insight into practical considerations for
calibrating this generative model to data, and how prior domain information
can be introduced to aid inference.

###  4.1 Model Inference Process with Generated Data

We will first show the outcome of the calibration process, then fill in
additional underlying technical and procedural details. Recall that this
generative model has four parameters: the length scale of the Gaussian process
intensity function $\rho$, the point density per square meter $\lambda$, and
the mean $\mu$ and standard deviation $\sigma$ of the Normal distribution for
disk radius. Generated patterns using values of these parameters drawn from
parameter prior distributions is shown in the panels on the top row of Fig. 2.
There is clear qualitative difference in the character of the fuels layouts
generated with different parameter settings. To validate that the metric set
used allows the identification of these settings through comparison with
observation, the panels in the middle row of Fig. 2 show an example of model-
generated data as synthetic observation that will be used in model
calibration, where all five are generated using the same parameter set. The
lower row panels show fuels generated from samples of calibrated parameter
distributions, which are qualitatively similar to the training data, and we
will show the quantitative connection in model parameters.

Figure 2: Impact of parameters on generative binary layouts. Top row: five
parameter sets from their prior ranges demonstrates breadth of generative
behavior; middle: synthetically generated layouts from one setting of
parameters; bottom row: the constrained model behavior associated with the
calibrated model, showing 5 generated outcomes from parameter sets drawn from
the calibrated posterior distribution.

####  4.1.1 Stochastic Likelihood Evaluation

In this application, find settings of parameters
$\theta=\\{\rho,\lambda,\mu,\sigma\\}$ is complicated by the model output
being non-deterministic, where the generated spatial binary map is a
stochastic outcome for a given a setting of parameters. As defined previously,
we denote a generated binary map as $\Xi$, and the metrics generated from that
binary map as $y=\eta(\Xi)$. The observation is $y^{obs}$ (as well as the
model-generated pseudo observation where used), and a generated instance from
the model is $y^{gen}_{\theta}$. For a deterministic model, MCMC samples the
posterior distribution of parameters, summarized in Eq. 6, with the sole model
outcome ($J=1$) $y_{gen}$. For a stochastic model, each realization has a
different spatial outcome and corresponding metrics set, resulting in a sample
from a corresponding stochastic likelihood. The likelihood in Eq. 5 acts to
estimate an expectation of the likelihood over multiple realizations, although
still resulting in a stochastic estimate. The need for multiple realizations
from a Cox process for robust parameter estimation has been documented
(Micheas, 2019). This model has the same need for multiple realizations,
especially for estimating the intensity function length scale $\rho$.

There is a trade-off between the cost associated with the number of model
realizations used, and the algorithmic performance associated with the noisy
estimate of $E[L(y^{obs}|\theta)]$. For such estimates of the likelihood,
algorithmic methods such as MCMC can encounter optimistic tail values through
random variation, making it difficult efficiently sample the posterior. Our
approach to addressing this issue with MCMC is to re-calculate the current
value for $L(y^{obs}|\theta)$ for each proposed MCMC step. In practice, we
have found that, over a range of test problems with different parameter
settings, 25 realizations are sufficient to effectively constrain the
posterior while keeping the computational cost as low as possible. Some
parameter values may require more samples due to larger variability in
$y^{gen}$; this will ultimately depend on the details of the data and domain.

Preliminary studies in optimization have shown a greater intractability for
optimizers to converge with stochastic model estimates as compared to MCMC.
For this reason, optimization is not explored as an inference tool in this
work.

####  4.1.2 Covariance Estimation

In Eq. 5 we defined a normal likelihood on the $m$ observed plots
$y^{obs}_{i}\;i=1,...,m$, and explicitly denote that the covariance matrix is
dependent on $\theta$ because of the stochastic nature of the data generating
process for some $\theta$. As $\Sigma_{y|\theta}$ is unknown, an obvious plug-
in estimate is the empirical covariance matrix of $y^{obs}_{i};\;i=1,...,m$.
This works well when $m$ is reasonably large. We have found that $m=25$ is
sufficient for estimation of domains consistent with Fig. 2, resulting in
posterior uncertainty is significantly reduced from prior uncertainty for all
parameters and posterior modes are consistent with their true settings. For
smaller observed datasets of around $m\leq 10$ we find that, all else equal,
posterior uncertainty can be quite large and posterior modes can miss the true
parameter settings. A representative example is shown in Fig. 3. Posterior
distributions are much more sharply peaked with $m=25$ observations, and for
all four parameters posterior modes are near the truth. For $m=10$ however,
posterior uncertainty is significantly greater and inference for $\sigma$ and
$\lambda$ may be biased towards smaller estimates.

Figure 3: Model prior distributions (dotted) and calibrated posterior
distributions for the parameters for the intensity function length scale, disk
radius mean and variance, and overall intensity scaling average density. The
parameters values used to generate the synthetic observations are shown in red
vertical lines. Two posteriors are shown to demonstrate the value of observed
data.

There are two factors at play when comparing posterior estimates from the
$m=10$ and $m=25$ case. First is that the sum likelihood has more information
for larger $m$. Second is that the estimate of $\Sigma_{y|\theta}$ is improved
with larger $m$. The former is not to be dismissed, however we have found that
the latter may be the more significant effect. For our representative example
in Figs. (2, 3) we found that inference with $m=10$ using the covariance
matrix estimated with $m=25$ is nearly identical to inference using $m=25$ for
both estimation of the covariance and calibration. Improving covariance
estimation can significantly improve inference without collecting more data.

To improve covariance estimation for these smaller data applications we have
observed benefit in an approach of augmenting the observed dataset $y_{obs}$
with simulated data generated from random samples of $\theta$ from informed
prior distributions. The augmented dataset is used only for covariance
estimation, rather than likelihood calculations. Informed prior distributions
can be developed naturally from the available observations, with empirical
parameter estimates
$\hat{\lambda}=n/(d_{X}*d_{Y}),\hat{\mu}=mean(r),\hat{\sigma}=sd(r)$.

To augment the observed data: draw $m^{*}$ samples
$\boldsymbol{\lambda}_{s}\sim N(\hat{\lambda},.1\hat{\lambda})$,
$\boldsymbol{\mu}_{s}\sim N(\hat{\mu},.1\hat{\mu})$,
$\boldsymbol{\sigma}_{s}\sim N(\hat{\sigma},.1\hat{\sigma})$. The $10\%$
standard error is designed to account for the fact that our empirical
estimates are subject to the random data $y^{obs}$. We cannot define a simple
empirical estimate for $\rho$ so we sample $\boldsymbol{\rho}_{s}$ from weakly
informative distribution such as a Uniform or Gamma with large variance. For
each of these $m^{*}$ samples
$\boldsymbol{\theta}_{s}=[\boldsymbol{\lambda}_{s}|\boldsymbol{\mu}_{s}|\boldsymbol{\sigma}_{s}|\boldsymbol{\rho}_{s}]$
generate $K$ sets of random metrics
$y_{gen}^{k}(\boldsymbol{\theta}_{m});\;k=1,...,K$. The covariance matrix is
then estimated as the empirical covariance over these $m^{*}\times K$ sampled
metrics together with the $m$ observations.

Fig. 4 shows $95\%$ posterior confidence intervals over three values of
$m\in\\{1,10,25\\}$ for different covariance matrix estimation approaches. In
orange we have posterior intervals where only the $m$ available observations
are used for covariance estimation. Blue and green intervals have incorporated
simulated data in the estimation procedure. Both of these intervals use
samples from the normal distributions defined above for $\lambda,\mu,\sigma$
but vary the prior for $\rho$. We find that sensitivity to the prior on $\rho$
is relatively small, with similar results from a Uniform distribution on
$[0,10]$ and a more informative $\Gamma(3,1)$ distribution which has a mean at
the true value. For this example, the conditional likelihood for $\rho$ with
other parameters fixed at the truth is fairly flat.

This comparison demonstrates the connection between observation data
availability and the ability to identify and constrain the parameter posterior
distributions. It also shows that information regarding domain details, i.e.
by a domain expert, can result in good outcomes even when relatively data
limited.

Figure 4: Posterior $95\%$ intervals over different observed dataset sizes and
different priors for simulated data used to improve covariance estimation.
Orange intervals are generated using only the data for covariance estimation.
Green and Blue intervals use the observed data and simulated data to estimate
the covariance. Augmenting observations with simulations for covariance
estimation can significantly reduce posterior uncertainty when data is
limited.

###  4.2 Calibration with Field Data

Point patterns in spatial domains may arise either from direct observation or
through post-processing of technical measurement. Direct observation is
recording the locations of events of interest, often aided by GPS. When
feasible, this form of survey produces high-fidelity measurements. However, it
requires considerably more effort than either plot-level summaries or
remotely-sensed surveys. A LiDAR scan produces a point cloud (i.e., a set of
points) in 3D, but the points do not have the context of higher-level features
such as plants, and the point clouds are processed to identify features of
interest. For mid-story fuels, features of interest can be thought of as units
of fuel organization, which, at different scales, could be bunch grass,
bushes, or trees, or, more generally thought of as regions of high fuel
density in the environment. The locations of the these features form a point
pattern.

There are several possible pathways for extracting a heterogeneous spatial
layout of fuels from observations of the environment. Monoculture ecosystems
of chaparral or bunch grasses might be easily extracted from remotely sensed
imagery or ALS as a coverage map that can be used directly as a binarized
domain. A similar analysis is the extraction of the location and size of trees
from imagery and ALS, for example (Stereńczak et al., 2020). However, the
extraction of mid-story fuels is in general more difficult in an environment
including an overstory. This section will show a simple method for extracting
fuels elements from TLS scans, and used for calibration of the fuels model.

####  4.2.1 Terrestrial LiDAR Survey (TLS)

TLS generates a pointcloud of the scanned forest environment. The data
collected here uses the Leica BLK360 instrument, deployed to sites in the
Santa Fe National Forest. The demonstration site is ponderosa forest, which
includes aspen at various maturity, gamble oak, and other typical mid-story
flora. In this demonstration, the generated pointcloud is clipped to a 15m
square domain centered at the scanner. In order to extract mid-story from the
example scan, the ground is normalized to a ground plane, and the domain is
segmented to the region from 0.1m to 3m. Fig. 5 is a visualization of the
resulting mid-story pointcloud data.

Figure 5: mid-story pointcloud used for example data extraction. This shows a
ponderosa forest lidar scan, normalized for ground level, and segmented to
exclude ground and canopy.

This report focuses on methodology detail rather than the environmental
conditions or data collection methodology, and so we present this example
dataset as representative data without defending its environmental properties.

####  4.2.2 Data representation

We follow a simple reduction of the pointcloud data to spatially compact fuel
units, which is the concept that motivated the spatial model representation.
Since we are working in 2D for this demonstration, the z value of the data is
disregarded, leaving the marginal projection to the ground plane. A Gaussian
mixture model (GMM) of the 2D data is built using the python scikit library’s
BayesianGaussianMixture procedure, which fits a Dirichlet process model. The
result a representation of the data as a number of Gaussian distribution, in
this case constrained to be circular. A prior on the GMM fit controls the fuel
data inferred to be in a reasonable range of radii. The inferred 2D fuels
density is then discretized into a binary model of disks of 2 fitted GMM
standard deviation radius centered at the corresponding GMM means. The fitted
disks are shown in Fig. 6. Although this is an incomplete representation of
the TLS mid-story information and there are alternative approaches for
automating inventory, for example (Loudermilk et al. (2023)), and creating a
characteristic representation, this description does reflect the variability
of the major fuels elements in the domain, and in a format consistent with
characterization of the mid-story as distinct units of fuel in the conceptual
framework of the generative model presented. Specific of the algorithm will be
highly dependent on the particulars of an application context of fuels
generation; this background is a brief overview to introduce a path to
generating the data shown in Fig. 6.

Figure 6: Pointcloud data projected down to the x-y plane, and the Gaussian
mixture model fit to that data used to create a binary dataset. For
presentation clarity, color distinguishes distributions and point membership.

####  4.2.3 Model Calibration and Fuels Generation

Calibrating the generative model with this disk-based BMM data corresponds
directly to the calibration done with generated data. Calibration to the data
in Fig. 6 results in the posterior distributions shown in Fig. 7. What we can
see from these views of the 4-parameter joint distribution is that there is a
significant degree of learning parameter posteriors from the presented example
dataset, particularly for $\mu,\sigma,\lambda$. Covariance estimation for this
calibration example was done using the data augmentation approach outlined in
Section 4.1.2 with $K=25$, $m^{*}=25$, and $\boldsymbol{\rho}_{s}$ drawn from
a Uniform distribution on $[0,10]$.

Figure 7: Parameter posterior distributions from model calibration against GMM
representation of TLS mid-story. Empirical estimates from the single
observation of $\mu,\sigma^{2},\lambda$ are shown in red vertical lines.

There is no ground truth in parameters in this case of data extracted from
real TLS. Instead the evaluation is the qualitative performance of the fuels
layout. This is shown in Fig. 8 which, as before, shows some examples of the
prior behavior of the model before calibration, the data used in calibration,
and examples of the generative output from the calibrated model. Qualitative
examination shows that the layouts of the calibrated model are reasonably
representative, particularly in the goal of capturing the heterogeneity of the
binarized maps.

Figure 8: In the top row, five generated fuel layouts using parameters drawn
from their prior distributions. The observed calibration data is shown in the
middle row. In the bottom row, five generated fuel layouts from the calibrated
model posterior distribution.

##  5 A Framework for Multiple Spatial Influences

The model developed can be used along with other spatial covariates to provide
a more contextually relevant layout. This can be accomplished by considering
the generative model as only one component of a more general intensity
function. Revisiting Eq. 1, we can expand the intensity function heterogeneity
component $W(s)$ to include other controlling terms, as:

$\omega(s)=f\big{(}\beta_{0}W(s)+\sum_{k}\beta_{k}X(s)_{k}\big{)}$ (7)

That is, additional terms are added, where $X$ are spatial covariate fields,
and $\beta$ are weights to control the relative impact of spatial field
components. We will briefly show a generative example where $X_{1}$ is a
canopy height model (CHM) over a domain, which we will take as an indicator of
greater likelihood of finding a mid-story fuel element, i.e. a shrub. The
canopy height model used for this example is shown in Fig. 9(a). In addition
to a canopy height covariate, we also incorporate an $X_{2}$ of the location
of a road, where fuel probability is zero. These are shown in Fig. 9(b).

One important consideration in this model is the need for a convention for the
semantics of the $X$ datasets. We will specify $X_{k}\in[-1,1]$, interpreted
as proportional to a likelihood of fuel. $\beta_{k}$ is a positive multiplier
interpreted as a relative scaling of the impact of indicator $k$. We will not
investigate the ability to infer $\beta_{k}$ parameters, and instead will take
them as supplied by a user guiding the system. To estimate $\beta$, data
requirements would be substantially increased, and likely the need for expert-
supplied guidance more critical, as well as complex. It is important to
recognize that while the relative scaling of the $\beta_{k}$ controls the
impact of indicators $X$ in the generated outcome, at the same time the
overall sum of $\beta_{k}$ will impact the layout as the total is projected
through the logistic function.

Fig. 9(c) shows an example of placing fuel element centers in the domain of
the CHM with the heterogeneity model only. In Fig. 9(d), canopy heights are
used to bias fuel placement toward areas of higher canopy. By setting $X_{2}$
to $-1$ at road locations and $0$ elsewhere, an effectively zero fuel
probability can be given to those locations by setting $\beta_{2}$ large.

This section is a demonstration example of a topic requiring further
development through field experience with these tools. It seems clear that a
number of (potentially) available spatial covariates may be important in fuels
layout, such as slope, aspect, canopy density, etc. We have seen some
preliminary results that the $\beta_{k}$ can be identified in some
circumstances, but in general increasing the problem complexity will make this
difficult, likely greatly increasing the need for observation data and/or
expert influence. Instead, our recommendation would be to make inferences of
the heterogeneity over controlled observations, and carry those settings
forward for expert tuning of the data.

(a) Canopy Height Model (CHM)

(b) Location of road.

(c) Generation without CHM influence.

(d) Generation with CHM influence.

Figure 9: In panel (a) a canopy height model is shown, scaled to [-1,1]. Panel
(b) shows the locations of the road, which should be set to zero probability
of fuel generation. Panel (c) shows a realization from the data generating
process given parameters $\tilde{\theta}$ with $\beta=(1,0,0)$ indicating no
effect from the canopy. Panel (d) is a realization using the same
$\tilde{\theta}$ but $\beta=(1,1,10)$ gives a strong correlation between
canopy height and fuel placement. $\beta_{2}=10$ gives effectively zero
probability of fuel placement at road locations indicated by $X_{2}$.

##  6 Discussion and Conclusion

This study shows the potential of a connection from LiDAR remote sensing of
mid-story to generating representative examples of the mid-story that manifest
key characteristics. This is accomplished through the definition of a
parameterized generative model that can capture these characteristics of the
layout, and the ability to calibrate, or learn, the parameters of this model.
The result is an interpretable method that can be used to compactly summarize
the relevant information from collected datasets, and populate unknown areas
with characteristic patterns of fuels. Although we emphasize the potential for
generative model calibration, in a practical application it is reasonable to
also include knowledge of the environment as either a more constrained prior
parameter distribution.

A key issue is the definition of the relevant features of fuels
characteristics. This work is motivated by prior investigations into
environmental generative models, and reinforced by the information considered
useful as embodied in historical plot survey approach to summarize the spatial
heterogeneity of fuels. Further, our approach is a reasonable method for
generating heterogeneous fuels layouts for 3D fire simulations by interpreting
the layout of disks as locations and scaling factors for a template 3D fuel
example, i.e., a shrub. The model presented is consistent with those
motivations.

The model presented should be regarded as an initial method investigated for
both its flexibility and its ability to be calibrated against both perfect
data and data collected in the field. There are clear extensions of this
model, as well as the method of connecting the model to observations (i.e.,
pointcloud processing, matching metrics, and likelihood evaluation) that would
be interesting, although their properties must be the subject for further
investigations. Interesting extensions could include fuel type (e.g., species
or coarse woody debris type) and a description of 3D properties including
shape and density. In addition, further exploration of spatial covariates for
fuel placement is appropriate for enhancing the ability to capture
environmental conditions related to, e.g., terrain, overstory, and
moisture/water. Finally, although we focus on mid-story, this generative
approach can be applied to fuels from surface to canopy. However, although the
model can be extended in many ways, one of the key findings is that model
complexity must be controlled to avoid models that are not feasible to infer,
to specify with respect to known environment principles, and/or to interpret.

The method of describing shape or topology, and the connection to kernels in
the stochastic intensity function appears to be an open area for
investigation. In application, there is limited guidance on what properties of
heterogeneity and shape are relevant to environmental assessment goals,
including the functional impact to wildland fire outcome. This suggests
several research directions both for the underlying statistical methodology,
and the connection to application pragmatic outcomes.

## Acknowledgements

This research was supported by the Los Alamos National Laboratory (LANL)
through its Laboratory Directed Research and Development (LDRD) program under
project number 20220024DR.

## References

  * Anderson et al. (2021) C. Anderson, S. Dietz, S. Pokswinski, A. Jenkins, M. Kaeser, J. Hiers, and B. Pelc.  Traditional field metrics and terrestrial lidar predict plant richness in southern pine forests.  _Forest Ecology and Management_ , 491:119118, 2021.  ISSN 0378-1127.  doi: https://doi.org/10.1016/j.foreco.2021.119118.  URL https://www.sciencedirect.com/science/article/pii/S0378112721002073.
  * Andrews (2018) P. L. Andrews.  The rothermel surface fire spread model and associated developments: A comprehensive explanation.  _Gen. Tech. Rep. RMRS-GTR-371. Fort Collins, CO: US Department of Agriculture, Forest Service, Rocky Mountain Research Station. 121 p._ , 371, 2018.
  * Atchley et al. (2021) A. L. Atchley, R. Linn, A. Jonko, C. Hoffman, J. D. Hyman, F. Pimont, C. Sieg, and R. S. Middleton.  Effects of fuel spatial distribution on wildland fire behaviour.  _International journal of wildland fire_ , 30(3):179–189, 2021.
  * Banerjee et al. (2003) S. Banerjee, B. P. Carlin, and A. E. Gelfand.  _Hierarchical modeling and analysis for spatial data_.  Chapman and Hall/CRC, 2003.
  * Diggle (1981) P. J. Diggle.  Binary mosaics and the spatial pattern of heather.  _Biometrics_ , pages 531–539, 1981.
  * Gelman et al. (2013) A. Gelman, J. B. Carlin, H. S. Stern, D. Dunson, A. Vehtari, and D. Rubin.  _Bayesian Data Analysis_.  Chapman and Hall/CRC, 2013.
  * Hiers et al. (2009) J. K. Hiers, J. J. O’Brien, R. J. Mitchell, J. M. Grego, and E. L. Loudermilk.  The wildland fuel cell concept: an approach to characterize fine-scale variation in fuels and fire in frequently burned longleaf pine forests.  _International Journal of Wildland Fire_ , 18:315–325, 2009\.
  * Jonko et al. (2021) A. K. Jonko, K. M. Yedinak, J. L. Conley, and R. R. Linn.  Sensitivity of grass fires burning in marginal conditions to atmospheric turbulence.  _Journal of Geophysical Research: Atmospheres_ , 126(13):e2020JD033384, 2021.  doi: https://doi.org/10.1029/2020JD033384.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020JD033384.  e2020JD033384 2020JD033384.
  * Knapp and Keeley (2006) E. E. Knapp and J. E. Keeley.  Heterogeneity in fire severity within early season and late season prescribed burns in a mixed-conifer forest.  _International Journal of Wildland Fire_ , 15(1):37–45, 2006.
  * Li et al. (2012) W. Li, Q. Guo, M. K. Jakubowski, and M. Kelly.  A new method for segmenting individual trees from the lidar point cloud.  _Photogrammetric Engineering & Remote Sensing_, 78(1):75–84, 2012.
  * Linn et al. (2002) R. Linn, J. Reisner, J. J. Colman, and J. Winterkamp.  Studying wildfire behavior using firetec.  _International journal of wildland fire_ , 11(4):233–246, 2002.
  * Linn et al. (2005) R. Linn, J. Wintercamp, J. J. Colman, C. Edminster, and J. D. Bailey.  Modeling interactions between fire and atmosphere in discrete element fuel beds.  _International journal of wildland fire_ , 14(1):37–38, 2005.
  * Linn et al. (2012) R. Linn, K. Anderson, J. Winterkamp, A. Brooks, M. Wotton, J.-L. Dupuy, F. Pimont, and C. Edminster.  Incorporating field wind data into firetec simulations of the international crown fire modeling experiment (icfme): preliminary lessons learned.  _Canadian Journal of Forest Research_ , 42(5):879–898, 2012.
  * Linn et al. (2020) R. Linn, S. Goodrick, S. Brambilla, M. Brown, R. Middleton, J. O’Brien, and J. Hiers.  Quic-fire: A fast-running simulation tool for prescribed fire planning.  _Environmental Modelling & Software_, 125:104616, 2020\.
  * Linn et al. (2013) R. R. Linn, C. H. Sieg, C. M. Hoffman, J. L. Winterkamp, and J. D. McMillin.  Modeling wind fields and fire propagation following bark beetle outbreaks in spatially-heterogeneous pinyon-juniper woodland fuel complexes.  _Agricultural and Forest Meteorology_ , 173:139–153, 2013\.
  * Linn et al. (2021) R. R. Linn, J. L. Winterkamp, J. H. Furman, B. Williams, J. K. Hiers, A. Jonko, J. J. O’Brien, K. M. Yedinak, and S. Goodrick.  Modeling low intensity fires: Lessons learned from 2012 rxcadre.  _Atmosphere_ , 12(2), 2021.
  * Loudermilk et al. (2023) E. L. Loudermilk, S. Pokswinski, C. M. Hawley, A. Maxwell, M. R. Gallagher, N. S. Skowronski, A. T. Hudak, C. Hoffman, and J. K. Hiers.  Terrestrial laser scan metrics predict surface vegetation biomass and consumption in a frequently burned southeastern u.s. ecosystem.  _Fire_ , 6(4), 2023.  ISSN 2571-6255.  doi: 10.3390/fire6040151.  URL https://www.mdpi.com/2571-6255/6/4/151.
  * McDanold et al. (2023) J. S. McDanold, R. R. Linn, A. K. Jonko, A. L. Atchley, S. L. Goodrick, J. K. Hiers, C. M. Hoffman, E. L. Loudermilk, J. J. O’Brien, R. A. Parsons, C. H. Sieg, and J. A. Oliveto.  Duet - distribution of understory using elliptical transport: A mechanistic model of leaf litter and herbaceous spatial distribution based on tree canopy structure.  _Ecological Modelling_ , 483:110425, 2023.  ISSN 0304-3800.  doi: https://doi.org/10.1016/j.ecolmodel.2023.110425.  URL https://www.sciencedirect.com/science/article/pii/S0304380023001564.
  * Micheas (2019) A. C. Micheas.  Cox point processes: Why one realisation is not enough.  _International Statistical Review_ , 87(2):306–325, 2019.  doi: https://doi.org/10.1111/insr.12308.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/insr.12308.
  * Molchanov (1997) I. Molchanov.  _Statistics of the Boolean model for practitioners and mathematicians_.  Wiley, 1997.
  * Møller and Helisová (2010) J. Møller and K. Helisová.  Likelihood inference for unions of interacting discs.  _Scandinavian Journal of Statistics_ , 37(3):365–381, 2010.
  * Møller and Waagepetersen (2007) J. Møller and R. P. Waagepetersen.  Modern statistics for spatial point processes.  _Scandinavian Journal of Statistics_ , 34(4):643–684, 2007.
  * Parsons et al. (2017) R. A. Parsons, R. R. Linn, F. Pimont, C. Hoffman, J. Sauer, J. Winterkamp, C. H. Sieg, and W. M. Jolly.  Numerical investigation of aggregated fuel spatial pattern impacts on fire behavior.  _Land_ , 6(2):43, 2017.
  * Pielou (1964) E. Pielou.  The spatial pattern of two-phase patchworks of vegetation.  _Biometrics_ , pages 156–167, 1964.
  * Pokswinski et al. (2021) S. Pokswinski, M. R. Gallagher, N. S. Skowronski, E. L. Loudermilk, C. Hawley, D. Wallace, A. Everland, J. Wallace, and J. K. Hiers.  A simplified and affordable approach to forest monitoring using single terrestrial laser scans and transect sampling.  _MethodsX_ , 8:101484, 2021.  ISSN 2215-0161.  doi: https://doi.org/10.1016/j.mex.2021.101484.  URL https://www.sciencedirect.com/science/article/pii/S2215016121002776.
  * Riley et al. (2016) K. L. Riley, I. C. Grenfell, and M. A. Finney.  Mapping forest vegetation for the western united states using modified random forests imputation of fia forest plots.  _Ecosphere_ , 7(10):e01472, 2016.
  * Rothermel (1972) R. C. Rothermel.  _A mathematical model for predicting fire spread in wildland fuels_ , volume 115.  Intermountain Forest & Range Experiment Station, Forest Service, US …, 1972.
  * Rowell et al. (2020) E. Rowell, E. L. Loudermilk, C. Hawley, S. Pokswinski, C. Seielstad, L. Queen, J. J. O’Brien, A. T. Hudak, S. Goodrick, and J. K. Hiers.  Coupling terrestrial laser scanning with 3d fuel biomass sampling for advancing wildland fuels characterization.  _Forest Ecology and Management_ , 462:117945, 2020.
  * Santner et al. (2018) T. J. Santner, B. J. Williams, and W. I. Notz.  _The design and analysis of computer experiments, second ed._ Springer, 2018.
  * Silva et al. (2016) C. A. Silva, A. T. Hudak, L. A. Vierling, E. L. Loudermilk, J. J. O’Brien, J. K. Hiers, S. B. Jack, C. Gonzalez-Benecke, H. Lee, M. J. Falkowski, and A. Khosravipour.  Imputation of individual longleaf pine (pinus palustris mill.) tree attributes from field and lidar data.  _Canadian Journal of Remote Sensing_ , 42(5):554–573, 2016.  doi: 10.1080/07038992.2016.1196582.
  * Stereńczak et al. (2020) K. Stereńczak, B. Kraszewski, M. Mielcarek, Żaneta Piasecka, M. Lisiewicz, and M. Heurich.  Mapping individual trees with airborne laser scanning data in an european lowland forest using a self-calibration algorithm.  _International Journal of Applied Earth Observation and Geoinformation_ , 93:102191, 2020.  ISSN 1569-8432.  doi: https://doi.org/10.1016/j.jag.2020.102191.  URL https://www.sciencedirect.com/science/article/pii/S0303243420301665.
  * Tinkham et al. (2018) W. T. Tinkham, P. R. Mahoney, A. T. Hudak, G. M. Domke, M. J. Falkowski, C. W. Woodall, and A. M. Smith.  Applications of the united states forest inventory and analysis dataset: A review and future directions.  _Canadian Journal of Forest Research_ , 48(11):1251–1268, 2018.
  * Toney et al. (2007) C. Toney, M. Rollins, K. Short, T. Frescino, R. Tymcio, and B. Peterson.  Use of fia plot data in the landfire project.  In _In: McRoberts, Ronald E.; Reams, Gregory A.; Van Deusen, Paul C.; McWilliams, William H., eds. Proceedings of the seventh annual forest inventory and analysis symposium; October 3-6, 2005; Portland, ME. Gen. Tech. Rep. WO-77. Washington, DC: US Department of Agriculture, Forest Service: 309-319._ , volume 77, 2007.