We’re going to look at some beach morphology data from Benson Beach, just north of the North Jetty on the Columbia River. The beach is accessible from Cape Disappointment State Park near Ilwaco, WA.

In 2010 and 2011, the USGS, OSU, WA State Department of Ecology, and the Army Corps of Engineers did a joint study of the impacts of beach nourishment. The Army Corps traditionally dumps dredging materials offshore, outside of the reach of littoral sediment transport. In the summer of 2010, they deposited these materials on high-energy, often erosive Benson Beach.

Today, we’ll do a difference analysis of beach elevations before and after the beach nourishment. Data was taken with foot- and Polaris side-by-side-borne GPS, and with jetskis equipped with GPS and sonar.

 

Start by copying the file S:\Geo599\Data\SEXTON_Difference_Tutorial into your own folder.

We’re going to start with the point data. It’s been cleaned, but otherwise has not been processed. I had some trouble with opening it, so please let me know if you have any issues – I’ve got a back-up plan.

The S01 data is the “before nourishment” baseline survey, performed on July 11, 2010. The S06 data is the “after nourishment” survey, performed September 22, 2010.

 

Open ArcMap.

In the catalog, navigate to where you saved the folder. Find the S01 folder and right click on the file: S01_20100711_topo.

Create a feature class from XY table.

Choose Field 1 for X, Field 2 for Y

Leave Z as “none” (we have Z data, but it got real mad when I told it here).

Click the Coordinate System of Input Coordinates button:

Select “Projected Coordinate Systems,” then “State Plane,” then “NAD 1983 (Meters).”

Navigate to and select the coordinate system called:

NAD_1983_StatePlane_Washington_South_FIPS_4602

Choose where to save your output feature class. I called mine: XYS01_20100711_topo.shp

Hit OK and cross your fingers – this is where your program might crash!

Once it’s done thinking (and maybe not responding for a bit – don’t panic if that happens, it kept working for me), navigate to the new shape file in your folder and drag it into your working area.

This should be pretty exciting, except that all the symbols are one color. To better visualize, right click on the layer, and find properties. Click the “Symbology” tab and then “Quantities” on the left. Change the “Value” to “Field3” and select a color ramp. Field three is your Z values, and it will color by elevation!

PS: when you tell it to use Field3, it gives you an angry message. I just hit OK, and it moved on and looked fine. If anyone has any suggestions, I’m all ears.

 

While this was super cool, I thought we’d short cut pulling in the rest of the XY data, and I’ve provided you with shape files of the rest.

In the S01 folder, there is a help folder. In there, you’ll find my version of the S01 topo shape file (which you likely just completed) and aS01 bathy shape file, called XYS01_20100711_bathy. Drag that into the work space. You can change the colors to see the elevation, but it’s not necessary to our next steps.

 

You’ll also need to drag in the S06 shape files. I still gave you the XYZ text files in the S06 folder, but again provided a “help” folder with shape files of each point set. Drag these files into your work space:

XYS06_20100922_topo.shp

XYS06_20100922_bathy.shp

You should now have four layers: two bathy, two topo; one from each survey.

Now we’re going to make some surfaces.

In the toolbox, open 3DAnalyst, then Data Management, then TIN, and find the “Create TIN” tool. (or search for it…)

Name your first TIN S01

You’ll need to choose the projection again: NAD_1983_StatePlane_Washington_South_FIPS_4602

Then choose the two S01 feature classes, bathy and topo.

When they appear in the list of input features, you’ll need to change the Height Field to Field3!

Hit OK and it will make a TIN.

Hopefully, it looks like this:

Repeat the process for the S06 data set. Don’t forget to specify the Height is in Field3.

Lastly, we’re going to perform a difference analysis. This will compare the baseline beach survey with the “after nourishment” beach survey and find where sand was deposited, and where erosion happened after the deposition (since the survey wasn’t directly after).

 

In the 3DAnalyst toolbox, expand the “Triangulated Surface” options, then open the “Surface Difference” tool.

Set the Input Surface to “S06”

Set Reference Surface to “S01”

Name your output feature, I called mine “Beach_Change”

Expand the Raster Options and name an output raster, I also called this one “beach_change”

Leave cell size as default (10)

