NCSU GIS/MEA582:
Geospatial Modeling and Analysis

Comments, questions and answers for assignments

Data acquisition and integration

Should we list and provide information about all the files in each Location and geodatabase, or a subset?

As the instructions for HW indicate you should be able to write the text of your report on 1 page (500 words), so rather than listing properties for each file you should provide a general description that will characterize the data. For example, after looking at most of the data in the data set, are you able to generalize some common properties? are all the data local? regional? do they have all same level of detail? do they represent natural or anthropogenic features? how are they georeferenced? The goal in this assignment is to ensure that the students understand the basic properties of digital geospatial data and can make general observations and conclusions based on broad range of specific technical information. The reports will get more meaningful for more complex assignments - the first ones are just a warm up. In the approach section - do not write about software issues, focus on general methodology, in the results sections focus on what is interesting in the data and the derived maps. Put software issues, including fixes and suggestions for improvements into the discussion section.

When I open GRASS, nc_spm_08 doesn't show up under "Project Locations", although it is downloaded and unzipped on my computer.

Make sure that when you unzip your nc_spm_08 is not hidden under another nc_spm_08. Some zip software tools create an extra directory by default, so you need to specify that you do not want a new directory created. Otherwise you end up with something like C:\grassdata\nc_spm_08\nc_spm_08 and if you then give GRASS path to GRASS GIS database as C:\grassdata it still does not see the correct nc_spm_08 (it checks for the folder PERMANENT under nc_spm_08 to find out whether it is a valid GRASS Location) and then it would not show up. To avoid future confusion, remove all nc_spm_08 folders and unzip nc_spm_08 again, making sure that no extra directory is created.

Is there a way to reorder the layers so that this information does not become hidden?

yes - just grab a layer with the mouse and move it where you want it to be.

data were in ... geographic coordinate system

there are no data in geographic coordinate system in the provided data set - everything is in State Plane coordinate system (datum:NAD83, projection:LCC)

I expected the data concerning projections, coordinate systems, etc. to be presented in GRASS GIS more like ArcGIS.

In your GRASS Location, everything is in the same coordinate system. You find the information about current Location by running g.proj -p.

I did not get the XY feature to work

Please be more specific. What did not work? Did you get an error message? More detailed information helps to to improve the instructions or, if there is a bug we can file a bug report. (And you will learn how to prepare bug reports too!)

I was not able to figure out how to create a copy of data in GRASS GIS.

In the g.copy dialog (which you get, e.g. by typing g.copy in Command Console and pressing Enter) Find the option (input field) for raster or vector, e.g. raster files to be copied. Then you need to select the "from" file using the browser and type name of "to" file.

Why there is no space after comma in g.copy raster=slope,myslope?

The comma syntax (as in to,from) is used in GRASS GIS a lot, so you should remember not to put space after comma in this case. Basically, the reason is that the same syntax is usable in command line where the space would have a meaning. Generally speaking, GRASS GIS is aiming at being easily used when scripting or programming in different languages, so this requires to be strict sometimes. Similarly, ArcGIS raster map calculator, it requires spaces here and there, e.g. around operators. On the other hand, GRASS GIS raster map calculator does not require spaces because there is just no advantage in requiring them. However they are considered as a good practice when used around operators and after comma.

I am getting Error "mygeology" map layer already exists.

Give your output a new name or you can force overwrite by adding --o to your command (this applies to practically all GRASS commands). There is an equivalent overwrite checkbox in each module's GUI.

Image rendering in GRASS GIS is slow.

check your region and make sure you have it set up as given in the instructions and then zoom to computational region - these are small data sets so they should be displayed instantly when added to the map layer. Also, do not keep large number of layers in the Layer Manager, start each assignment with empty list and if you accumulate many layers (10+) remove those that you don't need.

Do I have to use 3D views?

Yes, please save the image of your 3D view (aka NVIZ) so that I can see that you can work with it. We will be using 3D views a lot. Save the NVIZ image as TIFF file, make sure you have a TIFF support installed on your computer, otherwise when you try to open the file, it will be black. If you can't work with TIFFs, just take a screenshot.

Should I include my name inside the report or to the filename?

Include your name on the first page of your report and see Course logistics for filename conventions.

