NCSU GIS 714:
Geospatial Computing and Simulations

Introduction to human mobility pattern analysis

We will cover the following tasks:

  • per cell correlation maps between variables: binned or interpolated participant data with visualization over lidar DSM, computing viewsheds and distances to attractions
  • basic vector operations (clipping), computing and analyzing kernel density maps of participants travelling patterns, distances to attractions (Overijssel)
  • viewshed composition summaries and viewshed maps for cycklists, relations to participant data (CHIPS)
  • working with own data, identifying suitable data and tools, setting up a project location

Software

Data

We will use open data from Netherlands in two coordinate reference systems (CRS):

  • EPSG 4326: geographic coordinates longitude,latitude referenced to WGS84 elipsoid used by most GPS
  • EPSG:28992 (Amersfoort / RD New ): projected coordinate system used in NL (onlique stereographic projection with Bessel elipsoid)
Note that there are other official CRS used in NL, such as EPSG 7415 Amersfoort / RD New + NAP height, or EPSG 4289 which is lat/long registered to Bessel elipsoid 1841. Always use software with projection conversions support to process geospatial data.

Learn more (slides with brief description and related links):

Netherlands provides standardized geospatial data through PDOK web portal.

Before starting it is useful to go through GRASS GIS Quickstart and learn some GRASS GIS terminology for its data organization.

The workshop data should be downloaded using links provided through email

  • nl_epsg28992_amersfoort_workshp.zip and nl_epsg28992_amersf_overijssel.zip datasets
  • workshop_breda_external.zip

Now, let us organize our work:

  • Move the unzipped folder nl_epsg28992_amersfoort_worksh (this is our project location) into the grassdata directory (grass database) created automatically when installing GRASS GIS.
  • Then create a directory giswork outside grassdata. We will use it for external (non-GRASS) geospatial data, for example, data downloaded from on-line repositories or or your georeferenced participant data. You can also use it for various supporting files (color ramps, series lists, etc.) and your outputs. Move the workshop_breda_external.zip into your giswork directory.

Start GRASS GIS

Start GRASS - click on GRASS icon or type
grass

In data catalog select grassdata as your database and nl_epsg28992_amersfoort_workshp as your location (project) which has several mapsets (shown as subdirectories). Double-click on nuenen_outdoor and select Switch mapset - nuenen_outdoor will be shown in bold and we can start working with our data.

Note: You can work in GRASS GIS through the graphical user interface (GUI) or using command line interface (CLI) by typing (or copy-paste) commands in the Console or for linux/mac Grass Terminal You can write your workflpws as python scripts and share them in Jupyter notebooks (not covered yet here).

Display the available data

In data catalog (Data tab) display the data by double click on the listed map layer - it will be added to the display and the layer manager (Layers tab). To display the aerial photography and urban topography select (double click) ortho_vinc_patch05, dsm_05mfill, dsm_05mfillsh . To manage the displayed layers switch to Layers and you can change the opacity of dsm_05mfillsh (right click and select Change opacity level) to create a colored shaded relief. On the bottom of the Layer manager window you can see the command line (CLI) equivalents of displaying data: 
d.rast ortho_vinc_patch05
d.rast dsm_05mfill
d.rast dsm_05mfillsh
in a similar way we can display the points that represent our exhibits and the participant data, to use CLI you can paste these commands into Console or just double click the layer in the Data catalog: 
d.vect map=exhibits_all_om@nuenen_outdoor icon=basic/circle size=10
d.vect nu_all_hres size=2

To transform our point data into a continuous point density representation we can compute a kernel density map - it approximates the pattern of most used areas in terms of number of people and time. First we define the spatial extent and resolution for our analysis, then we run the kernel density tool and finally we adjust the output color ramp and display only values above min threshold: 

g.region raster=dsm_05mfill res=0.5
v.kernel nu_all_hres out=kd_nu_all_hres rad=7 kern=epanechnikov
r.null setnull=0 map=kd_nu_all_hres
r.colors kd_nu_all_hres co=gyr