Click OK! This might take a minute.

The resulting raster can be recolored to more readily show differences.

The feature class generated currently shows only where there sand was deposited (blue), eroded (green), or unchanged.  While this might be helpful, I wanted elevation changes in a feature class.

If you’d rather a shape file:

Open the attribute table for surface_dif1 feature class

Add a field, call it “ele_dif” and make it a “float”

Open field calculator and calculate the elevation change by multiplying volume and the code (negative means loss of sediment, positive means gain), and dividing by “SArea.” Then click okay.

Close the attribute table and open Properties. Choose the “Symbology” tab, then “Quantities.”

Choose the new “ele_dif” as the value, and increase the number of classes (I chose 9). When you hit okay, the feature class will show the elevation changes!

 

Here is a simple tutorial in ArcGIS on using spatial statistics tools. Not much new here but I found it useful and interesting because it uses the tools to help inform analysis with other tools in the toolset.

• Go to blog website:
  • Follow linkà http://www.arcgis.com/home/item.html?id=020e1d9e9d2d44f9ae7da8663cf3bac5
  • Googleà spatial statistics blog arcgis->Tutorial/Training->Spatial Pattern Analysis in ArcGIS 10
• Download tutorial data: Save to your student folder
• Follow along!

Here are the web links and full reference of the Isaak et al. (2014) paper I showed while going through the spatial stream network model tutorial in class on Monday. I (and many others) have been waiting a long time for the availability of methods such as these and now we have them!

SSN/STARS – Spatial Statistics on Stream Networks

OLD ESRI Blog post (tools are now updated for ArcGIS 10.2)

NorWeST Regional Stream Temperature

Daniel J. Isaak, Erin E. Peterson, Jay M. Ver Hoef, Seth J. Wenger, Jeffrey A. Falke, Christian E. Torgersen, Colin Sowder, E. Ashley Steel, Marie-Josee Fortin, Chris E. Jordan, Aaron S. Ruesch, Nicholas Som and Pascal Monestiez. 2014. Applications of spatial statistical network models to stream data. Wiley Interdisciplinary Reviews: Water. DOI: 10.1002/wat2.1023

I would like to show you how to use a script tool within the ArcGIS toolbox to call R remotely, have it perform some function on your data, and return the results into a new shapefile projected onto your workspace in ArcMap.  I won’t be doing that.

The scripting tool works.  I have used it myself.  But it does not work in the computer lab at school.  The problem has to do with how R handles packages.  When installing a package at school, R complains that students do not have permission to modify its main library, and automatically places new packages into a user-specific library.  R does not have any trouble finding this user library when running on its own, but when calling R inside ArcMap it nonsensically forgets where the user-library is located on the S: drive.  Considering all spatial statistics in R requires user-installed packages like maptools, there is no point to calling R and asking it to do anything if it can’t find its own library.  So instead I am going to give a plain vanilla tutorial on how to open a shapefile in R, lay some statistics smackdown on it, and save the results in a new shapefile and some tables you can open in ArcMap.  You might want to grab a cup of coffee, this is going to be boring.

If, for some reason, you actually want to know how to call R remotely using a python script tool within ArcGIS, I will leave a trail of breadcrumbs for you at the end of the tutorial.  It is not that hard.

First, locate the folder “demo_geo599” within Geo599/Data/ and copy the folder to your student folder.

Open the program R, press Ctrl + O and navigate to your copy of the demo_geo599 folder.  Open the file “demo_nb_GLM.r”.

Install the packages required for this demo using the code at the top of the file.  Once the packages are installed, you still need to load them.  Note that you can load packages using either the command “library” or “require”, why I don’t know.

The next step to getting this demo working requires you to change a line of code.  Where the file says “Set Working Directory”, below it you will find the command “setwd()”.  Inside those braces you need to replace the file path with the path to the demo folder in your student folder on the S: drive.  R interprets the file path as a string, so be sure to surround the file path with single or double quotation marks, which is how you declare a string in R.  By running this command, you direct all output to this file path, and you can manipulate file names within this path without specifying the

