When georeferencing using the GDAL georeferencer the mean error (root-mean-square error) is the standard value to evaluate how much error will be produced using these points and algorithm. However, how big this error can/should be is highly dependend on the quality of the input data. E. g. when I have a high-resolution aerial image with 1 cm ground resolution, I would not accept an RMSE of 1 m, while for a scanned old map with a ground resolution of 25 m this would be a very good result. But – as far as I can tell – how much error the user is willing to accept as tolerable is generally up to the user.
However, this is not a very comprehensible approach.
Is there any rule of thumb or actual formula or table where the RMSE-values acceptable for certain resolutions etc. are listed?
If not, I would be interested in how you approach this problem/what is a acceptable error for you?
I am currently trying to create an erosion model of the mound of Oymaagac Höyük (Samsun, Turkey) to investigate the relation between the distribution of the survey findings there collected and the potential erosion of the tell. After some bibliographic research, I decided to try two different methods:
- Combining the RUSLE formula with GIS tools as in Howland 2018 (Howland, M. D. et al. 2018, Quantifying the effects of erosion on archaeological sites with low-altitude aerial photography, structure from motion, and GIS: A case study from southern Jordan. JOURNAL OF ARCHAEOLOGICAL SCIENCE, 90 (Antiquity 84 2010) https://escholarship.org/uc/item/27q3f2vt);
- Calculating the Topographic Index with GRASS r.topidx.
As it concerns the 1st method, the Revised Universal Soil Loss Equation (RUSLE) is a mathematical model used to predict rates of soil erosion caused by rainfall and land use. I am still working on the calculation of the RUSLE factors (A= R K LS C P, where A = average annual soil loss, R= rainfall erosivity factor K = soil erodibility factor, L = slope length factor, S = slope steepness factor, C = cover-management factor, and P = supporting practices factor). If R can be easily computed via GRASS r.usler, and L and S gained through the DEM, K is more complicated than expected, as I do not have the data necessary for the calculation of the K factor (that would be otherwise easier to calculate - GRASS r.uslek). I only have the rough archaeological description "clay loam".
Does anyone have an idea about how to correctly calculate it?
Considering the difficulties in calculating the RUSLE factors, I decided to try with the Topographic Index (GRASS r.topidx). Assuming that topography drives the flow, the topographic index predicts the flow accumulation in topological low spots. As the slope increases, the potential of accumulation decreases. Here you can see the final result.
Does anyone know why I get negative values? I have read that I am supposed to gain only positive values (I am using high-resolution DTM - pixel size 15 cm-). And moreover, can I actually connect it with the erosion of the mound and use it to predict where the survey findings will most probably accumulate?
Coordinate reference systems and projections are at the heart of any GIS. Unless something went very wrong you usually don't see much of it, once you have settled for one system. However, I stumbled upon a CRS-related issue in QGIS, which is quite hard to spot on first glance and for which I currently don't find a solution.
1) If you use the Field Calculator to calculate certain properties of your geometries, you may have noticed, that you have two functions for most calculations, e.g. $area and area($geometry). On first glance the results appear to be the same, but if you compare the results of two calculations directly you will notice a difference, especially for larger polygons. For small features differ only by a few square centimeters. The same is the case for $length and length($geometry) etc.
The issue becomes clear, if you look into the description of the two functions:
$area: "Returns the area of the current feature. The area calculated by this function respects both the current project's ellipsoid setting and area unit settings. For example, if an ellipsoid has been set for the project then the calculated area will be ellipsoidal, and if no ellipsoid is set then the calculated area will be planimetric."
area: "Returns the area of a geometry polygon object. Calculations are always planimetric in the Spatial Reference System (SRS) of this geometry, and the units of the returned area will match the units for the SRS. This differs from the calculations performed by the $area function, which will perform ellipsoidal calculations based on the project's ellipsoid and area unit settings."
This isn't a problem per se, since it only means that you have to choose the right function for your calculation. For most problems $area should be the right option. Note that the results differ according to the CRS you have set (compare the results in the image above (EPSG 3857) to the ones below (EPSG 25833)).
2) However, this shows clearly how much the ellipsoidal and the planimetric calculations differ. This is an issue, if you want to construct features using the Advanced Digitizing Toolbar > CAD-Tool. As you may know, this tool does not work for geographic coordinates, so you will have to use a projected coordinate system. This would not be a problem if the CAD-Tool would be able to use ellipsoidal distances, but this is not the case.
The polygon in this image has been constructed using the CAD-Tool and with a fixed distance of exactly 100 m and a fixed angle of 90° in EPSG 25833 (note that the results differ a lot if you use e.g. EPSG 3857). The results of the calculations show that the construction was done planimetric and therefore differs from the real result on the curved surface of the earth.
Now my questions:
1) Is there a way to construct features in QGIS that take the ellipsoid into acount?
2) If not, which of the results can be taken to be "correct"?
3) Is the difference of in this case 0.05% a real problem for the daily archaeological practice?
In our last meeting Benjamin gave us a short overview of the open source software Survey2GIS provided by the Cultural Heritage Authority of the Federal State of Baden-Württemberg. Survey2GIS allows to import survey data from a total station directly into GIS without having to do a lot of manual work. It resembles TachyCAD in its functionality, but imports into GIS instead of CAD and it only offers post processing. Using codes on the total station allows to import e.g. a polygon data set directly into GIS where it will be represented as a correct polygon as well as offering a variety of attributes.
The software is available on GitHub currently with the stable version 1.5.0 and the 1.5.1Beta which mainly adds more functions of automatic topological corrections.
There is no installation needed, the software runs on Windows, Linux and MacOS, although the latter requires some additional installations. The download folder also contains a detailed user manual in German and English as well as test data.
The key feature of Survey2GIS is the parser file, that can be adapted very precisely to your survey project, including the measurement codes, skipping certain data etc. This allows to tailor the software exactly to your project instead of the other way around.
Although it is aimed at a usual 2D GIS since version 1.5.0 an additional feature allows to import profile data as well. To do this, Survey2GIS switches the Y- and Z-axis. This function is by now offered by other software as well and can also be done by hand in QGIS, but Survey2GIS allows you to do it automatically and does not require recalculations of coordinates by hand or additional fix points. To import profile data, go to "Extra" and switch the default Koordinatesystem "World" to "XZ-Orientation".
If further questions arise, don't hesitate to use the mailing list for them.
Is it possible to compare two versions of a shp-layer, mostly automatically?
I want to know whether two vector point layers are identical or at least: is layer A part of layer B, concerning not only the point coordinates but also their attributes. Have attributes (like objecttype, date, size, owner etc.) changed? The data consist of up to some hundred items with 10-20 attributes so that comparing by hand is not the first thing to do.
An extended solution would be a function similar to MS Word "compare two versions of a document".
Any ideas around there?
What is Blender?
Blender is a free and open source 3D Software. It excells at its field and managed to get funds and is a very good example for successfull FOSS software. It can be used for 3D Modelling, Rendering, Animation and Architecture. You can imagine that this software may also be interessting for Archaeologists for different tasks.
What can i do with Blender?
I want to show you two things you can do with blender as an archaeologist. First will be a workflow how to create a hillshade in 3D and second how to use the addon BlenderGIS to explore a region. Blender is also a great tool when you work with structure from motion, so its a good idea to take a look at this software.
Hilldshades
The visualisation of 3D surface forms on a 2D was a challenge since the beginning of mapmaking. In recent times the art of hillshading has been lost, since it was easier to rely on algorithms like r.relief in grass or the GDAL hillshade, which only take a few parameters input.
Blender allows the free movement of lights, the placement of lights, simulation of sunrays and the usage of additional lights. That can really spice up your hillshading.
The tutorial for that workflow is mainly this documentation: https://somethingaboutmaps.wordpress.com/2017/11/16/creating-shaded-relief-in-blender/
It also shows advacned techniques and what blender can do.
BlenderGIS
You can also use basic GIS tasks in Blender with the addon. Its really simple to use. Its also well documented:
https://github.com/domlysz/BlenderGIS/wiki/Quick-start
You can get a decent landscape render or even modern urban regions in only a few minutes.
This video might give you a better expression of what is possible: https://www.youtube.com/watch?v=Mj7Z1P2hUWk
The downside of BlenderGIS is that it only gives you some "eyecandy". But a nice looking render of a region might at least impress people that have potential funds for your project.
Maybe you will need help with your GIS Project and cant wait very long, so there are two options i can suggest for fast answers or input that might help you:
Discord - QGIS inoffical Server
Discord is kind of the next generation of irc chat, so you just click the link and can start using the chatroom, you dont even need a account: https://discord.gg/AjNcVcF
Reddit - r/qgis, r/gis, r/geospatial, r/askgis
Reddit can be very helpfull if you look at the right places. r/qgis is a really nice, but smal subreddit which can help you with most problems with qgis.
r/gis is more focused on other GIS Softwares like ArcGIS.
Stackexchange
Of course if you have enough time you will get the best help at https://gis.stackexchange.com/
As you all may know, Grass GIS is a great software for academic purposes and using it in an environment like R/RStudio is a very good idea, because you can share your code and other people can reproduce the tasks you did.
Unfortunatly the r-package used for this interactiuon between R and Grass can have some problems when using Windows as an Operating System.
Luckily the creator, Roger S. Bivand provided a safe solution, that should be used initialy so you can work without a problem. Other users at the support board also added some lines and i hope it this code can help maybe save you some time in the near future:
First of all, you need to find your GrassGIS installation folder. When using Windows i highly recomend using the OSGEO4W installation, because it provides you some other nice tools that may be needed. So your installation folder Grass is perhaps included in OSGEO4W64, but you better check if this is true. Otherwise you need to change the path in the following code. Also you need to check which version of grass is installed and which version of python is availabe in OSGEO4W, otherwise it may not work as intended. For me its Python37 and grass78. Be aware that grass may have different folder path for different versions like "grass-6.4.3", "grass78" or even "GRASS GIS 7.8" (having dots in filenames/paths is just not a good idea).
library(rgrass7)
use_sp()
osgeo4w.root<-"C:\\OSGEO4W64"
Sys.setenv(OSGEO4W_ROOT=osgeo4w.root)
# define GISBASE
grass.gis.base<-paste0(osgeo4w.root,"\\apps\\grass\\grass78")
Sys.setenv(GISBASE=grass.gis.base)
Sys.setenv(GRASS_PYTHON=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\bin\\python.exe"))
Sys.setenv(PYTHONHOME=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\Python37"))
Sys.setenv(PYTHONPATH=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\apps\\grass\\grass78\\etc\\python"))
Sys.setenv(GRASS_PROJSHARE=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
Sys.setenv(PROJ_LIB=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\proj"))
Sys.setenv(GDAL_DATA=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\gdal"))
Sys.setenv(GEOTIFF_CSV=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\share\\epsg_csv"))
Sys.setenv(FONTCONFIG_FILE=paste0(Sys.getenv("OSGEO4W_ROOT"),"\\etc\\fonts.conf"))
# call all OSGEO4W settings
system("C:/OSGeo4W64/bin/o-help.bat")
# create PATH variable
Sys.setenv(PATH=paste0(grass.gis.base,";",
"C:\\OSGEO4~1\\apps\\Python37\\lib\\site-packages\\numpy\\core",";",
"C:\\OSGeo4W64\\apps\\grass\\grass78\\bin",";",
"C:\\OSGeo4W64\\apps\\grass\\grass78\\lib",";",
"C:\\OSGeo4W64\\apps\\grass\\grass78\\etc",";",
"C:\\OSGeo4W64\\apps\\grass\\grass78\\etc\\python",";",
"C:\\OSGeo4W64\\apps\\Python37\\Scripts",";",
"C:\\OSGeo4W64\\bin",";",
"c:\\OSGeo4W64\\apps",";",
"C:\\OSGEO4~1\\apps\\saga",";",
paste0(Sys.getenv("WINDIR"),"/WBem"),";",
Sys.getenv("PATH")))
# initial again to be sure
use_sp()
#create a location so you can start with grass tasks
loc <- rgrass7::initGRASS(gisBase=grass.gis.base,
home=tempdir(),
mapset='PERMANENT',
override=TRUE)
I have also used Linux (Manjaro KDE) and it was not necessary to have this code for rgrass7 to work. rgrass7 just found the grass files on its own.
Run this code anytime you start up your R Environment.
Today we had our second GIS-meeting online. Attendance was still limited, but at least we have found a rhythm now. Today's meeting was short but informative.
Extracting islands from DEM
We discussed the May topic at length, which could not be covered last time. The solution is posted under the original question.
Raster reprojection in QGIS
Fabian asked an intriguing question. As most QGIS users know, there are at least two ways to reproject a raster layer in QGIS.
1) 'Raster > Projection > Warp' using this tool you can reproject a raster file to a new CRS and save it to a new file.
2) 'Export > Save as' when saving a layer you can assign a new CRS.
You can also use other external tools from GRASS or SAGA to reproject. However, when comparing the results of these ways, you will see that the reprojection results are slightly different. This effect is known in the community (https://gis.stackexchange.com/questions/214757/raster-reprojection-in-qgis-leads-to-different-results), but there appears to be no detailed explanation so far. Apparently the different ways use different tools for the reprojection. We do not know, which one gives the best results. Further research (e.g. QGIS-User-List) would be useful.
Ideas for next meetings
Fabian explained his troubles to connect GRASS 7 with R. He will prepare a little presentation on that topic for the next meeting. We decided to dedicate the next meeting to interfaces of GIS to other software packages.
Martin asked for an introduction to spatial network analyses. This will be our topic for the august meeting. We will prepare a data set and Benjamin agreed on a little intro to that topic. Benjamin also showed us the GAMA-Project, which is an open source software for agent based spatial modelling (https://gama-platform.github.io/wiki/Home). This tool can also be used for network analyses.
High quality test and open data sets are available for free from the Mexican governmental surveying office (https://www.inegi.org.mx/datos/). The problem is, that so far non of us has any archaeological data to combine with these.
We met on May 29th on the platform Big Blue Button. Attendance was a bit sparse, but we had interesting topics.
Link points to fotos
Wassim formulated his problem, that he wants to create a "fotopoints" layer for the ancient city of Bosra (Syria). For the DAIs digital map of Bosra the images of the city from the image archive Arachne are supposed to be linked to points on the map, representing the point of view of the camera. A similar linkage was earlier established for the digital map of Palmyra (https://arachne.dainst.org/project/palmyra-gis).
The solution we took for Palmyra was to add the URL of the image in the Arachne to the attribute table. Since every image in the arachne has its own ID and the URL always has the same structure, this is a quite straightforward matter. It is even possible to store just the ID of the foto and create the URL from that on the fly. For this example I created just a temporary layer representing some imaginary building and a second layer representing three viewpoints to this building (here represented by an SVG symbol for viewpoint from the QGIS standard svg-library).
To put in the Arachne-ID and/or the URL you will have to create a suitable field in the attribute table (the layer has to be in editing mode). For the ID a integer field will do, for the URL you will need a text field. The following image shows some random IDs and their corresponding URLs. The fields names are up to you of course.
To open this link directly from QGIS you have to create an "action" for the layer. To do so open the layer properties and select the tier "Actions". Here you will have create a new action by clicking the "+"-Button. Now you can enter the command you want to run with this action. In our case this would look like this.
Now you will have to select the "Run feature action"-Tool in the toolbar. You may have to select the dropdown menu first to select your "Open_URL"-Action. If you have the tool activated and click on a point with a valid URL, it should open this URL now.
The second question was, how to make the point rotate to a certain angles, as to represent the angle of view of the original foto on the map. Here we have to create a new attribute again. Let's call it "rotation". It is useful to make it a decimal number, with two decimal places, if you want to be extra precise. But a 3 digit integer should also do it. Your rotation is stored in degree with 0/360 pointing north.
Now all you have to do is to use this field to rotate your points. This is done in the layer properties in the tier "Symbology". Go to "Rotation" and activate "Data defined override". Here you have to select "Field type" and then select the field you want to use for the rotation value, in our case "rotation". Now the points will be rotated according to the values in your attribute table.
Corona-Satellite-Imagery
No pun intended, but we really discussed the images taken by the Corona/Keyhole-Satellite-Program. This is one of the earliest satellite programs dedicated to take images of more or less the entire earth originally for military purposes and espionage (see https://de.wikipedia.org/wiki/Keyhole). The images date back to the early 1960s and are of particular value to document the change of landscapes over the last 50 years. The images have a quite high ground resolution (about 10 m) and can be used for stereoscopic 3D-reconstructions. A lot of the early images are now declassified and can be obtained via the USGS (https://earthexplorer.usgs.gov/).
Benjamin showed us the results of a DAI program to collect and georeference these images for the entire north of Africa. He pointed out, that the original 70mm-films show transects of 200-300 km length. This means that apart from the perspective adjustment the georefencing has to adjust for the curvature of the earth. A tool that does that has now been developed for GRASS GIS and is available on Github (https://github.com/OSGeo/grass/tree/releasebranch_7_8). It will be part of the next release of GRASS GIS as part of the i.ortho-modules.
Finding islands in a DEM
Due to shortage of time we were not able to discuss the problem related to r.geomorphon posted earlier in May. We will do this at our next meeting.
As usual I am tinkering with ways of automatically identifying islands on vector maps and in DEMs. Currently I am mostly working on DEMs and actually have a solution for my problem, although in practice it doesn't exactly work.
Exemplary I am using a clip from the EU-DEM (Image 1) of an water rich area north of Berlin.
As one can see, there are some islands in the lakes, which should be easily identifiable due to their topographic/geomorphological properties, so I am using r.geomorphon to identify the geomorphological forms in the DEM (Image 2). The lake areas (which are plains in these SRTM data) are marked in grey, the islands appear as slopes, shoulders, ridges and peaks.
Now all I need to do is to identify all flats, find "rings" within them and identify the geomorpholocial forms in these rings. As a first step, I vectorize the output of r.geomorphon using GDAL Polygonize (Image 3).
The next logical step would be to select the plains and use e.g. the QGIS tool Difference to create polygons within the rings of these plains. Unfortunately this step does not work, because there are to many topological errors in the vectorized layer. Difference throws the following error message: "Feature (82) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option. Execution failed after 0.53 seconds". It appears, that the raster output of r.geomorphon creates a lot of single cells on the edges of larger features, which can not be vectorized properly (Image 4).
I tried several tools to clean the topological errors (e.g. v.clean), but I did not find one that works. I certainly could remove all the errors by hand, but that would be rather pointless since I could also spend the same time to identify and vectorize the islands by hand. I thought about smoothing the original DEM, but I am afraid I going to lose the islands if I do that.
I would appreciate ideas and thoughts on that procedure.
I am currently trying to reconstruct more or less convincing networks for the north-eastern part of Germany based on the EU-DEM using the tool r.walk in GRASS GIS. After showing talking about that during the last meeting, I tried out a data set of a more hilly region (around the nice mountain Brocken in central Germany).
I ran a few tests with several different start and stop points, trying to reconstruct a network between some of the larger settlements surrounding the mountain. The workflow was basically 1) r.walk with 1 of the points as the start and all the others as the targets iterated for all points and 2) r.drain on the resulting cost surfaces inverted, meaning the start point of r.walk is now the target point and vice versa. What I got is something like this.
My questions would be:
1) Is there a way to script this? I did not find a way to convince QGIS or GRASS GIS to select every point in a layer in an iterative process for the start points for r.walk and to invert this selection for the corresponding target points.
2) Since r.walk is nonisotropic, sometimes there are two paths connecting two points, according to which direction one walks. Is there a way to normalize this?
This is a brief explanation of the problem I presented in the 28. January 2020 GIS stammtisch.
I have a database in Excel. the larfer part of the data is nominal. The solumns represent categorical variables and the rows represent the sites.
I am looking for different methods that I can analyze the data. group the sites, or find correlations between variables.
I would like to find groups of sites that are more similar to eachother according to this data.
In total I have around 60 sites and 30 variables.
I created a contingency table and tried the Chi square test.
I also converterd the data into binary data and made the fisher exact test on it. but it is difficult to interpret the data.
I was thinking about doing other tests like uncertainty coefficient, etc.
I am looking for other suggestions. 
I got a point cloud representing the early medieval sites in Brandenburg (about 7000). They cluster roughly like displayed below.
There are some articles (e.g. Nakoinz 2015) about how to deduce networks from clusters like this, by generating cluster centroids and using a Delaunay-triangulation on those. However, they usually don't go into too much detail. Has anyone worked with something like this before or can give me hints on literature or methodology.
What I would like to get in the end is an approximation of the most likely spatial network, which I would like to compare to the cost surface of the area. My guess would be that network connections would be stronger along paths with lower costs.
Wir hatten das letzte mal schon über mein Anliegen gesprochen. Hier eine kurze Darstellung wie weit ich bisher bin.
Leider fand ich noch keine Zeit, mir eine ältere Version von Qgis runterzuladen und zu installieren um das Interpolationsplugin zu erhalten.
Hier eine kurze Darstellung, wie weit ich bisher bin:
Vorhanden sind die eingemessenen Punkte der Probenentnahme, in der Attributtabelle dazu die Höhe der Phosphatwerte.
Das Raster, die Quadranten ebenfalls mit den dazugehörigen Phosphatwerten in der Attributtabelle. Hier schon farbig abgestuft.
Vorläufige Kartierung des Phosphatgehaltes im Rastersystem der Probenentnahme (Software QGIS 3.4, Abstufung der Werte – Natürliche Unterbrechungen (Jenks)
Mittels des Plugins Contour konnte ich diese Linien (mit dem Punktlayer der Probenentnahme) erstellen. Kann sie aber nicht speichern.
Fragen in die Gruppe: Wie kann ich diese Konturen speichern?
Wenn ich das Interpolationsplugin habe (ich versuche es bis zum Termin): Bis zu welcher Interpolation ist es wissenschaftlich vertretbar?
Welche Abstufung der Werte ist am sinnvollsten oder am korrektesten? - Gleichmäßige Intervalle, logarithmisch oder wie oben natürliche Unterbrechung?
Viele Grüße, Susen
