Merge branch 'development-branch' into 'master'

Make more robust to setting location if location already exists

See merge request arsf/arsf_dem_scripts!26

arsf_dem/dem_utilities.py

@@ -237,13 +237,17 @@ def offset_null_fill_dem(in_demfile, out_demfile=None,
         # Fill Null values
         print('Filling Null values')
         null_filled_name = 'patched_elevated_filled'
-        grass.run_command('r.fillnulls',
-                          input=elevated_name,
-                          output=null_filled_name,
-                          tension=40,
-                          smooth=0.1,
-                          overwrite=True)
-
+        try:
+            grass.run_command('r.fillnulls',
+                              input=elevated_name,
+                              output=null_filled_name,
+                              tension=40,
+                              smooth=0.1,
+                              overwrite=True)
+        # If this fails, pass. Will check for file in following step and print
+        # warning there.
+        except Exception as err:
+            pass
         # Check file exists (to confirm command has run correctly
         if not grass_library.checkFileExists(null_filled_name):
             dem_common_functions.WARNING('Could not NULL fill DEM, possibly there are no NULL values to fill')
@@ -270,6 +274,7 @@ def offset_null_fill_dem(in_demfile, out_demfile=None,
                           input=smoothed_name,
                           output=out_demfile,
                           nodata=nodata,
+                          flags='fc',
                           overwrite=True)
         remove_gdal_aux_file(out_demfile)
 
@@ -420,7 +425,7 @@ def patch_files(in_file_list,
                      input=patched_name,
                      output=out_file,
                      nodata=nodata,
-                     flags='f',
+                     flags='fc',
                      overwrite=True)
         remove_gdal_aux_file(out_file)
 
@@ -536,6 +541,7 @@ def replace_nodata_val(in_demfile, out_demfile=None,
                           input=replace_nodata_name,
                           output=out_demfile,
                           nodata=outnodata,
+                          flags='fc',
                           overwrite=True)
         remove_gdal_aux_file(out_demfile)
 

arsf_dem/grass_library.py

@@ -624,10 +624,12 @@ def getGRASSProjFromGDAL(in_file):
              (projection.lower() == "transverse mercator")) and \
              (spheroid.lower().find("airy") > -1):
         grass_proj = "UKBNG"
+    elif projection == "" or projection is None:
+        grass_proj = "UNPROJECTED"
     else:
         gdaldataset = None
         raise Exception('Could not identify projection from file.\n' +
-                       'WKT string is: {}'.format(spatial_ref.ExportToPrettyWkt()))
+                        'WKT string is: {}'.format(spatial_ref.ExportToPrettyWkt()))
 
     gdaldataset = None
     return grass_proj
@@ -644,12 +646,17 @@ def setLocation(projection):
 
        Returns: The location switched to, or created
     """
+    print("Setting location to: {}".format(projection))
     if projection == "WGS84LL":
         location = projection
     elif projection == "UKBNG":
         location = projection
     elif 'UTM' in projection:
         location = newLocation(projection)
+    elif projection == "UNPROJECTED":
+        # Create a dummy coordinate system for unprojected data.
+        location = newLocation("+proj=tmerc +lat_0=0+lon_0=0 +ellps=WGS84 +units=m +no_defs",
+                               projection)
     else:
         # If the projection is not the expected GRASS location format
         # try as a proj4 string.
@@ -660,11 +667,11 @@ def setLocation(projection):
             location = None
     if location is not None:
         grass.run_command('g.gisenv',
-                    set="LOCATION_NAME=%s" % (location))
+                          set="LOCATION_NAME=%s" % (location))
     return location
 #end function
 
-def newLocation(projection):
+def newLocation(projection, grass_loc=None):
     """
     Function newLocation
 
@@ -679,15 +686,24 @@ def newLocation(projection):
     * projection: projection to create new location can be
                          proj4 string or GRASS style projection name
                          (e.g., UKBNG)