The next line of code uses the package maptools to import the shapefile using the command “readShapePoints()”.  This line of code will run so long as your file path specifies the demo folder.  The shapefile “Local_Farms.shp” contains covariate data on some of the farms serving Corvallis Farmer’s Market, the ones that had responded to my survey at the time of making this demo.  Examine the file contents using the “head()” command, and note that the shapefile data resembles an R data.frame.  Although it includes a few additional properties, in many ways you can manipulate the shapefile data the same way you would a data.frame, which is how we will add the results of our regression to the data and output a new shapefile.

I adapted this file from the script that python calls, as I mentioned above, so this next bit of code is a legacy from adapting that file.  There are easier ways to specify a formula.  In this case, I specify the names of the variables as strings, and then paste them together into a formula.  Formulas in R follow a syntax:  y ~ x1 + x2 + x3.   Before we fit the model, we are going to perform some initial diagnostics.  The package GGally include a function called “ggpairs()”.  This function makes a scatterplot of every variable in a data.frame against every other variable and presents the results in the lower triangle of a matrix.  In the upper triangle it prints the correlation between the two variables.  If you run the command alone, these results will appear in the graphic window.  Above the line is the command you use to print the graph as a portable network graphic (.png), and the command “dev.off()” tells R that you are done specifying the contents of the .png file.  While you likely have no interest in keeping a pretty picture of the output, this format is essentially when you are calling R remotely from ArcGIS.  This is because the R script runs in full and exits within the python script and you can’t see any of the graphs.  So after I run the file remotely in ArcGIS, I open the working directory and examine the graphs.  From the perspective of diagnostics, it is better to work in R directly.  If you need to change the scale of the response variable or, you know, respond to the diagnostics in any way, you can’t do that when running the script remotely in ArcGIS.

To make sure that we only compare variables of interest, I preface the ggpairs() call with a few lines of code that build a data.frame containing only the variables we want to compare.  The line “#var_tab$Acreage <- as.integer(var_tab$Acreage)” begins with a #, which is how you specify comments in R.  That means this line is “commented out”.  We will be running a negative binomial regression, and R assumes you will provide count data as the dependent variable.  The regression still operates with float data, but R will issue warnings, which you can safely ignore.  If you do not want to see warnings, simply round off all the decimals by uncommenting this line and running it.

Note that I have commented out another line in the script: #ct <- qt(.95, df=(nrow(var_tab) – length(independentVars))).  I didn’t end up using this because the package I used to report coefficient results did not allow me to specify different confidence intervals using critical values, but this is too bad.  I do not have enough farms to justify a .95 probability and should be using a t distribution to determine confidence intervals.

You may be wondering: why negative binomial regression?  Well, “number of acres” is essentially count data.  I have no reason for assuming it to be normally distributed around any particular mean acreage.  In fact, I expect smaller farms to be more common than larger, so it is not unreasonable to hypothesize that the distribution of farms would match a poisson or negative binomial distribution.  Turns out, the poisson regression was over-dispersed, so the negative binomial is the appropriate distribution in this case.  It seems to me that the poisson is always over-dispersed, and I am getting into the habit of skipping straight to negative binomial.  For an interesting analysis of the suitability of different distributions for analyzing ecological data, I recommend a paper by Walter Stroup called “Rethinking the Analysis of Non-Normal Data in Plant and Soil Science”.

Go ahead and run the regression and examine the results using the command “summary()”.  Note that the spatial component of my analysis “miles”, which represents the distance of the farm from the Corvallis Farmer’s Market in miles, is not significant.  This is a shame because I hypothesized that farms would become smaller the closer they were to the market.  In particular, I believe that there is a new business model for farms emerging based on small acreages producing high quality food intended for direct-to-consumer local markets.  These results show that older farms tend to have larger acreage, and the higher proportion of their production that goes to local markets, the smaller the acreage tends to be.

The package texreg is for formatting model results into the LaTeX language, which many academics use for publishing because it specializes in formatting model formulas and symbols.  The function “plotreg()” makes an attractive graphic of the coefficient results.  The next few lines of code simply print a Q-Q plot and a plot of the residual vs. fitted values, so we can test the model for normality and independence.  The errors are reasonably normally distributed.  You can see in the residuals vs. fitted plot that most of the farms are very small in size, but the errors seem independent of the farm size.  Hard to tell with so few data points!