We can import and rasterize our participant data directly from the csv files. To access the files without the need to specify full path we change working directory to giswork
Settings > GRASS working environment > Change working directory > select giswork
or type cd (stands for change directory) into the GUI Console and hit Enter: 

cd

Now you can use the commands requiring external supporting text files without the need to specify the full path to the file.

To perform processing and analysis of our participant data we set the spatial extent to existing raster dsm_05mfill and the grid (raster) resolution to 2m. In our example below we rasterize (bin) the point data upon import using mean phasic z computed from points in each grid cell for each participant: 

g.region raster=dsm_05mfill res=2
r.in.xyz in=opt11h.csv out=opt11_mean_2m meth=mean sep=, x=3 y=4 z=13 skip=0
r.in.xyz in=opt13h.csv out=opt13_mean_2m meth=mean sep=, x=3 y=4 z=13 skip=0
r.in.xyz in=opt14h.csv out=opt14_mean_2m meth=mean sep=, x=3 y=4 z=13 skip=0
r.in.xyz in=opt20h.csv out=opt20_mean_2m meth=mean sep=, x=3 y=4 z=13 skip=0
r.in.xyz in=opt21h.csv out=opt21_mean_2m meth=mean sep=, x=3 y=4 z=13 skip=0

number of participants per each grid cell 

r.series in=opt11_mean_2m,opt20_mean_2m,opt14_mean_2m,opt13_mean_2m,opt21_mean_2m out=opt_count_2m meth=count
Note that if our data are in csv format we can import them as points using v.in.ascii or v.in.csv commands and then rasterize them (perform the binning) using the add-on tool r.vect.stats. Use g.extension tool to install an addon. To compute the correlation between the survey data (recommendation level) and the physiological response (mean phasic z) we first order the participants by recommendation level 
11 - 0
20 - 7
14 - 8
13 - 9
21 - 9
and then we compute the mean and regression for a subset for grid cells with at least 3 participants 
r.to.vect opt_slope out=opt_slope
r.series in=opt11_mean,opt14_mean,opt20_mean,opt13_mean,opt21_mean out=opt_mean meth=mean
r.series in=opt11_mean,opt14_mean,opt20_mean,opt13_mean,opt21_mean out=opt_slope meth=slope
r.series in=opt11_mean0,opt20_mean7,opt14_mean8,opt13_mean9,opt21_mean9 out=opt_slope meth=slope
r.series in=opt11_mean0,opt20_mean7,opt14_mean8,opt13_mean9,opt21_mean9 out=opt_detcoef meth=detcoef
To compute weighted mean phasic we first change nulls to zeroes and compute the weighted mean using map algebra tool r.mapcalc: 
r.null opt11_mean null=0
r.null opt13_mean null=0
r.null opt14_mean null=0
r.null opt20_mean null=0
r.null opt21_mean null=0

r.mapcalc "opt_partwmeanph = (opt11_n*opt11_mean + opt13_n*opt13_mean + opt14_n*opt14_mean + opt20_n*opt20_mean + opt21_n*opt21_mean) \
/(opt11_n + opt13_n + opt14_n + opt20_n + opt13_n)"
To create continuous representation map for mean phasic z and correlation coefficient we first set the grid resolution to 0.5m, then convert our 2m resolution r.series outputs to points and interpolate the 0.5m resolution maps using smoothing spline with tension tool v.surf.rst: 
g.region res=0.5
r.to.vect -z in=opt_partwmeanph out=opt_partwmeanph type=point
v.surf.rst in=opt_partwmeanph elev=opt_partwmeanph_rstdef
r.colors opt_partwmeanph_rstdef rules=co_phasic_pmean.txt

r.to.vect -z in=opt_detcoef out=opt_detcoef type=point
v.surf.rst in=opt_detcoef elev=opt_detcoef
r.colors opt_detcoef rules=co_detcoef.txt
#overijssel 
switch to mapset in a different location
Workflows for this part and more Nuennen and CHIPS analysis are available as text files and selected subsets will be added here. To add the newly output pngs to a web map continue by following the instructions found in the next tutorial here https://github.com/ncsu-geoforall-lab/grass-mapbox-tutorial