Can I use data with different projections at the same time?

GRASS GIS manages integrity of the data by allowing to work with data in one projection at the time (i.e. it considers on-the-fly projection as such a bad practice, that it doesn't allow it at all) and it also requires you to import the data, to have all the data in one format, e.g. to store vector data in GRASS topological format.

How the geospatial data are organized in GRASS GIS?

The data are stored in GRASS Location and you can work with one Location at the time, it is similar concept to database (speaking, e.g. about relational databases). GRASS Location is divided into Mapsets. This is a way how to organize your data. Common data are in special PERMANENT Mapset. Other Mapsets contain the data related to one task, or project, for example one class assignment. You always work in one Mapset and you can access data in PERMANENT Mapset, however you cannot change them which protects you from deleting or overwriting the data accidentally. Optionally, you can access the data in other Mapsets within the same GRASS Location but again, only for reading. Obviously, you cannot access Mapsets from other GRASS Locations because they might not be in the same projection. This also means that it is always defined where the new data will be created and that's the current Mapset. The result is that once the data are imported, one can focus on analysis, not data formats. In fact, you can link the data from different sources such as GeoTIFF or PostGIS database but that's a different story.

Data display and visualization

I'm not sure what setting the grid parameters "size=0:02" did.

It draws the lat/long grid - you should read the numbers along the image broder to find out. Maps used to provide both x/y grid in meters and lat/long in degrees:min:sec so that both cartesian and geographic coordinates could be derived.

What is the ArcGIS equivalent to setting the computational region in GRASS GIS?

rather than clipping the data set, which would be a pain if you need to use many data sets with different spatial extent and if you need to work on different regions (you would be also creating a lot of extra files) is under: Environments > General Settings > Extent

When I run GRASS command for the second time, it gives me an error?

You probably already created the raster or vector maps by the first run and now GRASS don't want you to loose the already cleated data. Deleting or replacing the data usually requires additional action to be taken by the user (i.e., you) and it is the preferred behavior over simple removal of the previously created data. Check the behavior of your system's file manager when you try to copy a file somewhere and the file already exists and same would apply, e.g. to saving a document in a text editor.

Why we just don't use very high resolution for all raster data?

When you have low resolution (big cells) it does not necessary mean that it is worse. Sometimes we can even rescale hight resolution data because we don't need the precision because other data for analysis are in the lower resolution or because we are hitting some storage limits. Re-scaling to higher resolution (smaller cells) does not give use more precise information, although the image might be smoother. Similarly, when comparing raster and vector, it is important to realize that there both has limited (spatial, yx) precision and although the area boundary in vector looks nicer it might be much less accurate then the raster representation depending on how the data were acquired.

GRASS GIS Map Display is not using my system default font. Is there any way to change it?

You can change the default font used in Map Display by going to the main menu > Settings > Preferences > Map Display tab and click Set font button. Select name of any font you like.

The maps composed in GRASS GIS Map Display are not good enough for me. Is there something else I can use?

Advanced cartography in GRASS GIS is done in completely separate application/module called GRASS GIS Cartographic Composer (which is a GUI frontend for ps.map module). However, some users stay just with the GRASS GIS Map Display because it is easily scriptable thanks to display modules (d.*) or they combine it with Inkscape. Some other GRASS GIS users, however, use QGIS for high-quality cartography tasks.

I cannot change or create data in my other Mapsets, only in the current one.

That's how it should be. GRASS GIS is always connected to one GRASS Location inside a GRASS Database. GRASS Location contains multiple GRASS Mapsets and you can access data in any of these Mapsets. However, you can change (create, modify, delete) the data only in a current Mapset (that's why you start GRASS session with combination of database-location-mapset). Data can be copied into current Mapset using g.copy if you need to modify them (e.g. attributes of a vector map). By default, you see data only in the current Mapset and in the PERMANENT Mapset (that's the special thing about PERMANENT Mapset).

When using GRASS GIS, can I access data in different Mapset or Location?

If desired, you can show data from other Mapsets (in the current Location) using Settings > GRASS working environment > Mapset access. This helps to maintain safety of data and gives a system to organize the data. Data in other Locations are accessible only if you start GRASS session in the Location. Accessing data from multiple Location wouldn't be safe because data in different Locations are (can be) in different coordinate systems (while all Mapsets in one Location contains data in one particular coordinate system).

Global, zonal and focal operations, map algebra

What is happening when (computational) region is set?

It sets the spatial extent AND resolution at which all raster operations are done. If any of the rasters has different resolution, it is quietly resampled using nearest neighbor. It also sets the viewing area of the Map Display when the option Zoom into computational region is selected. It does not affect vector operations, so vector data need to be clipped if the operation is to be performed on a subset of the data. However, some vector modules can work with the region when it is advantageous or when rasters are involved.

what does the difference between SRTM and NED DEMs mean?

SRTM represents surface with vegetation and buildings, while NED is derived from bare earth lidar data, so it is lower in most areas. The issue is complicated by the fact, that in SRTM v1 surface water was assigned value 0 instead of NULL or actual value (there were no data), so lakes have the elevation value 0 leading the large negative difference. Both DEMs are properly georeferenced to the same vertical datum.

What is happening in the map algebra example with relative coordinates?

South-east corner is replaced.

How we got from the computer raster to the plane? How would I use it?

We are computing an equation for a plane here: ax+by+cz+d=0 which we can rewrite as
z = -1/c (ax + by) - d/c and use it to model, for example, a geologic fault.

NDVI computed with integers gives 0 in GRASS and -1, 0 in ArcGIS, what is going on?

The correct values for NDVI range between (-0.957, 0.979) so when the decimals are simply cut-off the result is 0, when the numbers are rounded down to nearest integer the result will be -1,0. Neither of them is correct.

Subsetting a raster using map algebra - what is happening there?

we set the spatial extent to orthophoto and resolution to 10m (that means that the 1m resolution of ortho is ignored). Then we compute a new raster by assigning it the same values as in the original raster only within the defined spatial extent. So the resulting map has the spatial extent of the orthophoto but the values and resolution are the same as in the original raster.

Buffers, cost surfaces, least cost path

Can I manually set region by drawing a rectangle on the map in display window?

yes, zoom to the area of interest using the Zoom in button and drawing a rectangle and then click Zoom options > Save display geometry to named region file, e.g. myregion. You can then set your region to this using g.region myregion or by running Config > Region > Set Region > Set from named region > myregion. Update for GRASS GIS 7: Find toolbar button 'Various zoom options' on Map Display toolbar, select 'Set computational region extent interactively' and then draw the rectangle.

What the internal structure is of a GRASS raster, where is category list stored?

We have talked about it in the first lecture, but you can find it in the GRASS GIS manual too. If you are really interested you can look to the directories yourself. Categories are under cats directory, color tables under colr, and history and other metadata in hist. These are all ASCII files so they are readable and editable as plain text. Never, ever edit these files unless you know what you are doing, you can corrupt your data, use the relevant commands instead.

I would have thought that I could have used the streets layer as a mask to contain all analysis to the roads system itself and still arrive at the same results. To my surprise, this was not the case. Why would areas outside the road system have such an impact on travel time and shortest path within this analysis?

I would have to see your results to say for sure, it should be the same if the firestations were all located on the roads, but they are not. They are located off-the road and they are not connected to the street network. So you need to get them from the points reprsenting the stations onto the road, either snaping them to the road (which may require manual aditing) or have the algorithm find the least cost path to the road by assigning off-road areas some low speed value. Depending on where the firestation is located in relation to the roads, the two approaches may lead to different paths.

Coordinates of firestations are not stored in the attribute table. Where do they come from?

v.out.ascii pulls them from the coor file that defines the geometry of each vector map and combines them with attributes from the SQLite table (or any other DBMS that you may use - GRASS supports several).

What are the friction values?

these are empirically selected values that measure how difficult it is to walk through cells with different types of cover, they increase the cost to pass a grid cell.

What is the best way to get the min and max raster value in GRASS?

Use r.info -r.

In GRASS I don't get any firestations for the given condition CVLAG < 0.1, I had to use CVLAG < 0.5 and the stations that I got were far away (S. East St. and W main St stations)

you may be searching in a wrong region, reset it with g.region swwake_30m -p

Terrain modeling

Is 2m or 6m resolution more suitable for DEM from the given lidar data?

if we are going to generate the DEM using binning, e.g. mean elevation for each cell, 6m provides more complete coverage, but as the range map shows we lose some detail. 2m resolution is not suitable for generating DEM by binning - the data are not dense enough. But if we interpolate (computationally more intensive task than binning) we will get a more detailed DEM.

mr lidar points in ArcGIS shape geodatabase were missing Z value

new version of the geodatabase with Z-value added is available on the Assignments web site

The image I see in 3D in GRASS GIS is just a lot of spikes or outliers, no digital elevation model.

When visualizing rasters in GRASS GIS 3D view, you have to switch off aspect and any other raster map layers which you don't want to show as surface. You can however, show any raster by draping it over a surface in 3D (by selecting method of coloring for a surface).

When I show the TIN in GRASS GIS 3D view, I get just large spheres or flat thick plate. How to view the 3D TIN.

When visualizing vectors in GRASS GIS 3D view, go through a Vector section of the Data tab and check values of all the fields. In case of TIN, you probably don't want to show any points or set the size of points to a small number. Note also that depending on vector type, you need to specify if you want to show the features on an existing raster surface (flat plane by default) or as actual 3D features. Switching off other layers beforehand in 2D is also a good idea.

When I show the raster map in GRASS GIS 3D view as a surface, I cannot see much detail and the view seems not precise, coarse or blurry.

Set fine mode resolution in 3D view tab, in the Surface part to 1 to get the surface in the actual (full, native) resolution of the raster.

What does the fine mode resolution in GRASS GIS 3D mean?

Fine mode resolution in 3D is resolution of resampling when using the raster in 3D, for example fine mode resolution 2 means that raster with resolution 3 is shown as if it would have resolution 6. This is for performance purposes. To see raster in its full (native) resolution, use fine mode resolution 1.

Spatial interpolation and approximation

How to design an experiment to test IDW with different parameters?

Run the command with the suggested parameters and then compare the results visually and statistically using shaded relief (hillshade), contours, univariate statistics and histograms. For example, you should be able to conclude that the lower the number of points, the rougher the surface The univariate statistics shows almost no difference illustrating the fact that summary statistics does not reveal subtle but significant differences in surface geometry.

Spatial interpolation and approximation: splines

GRASS defines spline method as regularized spline with tension, while ArcGIS defines spline method with two functions, regularized and tension. It seems the functions from these two software are different when both state it is a spline interpolation method.

There are many different spline interpolation methods based on what type of smoothness seminorm (also called roughness penalty) they minimize, see Table 1 in this paper (this paper has lot of typos as we were not given the proofs after the chapter was typeset, but it has a good overview of interpolation methods, note that in Table 1 regular should be regularized, and the spline by Hutchinson in Topo to grid is thin plate spline with tension, not just thin plate spline). Regularized spline with tension is equivalent to completely regularized spline but the others are different.

Another confusion I have is radial basis functions. Radial basis function methods are part of spline method? Or some are spline (thin-plate spline, spline with tension, and completely regularized spline), and some are not?

Splines are a special type of radial basis function methods, multiquadrics is another group RBF methods. But there are also splines that are smoothly connected piecewise polynomials that we have not covered at all (they are applied to meshes which are not very well supported in GIS).

And does spline with tension under radial basis functions create the same results from spline interpolation with tension type?

Spline with tension under radial basis functions and the one in Spatial Analyst are probably based on the same equations but the implementation may be different so unless it calls the same code, the results will be different ( e.g. the data are usually normalized to avoid numerical problems and each program does it differently, and the search for points may be different too - we covered that in the lecture.)

Spatial and Temporal Terrain analysis

In the time series analysis task my surfaces are fuzzy and the cutting plane does not work.

Make sure you set display to computational region and use fine resolution=1 in nviz. When using cutting planes both surfaces should have their resolution set to same value (e.g. fine=1), use constant plane for reference (-1m for the coast) and set its resolution to fine=1 too. When comparing two surfaces assign each of them different constant color and use top or bottom surface option to color the cutting plane. You can lower your light source and reduce its brightness to create contrast between the cross-section and the surface if you are coloring the cross-section using the top surface option.