Next we append the residuals and how many standard errors each residual deviates from the mean to the shapefile data.  This is the same syntax used to add data to a data.frame.  The function writeSpatialShape() resides in the maptools package and does exactly what it sounds like, depositing the new shapefile in the working directory.  Then, since we cannot see the output from R when calling the script remotely in ArcGIS, the code prints the results in a .dbf file.  The function write.dbf() is in the package foreign.

Information criteria are all the rage now, so the final part of the code makes a .dbf table that includes every possible combination of independent variables included in the original model fit against the dependent variable, with the resulting AIC and BIC.  Unless you have an abiding fascination with for and if loops, you can safely ignore the algorithm I use to permutate through every unique combination of independent variables.  It stores all the combinations as vectors of numbers in a list.  The code after it connects the variable names to the numbers and then pastes them together into formulas.  This part of the code is only interesting if you work with list type data and want to know how to extract and manipulate data contained in lists.  Hint: use the command “unlist()”.

Use write.dbf again to produce a table you can load in ArcGIS, and feel free to run the line containing just the name of the table “diag_tab” to view the diagnostics table in R.

This concludes a poor tutorial in R.

If you are interested in how ArcGIS and R can talk to each other, here are a few tips to get you started.

I found this video by Mark Janikas describing a python script that calls R to perform a logit regression and return the results to a shapefile in ArcGIS.  The video is not super helpful in explaining how this process works, its intended audience is a conference full of pro programmers who are hip to hacking lingo, which is not me.  At the end, however, he refers the interested to this github site that hosts the python and r scripts.

One important note: certain ArcGIS directories have to be added to the search path for Python.  You can manually add these directories to the search path by modifying the path variable in the environmental variables according to the instructions in the folder demo_geo599/R_Tools/Doc.  However, I did this and it did not work, so I used a workaround.  I manually add the paths using the code in the python file “pathlib.py” which I include in the demo folder.  This file and the file “arcpyWithR.py” need to be in the python search path somewhere for my example file to work.  And yes it will work on your home machine, just not in the lab.  I have added a python file and R file titled “GLMwithR” to the folder R_Tools/Scripts.

If you open the python file, you will note the first line reads “import pathlib”, which adds the necessary directories to the search path.  Scrolling down, you will see the function “ARCPY.GetParameterAsText()” called six times.  Each of these commands receives user input through the ArcGIS GUI interface.  You can add the toolbox inside the R_Tools folder to your ArcGIS toolbox by right-clicking on the directory in the toolbox, selecting “Add Toolbox…” and navigating to the “R Tools” toolbox in the R_Tools folder.  I added my script to this toolbox by right-clicking on the box and selected “Add -> Script…”.  I would refer you to the help file on creating script tools.  The hardest part is actually configuring the parameters.  If you right-click on the “Negative Binomial GLM” tool and select “Properties”, you can view the parameters by selecting the “Parameter” tab.

Note that there are six parameters.  These parameters need to be specified in the same order as they appear on the python script, and each parameter needs to be configured properly.  In the Property field, you have to select the proper values for Type, MuliValue, Obtained from, etc., or the script will fail to work.  It’s pretty annoying.

I don’t understand all the python code myself, but the middle portion passes all the inputs to R in the list called “args”.  Feel free to open the r script “GLMwithR.r” in the folder R_Tools/Scripts and note the differences between it and the tutorial file we just worked through.

Python uses the functions “ARCPY.gp.GetParameterInfo()” and “ARCPY.SetParameterAsText()” to read the new shapefile and the two .dbf tables respectively that come from the R script.  The whole process of sending info in between Python and R is still mysterious to me, so I can’t really shed any light on how or why the python script works.

That concludes my trail of breadcrumbs.  Good luck!

 

 

The goal of this project was to use R to run statistical operations on data from ESRI shapefiles, and see about getting the results of these operations back into a format that can be used in ArcGIS.  This was my first experience using R, and I was surprised at how easy it was (or at least parts of it), though having experience with Python and Java helps.  This may be pretty basic for those who already have used R, but hopefully it will be a useful introduction for neophytes like myself.