+    * location_name: name for projection, if not provided will get from projection
 
     Returns:
 
     * GRASS style projection name for use with g.gisenv
 
     """
     try:
-        grass_loc = proj4_to_grass_location(projection)
-        proj4_str = projection
+        if grass_loc is not None:
+            try:
+                proj4_to_grass_location(projection)
+                proj4_str = projection
+            except Exception as err:
+                raise Exception("GRASS location name provided but projection "
+                                "not provided as proj4 string")
+        else:
+            grass_loc = proj4_to_grass_location(projection)
+            proj4_str = projection
     except Exception:
         grass_loc = projection
         proj4_str = grass_location_to_proj4(grass_loc)
@@ -920,13 +936,8 @@ def proj4_to_grass_location(in_proj4):
              (spheroid.lower().find("airy") > -1):
         grass_proj = "UKBNG"
     else:
-        dem_common_functions.WARNING('''Could not determine GRASS location name from:
-
-  {}
-
-  Setting to 'UNKNOWN'
-  '''.format(in_proj4))
-        grass_proj = "UNKNOWN"
+        raise Exception('Could not determine GRASS location name from: \n'
+                        '{}'.format(in_proj4))
 
     return grass_proj
 
@@ -1347,16 +1358,9 @@ def locationFromFile(filename, locationname=None):
     Create a new location that matches the given file. Can specify the new location name.
     """
     if locationname is None:
-        try:
-            locationname = getGRASSProjFromGDAL(filename)
-        except Exception:
-            pid = os.getpid()
-            t = time.strftime("%H%M%S")
-            locationname="custom-%s-%s" % (pid,t)
+        locationname = getGRASSProjFromGDAL(filename)
 
-    if grass.run_command('g.proj',georef=filename,flags='c',location=locationname) != 0:
-        #An error occurred
-        locationname=None
+    setLocation(locationname)
 
     return locationname
 

doc/source/tutorial_lidar.md

@@ -80,8 +80,9 @@ create_dem_from_lidar.py --in_projection UTM33N \
 
 This will create a DSM from files 1 and 2.
 Note to use all lines in the folder just provide the directory (las1.0) rather than individual lines. It is recommended you only use two for the tutorial so it runs through quicker.
-The flag '--in\_projection' is used to specify the projection of the LAS files is
-UTM33N. The resolution can be specified using `--resolution 1`, if this is not supplied the default resolution (2 m) will be used. Note that the command options (starting with --) do not need to be written in any specific order.
+The flag `--in_projection` is used to specify the projection of the LAS files, which is UTM33N.
+The resolution can be specified using `--resolution <value>`, where you replace `<value>` with the required resolution in metres. If this is not supplied the default resolution (2 m) will be used
+Note that the command options (starting with `--`) do not need to be written in any specific order.
 
 As part of the process of creating a DSM, only points which are the
 first return are kept and any points flagged as noise (class 7) are
@@ -164,7 +165,7 @@ Processing Library (APL) to geocorrect hyperspectral data, some extra
 consideration are needed:
 
 * The DSM needs to use WGS-84 Lat/Long projection and heights need to be relative to the WGS-84 ellipsoid.
-* Areas of no-data need to be filled (e.g., with a courser resolution DEM).
+* Areas of no-data need to be filled (e.g., with a coarser resolution DEM).
 * The format needs to be ENVI Band Interleaved by Line (BIL) or Band Sequential (BSQ).
 
 The same `create_dem_from_lidar.py` script can be used to generate a DSM