An easy way to import ESRI shapefile data into R is through the use of the maptools package.  You can install the package in the R gui by selecting Packages>Install Package(s), selecting your mirror site (we have one at USA (OR)), and finding maptools in the list.  You’ll also need the rgl, spatstat, sm, and misc3d packages to run all the commands here.  Note that Digital Earth uses R 3.0.1, for which spatstat is not available, so some of these commands are not possible in the lab.

For this exercise, we’ll be using the Mount Hood earthquake data (hood3d.shp):

To read a shapefile, you can use one of the following commands in R, depending on the type of shapefile:

readShapeLines(file.choose())
readShapePoints(file.choose())
readShapePoly(file.choose())

The file.choose command is an alternative to hard-coding a file path and name; it will bring up a file selection dialog that you can use to interactively select your file.

Run the following command to create a SpatialPointsDataFrame object called shape, which will contain the geometry and attributes of the shapefile:

shape <- readShapePoints(file.choose())

You can view the data with the summary or attributes commands:

Note that the X, Y, and Z coordinates are named here as coords.x1, coords.x2, and coords.x3.  We can assign each to a more intuitively-named vector (or list for those who are more used to Python):
x <- shape$coords.x1
y <- shape$coords.x2
z <- shape$coords.x3

Now we can do a number of things with these points.  One of the more fun options is to compute the 3d density of these points and create density isosurfaces.  This uses the kde3d (kernel density in 3d) and contour3d functions in the misc3d package.

Our data is in a UTM projection, which means that the units are meters, which means that the distances between points are large numbers.  I haven’t yet been able to get good results using the raw data in the kde3d function, but do get good results when dividing the input values by 1000:

mx <- x/1000
my <- y/1000
mz <- z/1000

Through trial and error, I got decent results using the following parameters:

with(shape, {
d <- kde3d(mx, my, mz, h=0.2, n = 150, lims=c(range(mx),range(my),range(mz)))
contour3d(d$d, 0.02, d$x, d$y, d$z,
color = "red", color2 = "gray", scale=TRUE
)
})

The result displays in an interactive window that can be zoomed and rotated:

There is another package that depends on the misc3d package, called sm.  It has a density function that calls contour3d and configures some more advanced parameters behind the scenes.  It’s also somewhat simpler to use:

m <- cbind(x,y,z)
sm.density(m, panel=TRUE, panel.plot = TRUE)

Note that  the raw x, y, and z coordinates worked here, though it is necessary to create a matrix with them first.  The result:

Notice the slider for the size of kernel used in the density calculation, as well as the bounding frame and 3d axes.  There are a number of arguments that can be used in both the density and contour3d functions to customize the appearance of the data and the layout.  By either using them or by writing your own R code, you could conceivably make 3d density plots to whatever level of sophistication you wish.

In addition to visualizing your data, you can perform whatever statistical analysis you would use with any other dataset in R.  For example, we can easily make histograms of depth or magnitude:

layout(matrix(c(1,2),1,2,byrow=TRUE))
hist(shape$Depth, breaks=20)
hist(shape$Magnitude, breaks=40)

The layout function allows for the creation of two side-by-side plots when given a matrix of one row and two columns; otherwise the second plot would overwrite the first:

Histograms can just as easily be made in ArcGIS, but there are much more complex data operations available in R.  For example, it is possible to  approximate the above density surfaces and get them back to ArcScene with nearest neighbor distances.

First, we make a 3d points object using the x, y, and z values as well as their ranges.  In this case, this object is called p3.  Then, use spatstat’s nndist.pp3 function with p3 as an argument to get a list of distances.  You can then send the result back to shape; simply referring to a column that doesn’t exist yet will create it:

p3 <- pp3(x,y,z, as.box3(range(x), range(y), range(z)))
dists <- nndist.pp3(p3)
shape$nndist<-dists

The final bit of R code is very similar to the first:

writePointsShape(shape, file.choose())

This brings up the same dialog as before, which allows you to name the output file.  R will warn you that the file does not exist and ask you to make it, though you’ll need to make sure that you give the file a .shp extension.

If you look at the results of the summary function from earlier, you can see that the projection information did not get imported with the attribute data.  However, the coordinates are all the same as the original shapefile, so you can simply define the coordinate system to be the same as the original.  If you have the original loaded, its coordinate system will appear in the Layers list in the Spatial Reference selection dialog:

The attribute table for this shapefile includes new columns, including nndist, the one we just added:

You could use this new column to adjust the color or size of the points.  Although it is not possible, as far as I know, to create a 3D isosurface in ArcScene, you can approximate it by creating 3D buffers, which can use a field value as the radius.  Since nndist is small in clustered data, and we would prefer larger spheres for our clusters, it’s best to create a new field calculated from the first.  In the table options menu, add a field.  Make it a float value.  Right-click on the header after it is created, select Field Calculator, and enter your conversion.  After some trial and error, I found that the conversion that I liked best was 1/nndist*20000.

Once you create the column, you can use the Buffer 3D tool found in 3D Analyst Tools> 3D Features to create the buffers.  You’ll want to use the size field rather than a constant value for distance:

Now the clusters can be viewed in ArcScene:

So, that’s about it for this project.  These operations were accomplished with very little code; most of the work is done by the packages and functions already in R, and it’s possible to make your own as well.  The most difficult part of this process is learning how to make the various functions work.  A great resource is the help(“f”) function, where f is a function you are trying to use.  Running this command will bring up the documentation for that function, usually with examples.

Dr. Narayanaraj kindly provided some example R code that served as the starting point for my research in this project.  I’ll upload the files to the Data folder in the S drive.

The seismic data used in this exercise is from the Advanced National Seismic System catalog, maintained by the University of California Berkely.

 

 

If you read my posts from last week you would have noticed that I had a mini-breakthrough with regard to coming up with a method to seasonally adjust my input air pollution concentrations for my LUR model. This week I proceeded with using these seasonally adjusted annual concentrations in my LUR model. I predicted naphthalene concentrations at two sites in Eugene and compared the seasonally adjusted model with the no seasonal adjustment model. The results are below (and exciting!). The table demonstrates how, for each monitoring location, my predictive model’s accuracy is improved markedly with the seasonal adjustments (e.g. 14.6% vs. 22.6% and 1.1% vs. 10.3%). This encouraging and provides preliminary evidence that modeling temporal variation will improve annual estimates of air pollution exposures.

Monitoring Site	LRAPA Observed	OLS + Seasonally Adjusted Estimated Concentration (% Difference)	OLS Without Seasonal Adjusted Estimated concentration (% Difference)
Bethel Site	0.0789	0.0674 (-14.6%)	0.06112 (-22.6%)
Amazon Site	0.0533	0.0527 (-1.1%)	0.0478 (-10.3%)

I proceeded with performing the seasonal adjustments I proposed at the end of my last blog post (see “Temporal Data” post). Briefly, I used the ratios of mean monthly naphthalene concentrations over the annual naphthalene concentrations for all air monitoring data (i.e. both sites combined). These ratios were then used as adjustment factors to simulate monthly observations for the NATA data that I have. The monthly simulated data was then averaged over a year-long period to obtain the “seasonally adjusted” annual estimates.

The graph below displays the seasonally adjusted versus the unadjusted naphthalene annual concentrations. The adjustment factors caused the annual estimates to increase slightly. These seasonal estimates will now be used as my simulated input data for estimation of annual concentrations for my LUR. model.

 

I have been busy trying to figure out how exactly to incorporate seasonal variation of naphthalene ambient air concentrations into my Land Use Regression (LUR) model. To start off, I obtained time series data of naphthalene concentrations from the Lane Regional Air Protection Agency (LRAPA), just to see if there is indeed a significant amount of variation temporally. The graph below portrays the variation of mean monthly naphthalene concentrations for two air monitoring sites in Eugene (Amazon Park and Bethel) between April 2010 through April 2011.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

While this graph above is somewhat informative, it does not quite describe the relationship in full. The graph below presents the data disaggregated by month without averaging over the months, which indicates a great deal of variation within months and between months. This suggests a complex temporal variation that would likely require a non-linear type model. I therefore decided to use all of the data, rather than aggregated means, to better characterize the temporal relationship (i.e. still retaining the within month variations).