@@ -181,11 +182,11 @@ create_dem_from_lidar.py --in_projection UTM33N \
 ```
 As previously only two lines are specified to reduce processing time.
 
-This will create a DSM mosaic from the LAS files reproject to WGS84 Lat/Long and patch with 'EUFAR11\_02-2011-187\_ASTER.dem' (as provided with the NERC-ARF hyperspectral delivery), cropped to the bounding box of all LiDAR data plus a buffer of 2 km. This assumes the vertical datum of the DEM mosaic is the same as that required for the output projection.
+This will create a DSM mosaic from the LAS files, reproject to WGS84 Lat/Long and patch with 'EUFAR11\_02-2011-187\_ASTER.dem' (as provided with the NERC-ARF hyperspectral delivery), cropped to the bounding box of all LiDAR data plus a buffer of 2 km. This assumes the vertical datum of the DEM mosaic is the same as that required for the output projection.
 
 ## Create DSM / DTM using additional programs ##
 
-So far only a simple DSM has been created by taking the average of all first-return points within a pixel. Pixels which do not contain a return are left as no-data. There are more advanced methods of interpolation in GRASS (which we will come onto later), there are also other programs which can be used to produce a DSM. These can be accessed by two utility programs, described below. For this section the Montrose Bay data will be used, navigate to the directory the data are stored using:
+So far only a simple DSM has been created by taking the average of all first-return points within a pixel. Pixels which do not contain a return are left as no-data. There are more advanced methods of interpolation in GRASS (which we will come onto later), there are also other programs which can be used to produce a DSM. These can be accessed by two utility programs, described below. For this section the Montrose Bay data will be used, navigate to the directory the data are stored in using:
 ```
 cd ~/nerc-arf-workshop/lidar_practical/GB13_08-217
 ```
@@ -214,8 +215,7 @@ las_to_dsm.py -o LDR-GB13_08-2014-217-02_subset_dsm_grass.tif \
 
 The format of the output file is set using the extension, using '.tif' will create a GeoTIFF.
 
-Other programs such as [LAStools](http://rapidlasso.com/lastools/) [^12], [SPDLib](http://spdlib.org) [^13], [FUSION](http://forsys.cfr.washington.edu/fusion/fusion_overview.html) [^14] and [points2grid](https://github.com/CRREL/points2grid) [^15] can be used by
-if they have been installed by setting the --method flag. Note, unlike the other methods, the tools required to create a DEM in LAStools are not free and binaries are only provided for Windows (although work on Linux through Wine). It is possible to run without a license for non-profit and personal work but they will introduce noise and artifacts (such as black diagonal lines) in the data. For more details see [LAStools License](http://www.cs.unc.edu/~isenburg/lastools/LICENSE.txt) [16^]
+Other programs such as [LAStools](http://rapidlasso.com/lastools/) [^12], [SPDLib](http://spdlib.org) [^13], [FUSION](http://forsys.cfr.washington.edu/fusion/fusion_overview.html) [^14] and [points2grid](https://github.com/CRREL/points2grid) [^15] can be used, if they have been installed, by setting the --method flag. Note, unlike the other methods, the tools required to create a DEM in LAStools are not free and binaries are only provided for Windows (although work on Linux through Wine). It is possible to run without a license for non-profit and personal work but they will introduce noise and artifacts (such as black diagonal lines) in the data. For more details see [LAStools License](http://www.cs.unc.edu/~isenburg/lastools/LICENSE.txt) [16^]
 
 If you have installed SPDLib or the non-free LAStools (running through wine on Linux) try running the same command but setting `--method SPDLib` or `--method LAStools`. For example:
 
@@ -228,7 +228,7 @@ las_to_dsm.py -o LDR-GB13_08-2014-217-02_subset_dsm_spdlib.tif \
 
 ### Digital Terrain Model (DTM) ###
 
-In addition to producing a DSM another common product from LiDAR data is a Digital Terrain Model (DTM), this represents the 'bare-earth', that is the elevation excluding buildings and vegetation. Producing a DTM is more complicated than a DSM, a first order approximation is to take only last returns, which assumes multiple returns are received from a pulse the last is from the ground. However, many last returns are not from the ground for example dense vegetation where the last return does not reach the ground or buildings where there are not multiple returns. Producing an accurate DTM first requires ground returns to be selected from the point cloud and gaps where there are no ground returns (e.g., buildings) to be interpolated.
+In addition to producing a DSM another common product from LiDAR data is a Digital Terrain Model (DTM), this represents the 'bare-earth', that is the elevation excluding buildings and vegetation. Producing a DTM is more complicated than a DSM, a first order approximation is to take only last returns, which assumes multiple returns are received from a pulse and that the last is from the ground. However, many last returns are not from the ground for example dense vegetation where the last return does not reach the ground or buildings where there are not multiple returns. Producing an accurate DTM first requires ground returns to be selected from the point cloud and gaps where there are no ground returns (e.g., buildings) to be interpolated.
 
 Packages such as [LAStools](http://rapidlasso.com/lastools/) [^12], [SPDLib](http://spdlib.org) [^13] and [FUSION](http://forsys.cfr.washington.edu/fusion/fusion_overview.html) [^14] have more advanced methods for classifying ground returns (see their
 respective manuals for more details).
@@ -246,7 +246,7 @@ As shown earlier the `gdaldem` command can be used to produce hillshade images f
 
 ![Comparison of Digital Terrain Model (DTM) and Digital Surface Model (DSM)](figures/dtm_dsm.png)
 
-Note, depending on the cover the default classification and interpolation parameters for each program used by las\_to\_dtm.py may not provide the best possible results.
+Depending on the cover, the default classification and interpolation parameters for each program used by las\_to\_dtm.py may not provide the best possible results.
 In these cases it is recommended you access the programs directly, as this will provide more control over the available options.
 
 ## More Advanced Manipulation in GRASS ##