Next I modeled the change in monthly log transformed naphthalene concentrations (since naphthalene was clearly not normally distributed and the transformation lead to a log-normal distribution) using month as a predictor and air monitoring location as an additional predictor. The seasonal predictor was fit with a smoothing function using generalized additive model in order to allow for the obvious non-linear trend between month and naphthalene concentrations.

 

 

 

 

 

The following graph displays the relationship between month and log naphthalene concentration using the GAM function in R to allow for a non-linear trend. The use of a smooth function proved to be a better fit given the obvious visual trend and that seasonal adjustment was only significantly predictive once treated as a non-linear trend (the linear model was not significant for season [again not surprisingly!]). This work confirmed the notion of a complex relationship between pollutant concentration and month, even after adjusting for the air monitoring location.

As far as next steps is concerned, I was thinking of using this temporal GAM modelling approach for my LUR model. I would use the historical dataset from LRAPA to simulate temporal observations (i.e. monthly)  for my NATA data set. If I’m able to simulate this dataset, it would enable me to include time in my LUR model (using a GAM for month), and therefore enable me to see if adding temporal relationships will improve my predictive LUR model. I view this as a novel way of incorporating temporal variation as an explanatory variable. One approach could be taking the ratios of monthly to annual observed LRAPA naphthalene concentrations, and multiplying these monthly ratio with the NATA annual estimates. The simulated data would then be inputs for performing the LUR GAM model.

LiDAR point information is usually available as a set of ASCII or LAS data files:

ArcGIS only supports LAS data files; to use ASCII LiDAR data with Arc, you’ll need to use an external tool to convert to LAS.

LAS files cannot be added directly; they must be combined into an LAS dataset that sets a consistent symbology and spatial reference for the entire collection.  To create a LAS dataset, go to ArcCatalog and right-click the folder you want to store the dataset in, and select new>LAS Dataset:

Note that you will need 3D Analyst or Spatial Analyst activated to do this.  I recommend checking all the extensions to be sure your tools run the first time.

Right-click the dataset in ArcCatalog and choose Properties.

Ensure your dataset has an appropriate coordinate system in the XY Coordinate System tab, which you’ll need to get from the metadata for the LiDAR.  Next, activate the LAS Files tab and click Add Files or Add Folders.   Once you are done adding files, activate the Statistics tab and press the Calculate button.

At this point, you can import your data into either ArcMap or ArcScene.  There are pros and cons to both.  As far as I’ve been able to determine, it is impossible to plot 3D point clouds in ArcMap with a DEM or map base.  This is possible in ArcScene, and it is also possible to color points according to intensity in 3D view, but unlike in ArcMap there is no ability to adjust point size, and very limited ability to adjust colors of points, at least as far as I’ve been able to determine over the last few days.

Some LAS datasets will include RGB information in addition to intensity, which allows 3D true-color visualizations.

This image shows the Middle Lookout Creek site in the HJ Andrews Experimental Forest as a point cloud colored by intensity in ArcScene.  The creek and some log jams are visible on the right, and a road is visible on the left.

To convert LiDAR points to a DEM, you’ll need to convert the dataset to a multipoint feature class first.  Open ArcToolbox and go to 3D Analyst Tools> Conversion>From File>LAS to Multipoint.  Select the specific LAS files you want.  It’s likely that you’ll want to use the ground points only.  The usual class code for ground points is 2, but you’ll want to check this by coloring the points by class.

Once you’ve created the multipoint feature, you need to interpolate the values between the points.  There are several ways to do this, with advantages and disadvantages.  Popular methods are Spline and Kriging.  Spline generates a smooth surface but can create nonexistent features, and doesn’t handle large variations over shorter than average distances very well.  Kriging is hard to get right without experience in what the parameters do, and can take a long time to achieve the best results, but attempts to take spatial autocorrelation into account.  In general, Kriging is better for relatively flat areas, and Spline is better for sloped areas.  Inverse Distance Weighting is popular, but produces results similar to a very poorly configured (for terrain anyway) Kriging interpolation.  I find that a safe but time-consuming bet is Emprical Bayesian Kriging, which can be found in the toolbox under Geostatistical Analyst Tools> Interpolation, along with a few other advanced interpolation methods that I am not as experienced with.  If anyone else is familiar with these, I’d welcome a post explaining how best to use them.