Category Archives: data

The EuroLST seamless and gap-free daily European maps of land surface temperatures

The EuroLST dataset is seamless and gap-free with a temporal resolution of four records per day and enhanced spatial resolution of 250 m. This newly developed reconstruction method (Metz et al, 2014) has been applied to Europe and neighbouring countries, resulting in complete daily coverage from 2001 onwards. To our knowledge, this new reconstructed LST time series exceeds the level of detail of comparable reconstructed LST datasets by several orders of magnitude. Studies on emerging diseases, parasite risk assessment and temperature anomalies can now be performed on the continental scale, maintaining high spatial and temporal detail. In their paper, the authors provide examples for implications and applications of the new LST dataset, such as disease risk assessment, epidemiology, environmental monitoring, and temperature anomalies.

Reconstructed MODIS Land Surface Temperature Dataset, at 250m pixel resolution (click figure to enlarge):
MODIS lst time series reconstructed

Section 1. Article and data citation:

EuroLST has been produced by the former PGIS group at Fondazione Edmund Mach, DBEM based on daily MODIS LST (Product of NASA) maps.

Metz, M.; Rocchini, D.; Neteler, M. 2014: Surface temperatures at the continental scale: Tracking changes with remote sensing at unprecedented detail. Remote Sensing. 2014, 6(5): 3822-3840 (DOI | HTML | PDF)

Section 2. Used software

Open Source commands used in processing (GRASS GIS 7):
links to the related manual pages involved in the data preparation

  • i.pca: Principal Components Analysis (PCA) for image processing.
  • r.regression.multi: it calculates multiple linear regression from raster maps
  • v.surf.bspline: it performs bicubic or bilinear spline interpolation with Tykhonov regularization.

Furthermore:

  • r.bioclim: calculates various bioclimatic indices from monthly temperature and optional precipitation time series (install in GRASS GIS 7 with “g.extention r.bioclim”)
  • pyModis: Free and Open Source Python based library to work with MODIS data

Section 3. Metadata

Map projection: EPSG 3035, prj file
PROJCS["Lambert Azimuthal Equal Area",
    GEOGCS["grs80",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",52],
    PARAMETER["longitude_of_center",10],
    PARAMETER["false_easting",4321000],
    PARAMETER["false_northing",3210000],
    UNIT["Meter",1]]

1. Selected open data derived from EuroLST

Section 1. BIOCLIM derived from reconstructed MODIS LST at 250m pixel resolution

BIO1: Annual mean temperature (°C*10) BIO2: Mean diurnal range (Mean monthly (max - min tem)) BIO3: Isothermality ((bio2/bio7)*100) BIO4: Temperature seasonality (standard deviation * 100) BIO5: Maximum temperature of the warmest month (°C*10) BIO6: Minimum temperature of the coldest month (°C*10) BIO7: Temperature annual range (bio5 - bio6) (°C*10) BIO10: Mean temperature of the warmest quarter (°C*10) BIO11: Mean temperature of the coldest quarter (°C*10)

BIOCLIM-like European LST maps following the “Bioclim” definition (Hutchinson et al., 2009) – derived from 10 years of reconstructed MODIS LST (download to be completed) as GeoTIFF files, 250m pixel resolution, in EU LAEA projection:

Each ZIP file contains the respective GeoTIFF file (for cell value units, see below), the color table as separate ASCII file and a README.txt with details.

Section 2. WMS/WCS Server

Using this URL, you can read the EuroLST BIOCLIM data directly via OGC WMS and WCS protocol:

http://geodati.fmach.it/production/ows_europe_lst

Section 3. OpenData License

The data published in this page are open data and released under the ODbL (Open Database License).

The full EuroLST dataset is not released online as open data (size: 18TB), please ask Luca Delucchi or Roberto Zorer for more info


2. Acknowledgments

The MOD11A1.005, MYD11A1.005 were retrieved from the online web site, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, http://e4ftl01.cr.usgs.gov/

Importing and visualizing LAS LiDAR files in GRASS GIS 7: r.in.lidar

banner_lidar

LiDAR point cloud data are commonly delivered in the ASPRS LAS format. The format is supported by libLAS, a BSD-licensed C++ library for reading/writing these data. GRASS GIS 7 supports the LAS format directly when built against libLAS (as the case for most binary packages being available for download).

In this exercise we will import a sample LAS data set covering a tiny area close to Raleigh, NC (USA), belonging to the North Carolina sample data set. Sample LAS data download: http://grass.osgeo.org/sampledata/north_carolina/  (25MB).

For a full exercise, we will, however, assume that no GRASS GIS location is ready so far (so: newbies are welcome!) and create a new one initially.

1. 1. Having the LAS file: now what?

In the first place, check the metadata of the LAS file using the lasinfo command (comes with libLAS; here only parts of the output shown):

lasinfo points.las
---------------------------------------------------------
  Header Summary
---------------------------------------------------------
  Version:  1.2
  Source ID:  0
  Reserved:  0
  Project ID/GUID:  '00000000-0000-0000-0000-000000000000'
  System ID:  'libLAS'
  Generating Software:  'libLAS 1.2'
[...]
  Spatial Reference:
None
[...]
---------------------------------------------------------
  Point Inspection Summary
---------------------------------------------------------
  Header Point Count: 1287775
  Actual Point Count: 1287775

  Minimum and Maximum Attributes (min,max)
---------------------------------------------------------
  Min X, Y, Z:   6066629.86, 2190053.45, -3.60
  Max X, Y, Z:   6070237.92, 2193507.74, 906.00
  Bounding Box:  6066629.86, 2190053.45, 6070237.92, 2193507.74
  Time:  0.000000, 0.000000
  Return Number:  1, 3
  Return Count:  1, 7
  Flightline Edge:  0, 0
  Intensity:  0, 256
  Scan Direction Flag:  0, 0
  Scan Angle Rank:  0, 0
  Classification:  2, 7
  Point Source Id:  0, 0
  User Data:  0, 0
  Minimum Color (RGB):  0 0 0 
  Maximum Color (RGB):  0 0 0 

  Number of Points by Return
---------------------------------------------------------
  (1) 1225886  (2) 61430  (3) 459

  Number of Returns by Pulse
---------------------------------------------------------
  (1) 30877  (2) 153  (5) 1225886  (6) 30706  (7) 153

  Point Classifications
---------------------------------------------------------
  647337 Ground (2) 
  639673 Low Vegetation (3) 
  740 Building (6) 
  25 Low Point (noise) (7) 
  -------------------------------------------------------
  0 withheld
  0 keypoint
  0 synthetic
  -------------------------------------------------------

We see: no spatial reference system indicated!
Luckily we know from here that the projection is  NAD83(HARN) / North Carolina, LCC 2SP metric, EPSG code 3358). Furthermore we see:

  • Number of Points by Return: 3 (i.e., first, mid, last)
  • Point Classifications: the points are already classified as “Ground” (class 2),  “Low Vegetation” (3),  “Building” (6), and Low Point (noise) (class 7). Something to play with later.

Time to create a GRASS GIS location and import the LAS file.

2. 2. Creating a GRASS GIS location for the LAS file

Since we know the EPSG code of the projection, that’s an easy task. Please note that GRASS GIS can generate locations directly from SHAPE files (with .prj file), GeoTIFF and more.

We fire up GRASS GIS 7 and open the Location Wizard:

grass7_startup

In the Location Wizard, we first define a name for the new location:

grass7_loc_wizard1

We select the “EPSG code” method for creating a new location:

grass7_loc_wizard2

You can search for “North Carolina” and select the EPSG code 3358 from the list:

grass7_loc_wizard3

Next summary should show up as follows (be sure to have the metric projection shown!):

grass7_loc_wizard4

With “Finish” you reach this notification (indeed, nothing to change! It is already fine):

grass7_loc_wizard5

Since we want to import the LAS file, no need to manually define any region extent here – just say “No”:

grass7_loc_wizard6

While we could import the data also into the PERMANENT mapset, we prefer to create an own mapset “lasdata” for our LAS data (once you reach hundreds of maps to manage, you will be happy about the concept of mapsets):

grass7_loc_wizard7

Voilà, we get back to the initial startup screen and can now start our GRASS GIS session with our “nc_nad83_lcc” location and “lasdata” mapset within the location: “Start GRASS”!

grass7_loc_wizard8

3. 3. Import of the LAS file

When creating a new location from a GeoTIFF or SHAPE file (or other GDAL supported format), then the data set is imported right away. This is not the case for LAS files, also due to the fact that we can directly apply binning statistics during import of the LAS file (e.g. percentiles, min or max) and create a raster surface from the points right away rather than importing them as vector points.

3. a) Creating a raster surface from LAS during import

The LAS import into a raster surface is available through r.in.lidar:

grass7_las_import1

First the LAS file needs to be selected and an output file name specified (in this example, we want to extract the 95th percentile as binning method, hence a reasonable map name):

grass7_las_import2

In the “Statistic” tab, we select the “percentile” method along with 95 as value:

grass7_las_import3

In the “Optional” tab we activate to extend the computational region from the LAS file and, since the spatial reference system metadata are lacking from the LAS file, also “override dataset projection” to use that predefined in the location. Finally, we define 5m as desired raster resolution for the resulting raster map:

grass7_las_import4

Upon conpletion of the import/binning, the new raster elevation map is shown after zooming to the map (r.in.lidar -e … restores upon completion the previous region settings, hence we may have to zoom):

grass7_las_import5

Now we can start to analyze or visually explore the imported LAS file.

4. 4. Visual LiDAR data exploration

Using the wxNVIZ 3D viewer, we can easily fly through the new DEM. Switch in the Map Display to “3D view” (1). Note that the default rendering is initially done at low visual resolution for speed reasons. You can switch to “Rotate mode” as well to easily navigate with the mouse. In the “Data” tab (2) you can increase the visual resolution (3) to obtain a crisp view:

grass7_las_viz1

 

Now all kinds of analysis steps may follow.

5. Outlook

For true LiDAR processing as points, see the following GRASS GIS 7 modules: v.in.lidar (for point import), v.lidar.correction, v.lidar.edgedetection, v.lidar.growing, v.outlier, and v.surf.bspline.

Using the 25m EU-DEM for shading OpenStreetMap layers

Inspired by Vaclav Petras posting about “Did you know that you can see streets of downtown Raleigh in elevation data from NC sample dataset?” I wanted to try the new GRASS GIS 7 Addon r.shaded.pca which creates shades from various directions and combines then into RGB composites just to see what happens when using the new EU-DEM at 25m.

To warm up, I registered the “normally” shaded DEM (previously generated with gdaldem) with r.external in a GRASS GIS 7 location (EPSG 3035, LAEA) and overlayed the OpenStreetMap layer using WMS with GRASS 7’s r.in.wms. An easy task thanks to University of Heidelberg’s www.osm-wms.de. Indeed, they offer a similar shading via WMS, however, in the screenshot below you see the new EU data being used for controlling the light on our own:

OpenStreetMap shaded with EU DEM 25m

OpenStreetMap shaded with EU DEM 25m (click to enlarge)

Next item: trying r.shaded.pca… It supports multi-core calculation and the possibility to strengthen the effects through z-rescaling. In my example, I used:

r.shaded.pca input=eu_dem_25 output=eu_dem_25_shaded_pca nproc=3 zmult=50

The leads to a colorized hillshading map, again with the OSM data on top (50% transparency):

eu_dem_25m_PCA_shaded_OSM_trento_rovereto_garda_lake

OpenStreetMap shaded with r.shaded.pca using EU DEM 25m (click to enlarge)

Yes, fun – I like it 🙂

Data sources:

EU-DEM: new Digital Surface Model at 25m

eu_dem_upper_garda_lake_riva_arco_italy

EU DEM 25m upper Garda Lake area with Riva del Garda and Arco (Italy). 3D view in wxNVIZ – GRASS GIS 7

The 25m European Digital Elevation Model (EU-DEM, Version 1) is a Digital Surface Model (DSM) representing the first surface as illuminated by the sensors:

eu_dem_s_michele_rotaliana_italy

EU DEM Rotaliana with Mezzocorona and S. Michele (Italy). Produced using Copernicus data and information funded by the European Union – EU-DEM layers.

Its elevations were captured at 1 arc second postings (2.78E-4 degrees). The tiles are provided at 25m resolution in EU-LAEA (EPSG. 3035) projection, temporal coverage: 2000, published in Oct 2013. It is a realisation of the Copernicus programme, managed by the European Commission, DG Enterprise and Industry. Metadata are provided here. According to their “Methodology” page it is a hybrid product based on SRTM and ASTER GDEM data fused by a weighted averaging approach and it has been generated as a contiguous dataset divided into 1 degree by 1 degree tiles, corresponding to the SRTM naming convention. In addition to the DEM data, a colour shaded relief image over Europe is provided.

From the metadata page: “The EU-DEM data are provided as is, i.e. without a formal validation yet. An independent statistical validation is scheduled as part of the GIO land monitoring service activities, and will be made available in the course of 2014.

1. Data download

Note that the GeoTIFF files are of major size, up to 5 GB:

2. Data import

The data come as ZIP compressed files, hence unzipping occurs (or simply use the fancy “vsizip” driver in GDAL).

Hint for GRASS GIS users: instead of importing the data, you can use the r.external command to register the GeoTIFF DEM file instead of imorting it within a EU LAEA projected location.

Enjoy!

eu_dem_trento_adige_s_michele_italy

Processing Landsat 8 data in GRASS GIS 7: Import and visualization

banner_landsat_rgb

The Landsat 8 mission is a collaboration between the U.S. Geological Survey (USGS) and National Aeronautics and Space Administration (NASA) which continues the acquisition of high-quality data for observing land use and land cover change.

The Landsat 8 spacecraft which was launched in 2013 carries they following key instruments:

  • OLI: the Operational Land Imager which collects data in the visible, near infrared, and shortwave infrared wavelength regions as well as a panchromatic band. With respect to Landsat 7 two new spectral bands have been added: a deep-blue band for coastal water and aerosol studies (band 1), and a band for cirrus cloud detection (band 9). Furthermore, a Quality Assurance band (BQA) is also included to indicate the presence of terrain shadowing, data artifacts, and clouds.
  • TIRS: The Thermal Infrared Sensor continues thermal imaging and is also intended to support emerging applications such as modeling evapotranspiration for monitoring water use consumption over irrigated lands.

The data from Landsat 8 are available for download at no charge and with no user restrictions.

For our analysis example, we’ll obtain (freely – thanks to NASA and USGS!) a Landsat 8 scene from http://earthexplorer.usgs.gov/

First of all, you should register.

1. Landsat 8 download procedure

1. Enter Search Criteria:

  • path/row tab, enter Type WRS2: Path: 16, Row: 35
  • Date range: 01/01/2013 – today
  • Click on the “Data sets >>” button

2. Select Your Data Set(s):

  • Expand the entry + Landsat Archive
    [x] L8 OLI/TIRS
  • Click on the “Results >>” button

(We jump over the additional criteria)

4. Search Results

From the resulting list, we pick the data set:

earthexplorer_selection_lsat8Entity ID: LC80160352013134LGN03
Coordinates: 36.04321,-79.28696
Acquisition Date: 14-MAY-13

Using the “Download options”, you can download the data set (requires login). Select the choice:
[x] Level 1 GeoTIFF Data Product (842.4 MB)

You will receive the file “LC80160352013134LGN03.tar.gz”.

2. Unpacking the downloaded Landsat 8 dataset

To unpack the data, run (or use a graphical tool at your choice):

tar xvfz LC80160352013134LGN03.tar.gz

A series of GeoTIFF files will be extracted: LC80160352013134LGN03_B1.TIF, LC80160352013134LGN03_B2.TIF, LC80160352013134LGN03_B3.TIF, LC80160352013134LGN03_B4.TIF, LC80160352013134LGN03_B5.TIF, LC80160352013134LGN03_B6.TIF, LC80160352013134LGN03_B7.TIF, LC80160352013134LGN03_B8.TIF, LC80160352013134LGN03_B9.TIF, LC80160352013134LGN03_B10.TIF, LC80160352013134LGN03_B11.TIF, LC80160352013134LGN03_BQA.TIF

We may check the metadata with “gdalinfo“:

gdalinfo LC80160352013134LGN03_B1.TIF
Driver: GTiff/GeoTIFF
Files: LC80160352013134LGN03_B1.TIF
Size is 7531, 7331
Coordinate System is:
PROJCS["WGS 84 / UTM zone 17N",
  GEOGCS["WGS 84",
  DATUM["WGS_1984",
  SPHEROID["WGS 84",6378137,298.257223563,
...
Pixel Size = (30.000000000000000,-30.000000000000000)
...

3. Want to spatially subset the Landsat scene first?

If you prefer to cut out a smaller area (subregion), check here for gdal_translate usage examples.

4. Import into GRASS GIS 7

Note: While this Landsat 8 scene covers the area of the North Carolina (NC) sample dataset, it is delivered in UTM rather than the NC’s state plane metric projection. Hence we preprocess the data first in its original UTM projection prior to the reprojection to NC SPM.

Using the Location Wizard, we can import the dataset easily into a new location (in case you don’t have UTM17N not already created earlier):

grass70 -gui

grass7_loc_wizard1
grass7_loc_wizard2
grass7_loc_wizard3
grass7_loc_wizard4
grass7_loc_wizard5
grass7_loc_wizard6
grass7_loc_wizard7
grass7_loc_wizard8
grass7_loc_wizard9

 

 

 

Now start GRASS GIS 7 and you will find the first band already imported (the others will follow shortly!).

For the lazy folks among us, we can also create a new GRASS GIS Location right away from the dataset on command line:

grass70 -c LC80160352013134LGN03_B10.TIF ~/grassdata/utm17n

5. Importing the remaining Landsat 8 bands

The remaining bands can be easily imported with the raster import tool:

grass7_import1

The bands can now be selected easily for import:

grass7_import2

  • Select “Directory” and navigate to the right one
  • The available GeoTIFF files will be shown automatically
  • Select those you want to import
  • You may rename (double-click) the target name for each band
  • Extend the computation region accordingly automatically

Click on “Import” to get the data into the GRASS GIS location. This takes a few minutes. Close the dialog window then.

In the “Map layers” tab you can select the bands to be shown:

grass7_visualize1

6. The bands of Landsat 8

(cited from USGS)

Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) images consist of nine spectral bands with a spatial resolution of 30 meters for Bands 1 to 7 and 9. New band 1 (ultra-blue) is useful for coastal and aerosol studies. New band 9 is useful for cirrus cloud detection. The resolution for Band 8 (panchromatic) is 15 meters. Thermal bands 10 and 11 are useful in providing more accurate surface temperatures and are collected at 100 meters. Approximate scene size is 170 km north-south by 183 km east-west (106 mi by 114 mi).

Landsat 7 Wavelength (micrometers) Resolution (meters) Landsat 8 Wavelength (micrometers) Resolution (meters)
Band 1 – Coastal aerosol 0.43 – 0.45 30
Band 1 – Blue 0.45 – 0.52 30 Band 2 – Blue 0.45 – 0.51 30
Band 2 – Green 0.52 – 0.60 30 Band 3 – Green 0.53 – 0.59 30
Band 3 – Red 0.63 – 0.69 30 Band 4 – Red 0.64 – 0.67 30
Band 4 (NIR) 0.77 – 0.90 30 Band 5 – Near Infrared (NIR) 0.85 – 0.88 30
Band 5 (SWIR 1) 1.55 – 1.75 30 Band 6 – SWIR 1 1.57 – 1.65 30
Band 7 (SWIR 2) 2.09 – 2.35 30 Band 7 – SWIR 2 2.11 – 2.29 30
Band 8 – Panchromatic 0.52 – 0.90 15 Band 8 – Panchromatic 0.50 – 0.68 15
Band 9 – Cirrus 1.36- 1.38 30
Band 6 – Thermal Infrared (TIR) 10.40 -12.50 60* (30) Band 10 – Thermal Infrared (TIRS) 1 10.60 – 11.19 100* (30)
Band 11 – Thermal Infrared (TIRS) 2 11.50- 12.51 100* (30)
* ETM+ Band 6 is acquired at 60-meter resolution. Products processed after February 25, 2010 are resampled to 30-meter pixels. * TIRS bands are acquired at 100 meter resolution, but are resampled to 30 meter in delivered data product.

7. Natural color view (RGB composite)

Due to the introduction of a new “Cirrus” band (#1), the BGR bands are now 2, 3, and 4, respectively. See also “Common band combinations in RGB” for Landsat 7 or Landsat 5, and Landsat 8.

From Digital Numer (DN) to reflectance:
Before creating an RGB composite, it is important to convert the digital number data (DN) to reflectance (or optionally radiance). Otherwise the colors of a “natural” RGB composite do not look convincing but rather hazy (see background in the next screenshot). This conversion is easily done using the metadata file which is included in the data set with i.landsat.toar:

grass7_landsat_toar0
grass7_landsat_toar1
grass7_landsat_toar2
grass7_landsat_toar3

Now we are ready to create a nice RGB composite (hint 2015: i.landsat.rgb has been renamed to i.colors.enhance):

grass7_landsat_rgb0

grass7_landsat_rgb1

Select the bands to be visually combined:

grass7_visualize2

… and voilà !

grass7_landsat_rgb2

8. Applying the Landsat 8 Quality Assessment (QA) Band

One of the bands of a Landsat 8 scene is named “BQA” which contains for each pixel a decimal value representing a bit-packed combination of surface, atmosphere, and sensor conditions found during the overpass. It can be used to judge the overall usefulness of a given pixel.

We can use this information to easily eliminate e.g. cloud contaminated pixels. In short, the QA concept is (cited here from the USGS page):

Cited from http://landsat.usgs.gov/L8QualityAssessmentBand.php‎

For the single bits (0, 1, 2, and 3):
0 = No, this condition does not exist
1 = Yes, this condition exists.

The double bits (4-5, 6-7, 8-9, 10-11, 12-13, and 14-15) represent levels of confidence that a condition exists:
00 = Algorithm did not determine the status of this condition
01 = Algorithm has low confidence that this condition exists (0-33 percent confidence)
10 = Algorithm has medium confidence that this condition exists (34-66 percent confidence)
11 = Algorithm has high confidence that this condition exists (67-100 percent confidence).

Detailed bit patterns (d: double bits; s: single bits):
d – Bit 15 = 0 = cloudy
d – Bit 14 = 0 = cloudy
d – Bit 13 = 0 = not a cirrus cloud
d – Bit 12 = 0 = not a cirrus cloud
d – Bit 11 = 0 = not snow/ice
d – Bit 10 = 0 = not snow/ice
d – Bit 9 = 0 = not populated
d – Bit 8 = 0 = not populated
d – Bit 7 = 0 = not populated
d – Bit 6 = 0 = not populated
d – Bit 5 = 0 = not water
d – Bit 4 = 0 = not water
s – Bit 3 = 0 = not populated
s – Bit 2 = 0 = not terrain occluded
s – Bit 1 = 0 = not a dropped frame
s – Bit 0 = 0 = not fill

Usage example 1: Creating a mask from a bitpattern

We can create a cloud mask (bit 15+14 are set) from this pattern:
cloud: 1100000000000000

Using the Python shell tab, we can easily convert this into the corresponding decimal number for r.mapcalc:

Cited from http://landsat.usgs.gov/L8QualityAssessmentBand.php‎

Welcome to wxGUI Interactive Python Shell 0.9.8

Type "help(grass)" for more GRASS scripting related information.
Type "AddLayer()" to add raster or vector to the layer tree.

Python 2.7.5 (default, Aug 22 2013, 09:31:58) 
[GCC 4.8.1 20130603 (Red Hat 4.8.1-1)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> int(0b1100000000000000)
49152

Using this decimal value of 49152, we can create a cloud mask:

# set NULL for cloudy pixels, 1 elsewhere:
r.mapcalc "cloudmask = if(LC80160352013134LGN03_BQA == 49152, null(), 1 )"

# apply this mask
r.mask cloudmask

In our sample scene, there are only tiny clouds in the north-east, so no much to be seen. Some spurious cloud pixels are scattered over the scene, too, which could be eliminated (in case of false positives) or kept.

Usage example 2: Querying the Landsat 8 BQA map and retrieve the bitpattern

Perhaps you prefer to query the BQA map itself (overlay the previously created RGB composite and query the BSA map by selecting it in the Layer Manager). In our example, we query the BQA value of the cloud:

Using again the Python shell tab, we can easily convert the decimal number (used for r.mapcalc) into the corresponding binary representation to verify with the table values above.

>>> x=61440
>>> print(bin(x & 0xffffffff))
0b1111000000000000

Hence, bits 15,14,13, and 12 are set: cloudy and not a cirrus cloud. Looking at the RGB composite we tend to agree 🙂 Time to mask out the cloud!

wxGUI menu >> Raster >> Mask [r.mask]

Or use the command line, as shown already above:

# remove existing mask (if active)
r.mask -r

# set NULL for cloudy pixels, 1 elsewhere:
r.mapcalc "cloudmask = if(LC80160352013134LGN03_BQA == 61440, null(), 1 )" --o

# apply the new mask
r.mask cloudmask

The visual effect in the RGB composite is minimal since the cloud is white anyway (as NULL cells, too). However, it is relevant for real calculations such as NDVI (vegetation index) or thermal maps.

We observe dark pixels around the cloud originating from thin clouds. In a subsequent identification/mask step we may eliminate also those pixels with a subsequent filter.

See also Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening

Landsat free data downloads

In a recent USGS press release I found:

“Although the USGS does not have detailed records since the mission’s inception in 1972, there is good evidence that more data have been distributed in the last 6 months than in the entire first 36 years of the Landsat missions combined.”

Amazing, no? Don’t have to say more when blinking at EU data sources…


Thermal map of Valsugana near Trento, Italy
(30 July 2003, ~ 9:30)

ASTER GDEM 30m quality assessment

The recent publication of the ASTER GDEM (30m) is a great step towards a worldwide high resolution elevation model. I have done some test in the Southern Alps around Trento which is among the most complex terrains in Europe. The scope was to calculate a difference map to the local high resolution DEM.

The steps

1. Mosaicking and reprojection to UTM32/WGS84 of ASTER GDEM tiles of the area (using gdalwarp).
2. Import of the ASTER GDEM map into GRASS (r.in.gdal)
3. Creation of a difference map to the provincial DEM (r.mapcalc, r.colors)
4. Histogram and univariate statistics (d.histogram, r.univar)

The results
(click for higher resolution or download slides as PDF):

Comments

The hydrological pattern appears to coincide with the provincial rivers map. It can be seen that lakes weren’t masked during the ASTER DEM preparation, howver, these could be identified for many areas in the world with OpenStreetMap layers.
Despite spikes (including unexpected craters), the overall quality appears to be acceptable for this first version of ASTER GDEM – the standard deviation is 18m for the test area with outliers predominantely found in the areas of complex terrain rather than in the valley floors. A shift of -2.3m is also observed (reprojection artifact?).
It will be worth to check if SRTM data (90m) could help to identify and remove the spikes from the DEM. More hopefully in the near future…

Acknowledgments: ASTER GDEM is a product of METI and NASA.

OGR vector data tips and tricks

Get vector map extent


You can easily grep the map extent of a vector map (bounding box):

ogrinfo /cdrom/ITALY_GC.SHP ITALY_GC | grep Extent
Extent: (9.299460, 43.787362) – (13.911350, 47.070188)

Merge of two SHAPE files

Merge of two SHAPE files ‘file1.shp’ and ‘file2.shp’ into a new file ‘file_merged.shp’ is performed like this:

ogr2ogr file_merged.shp file1.shp
ogr2ogr -update -append file_merged.shp file2.shp -nln file_merged file2

The second command is opening file_merged.shp in update mode, and trying to find existing layers and append the features being copied. The -nln option sets the name of the layer to be copied to.

Vector map reprojection

We reproject from the source projection (as defined in .prj file) to WGS84/LL:

ogr2ogr vmap0rd_ll.shp -t_srs “EPSG:4326” vmap0rd.shp

If the .prj file is missing, you can use the ‘epsg_tr.py’ utility to create it if you know the EPSG code:

epsg_tr.py -wkt 4326 > cities.prj

Reproject to current GRASS location projection:

ogr2ogr -t_srs “`g.proj -wf`” polbnda_italy_GB_ovest.shp polbnda_italy_LL.shp

Cut out a piece of a vector map
Use spatial query extents: -spat xmin ymin xmax ymax (W S E N)

ogr2ogr ARC_BZ.shp -spat 10 45 13 47 ARC.shp

Get VMAP0 metadata info:

ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr | grep bnd
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’
ogrinfo -ro gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘roadl@trans(*)_line’

MAP0: Extract spatial subregion, reproject from NAD83 to WGS84

# coordinate order: W S E N
ogr2ogr -spat 19.95035 -26.94755 29.42989 -17.72624 -t_srs ‘EPSG:4326’ \
polbnda_botswana.shp gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \
‘polbnda@bnd(*)_area’

OGR and SQL

Sample ‘where’ statements (use -sql for PostgreSQL driver):

# -where ‘fac_id in (195,196)’
# -where ‘fac_id = 195’
ogrinfo -ro -where ‘fac_id in (195,196)’ \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

VMAP0 examples

Find out the Countries VMAP0 coding:

ogdi_info -u gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr \
-l ‘polbnda@bnd(*)_area’ -f area | grep Botswana

or read the VMAP0 Military specs, page 75

Extract Botswana, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘BC'” polbnda_botswana.shp \
gltp:/vrf/grass0/warmerdam/v0soa/vmaplv0/soamafr ‘polbnda@bnd(*)_area’

Extract Germany, reproject on the fly from NAD83 to WGS84, store as SHAPE:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘GM'” polbnda_germany.shp
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘polbnda@bnd(*)_area’
ogrinfo -summary polbnda_germany.shp polbnda_germany | grep Extent
# Extent: (5.865639, 47.275776) – (15.039889, 55.055637)
# W S E N

VMAP0 Contour lines for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
contour_lines.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘contourl@elev(*)_line’

VMAP0 elevation spots (points) for Germany:

ogr2ogr -t_srs “EPSG:4326” -spat 5.865639 47.275776 15.039889 55.055637 \
elevation_spots.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘elevp@elev(*)_point’

VMAP0 lakes of Trentino province in Italy:

ogr2ogr -t_srs “EPSG:4326” -where “na2 = ‘IT'” \
-spat 10.340029 45.261888 10.98727 45.98993 \
lakes_italy.shp \
gltp:/vrf/grass0/warmerdam/v0eur/vmaplv0/eurnasia ‘inwatera@hydro(*)_area’

Connect OGR and PostgreSQL/PostGIS

ogrinfo PG:’host=grass.itc.it user=postgres dbname=ogc_simple’
ogr2ogr out.shape PG:’host=grass.itc.it user=postgres dbname=ogc_simple’ lake_geom

GRASS 6 and OGR

Convert GRASS 6 vector map to SHAPE (needs GDAL-OGR-GRASS plugin):

# -nln is “new layer name” for the result:
ogr2ogr archsites.shp grassdata/spearfish60/PERMANENT/vector/archsites/head 1 \
-nln archsites

Using WKT files with ogr2ogr

The definition is in ESRI WKT format. If you save it to a text file called out.wkt you can do the following in a translation to reproject input latlong points to this coordinate system:

ogr2ogr -s_srs WGS84 -t_srs ESRI::out.wkt out_dir indatasource

Most comand line options for GDAL/OGR tools that accept a coordinate system will allow you to give the name of a file containing WKT. And if you prefix the filename with ESRI:: the library will interprete the WKT as being ESRI WKT and convert to “standard” format accordingly. The -s_srs switch is assigning a source coordinate system to your input data (in case it didn’t have this properly defined already), and the -t_srs is defining a target coordinate system to reproject to.

TIGER files in OGR

# linear features:
ogr2ogr tiger_lines.shp tgr46081.rt1 CompleteChain

# area features:
export PYTHONPATH=/usr/local/lib/python2.5/site-packages
tigerpoly.py tgr46081.rt1 tiger_area.shp

OGR CSV driver: easily indicate column types

You can now write a little csv help file to indicate the columns types to OGR. It works as follows. Suppose you have a foobar.csv file that looks like this:

“ID”,”X”,”Y”,”AREA”,”NAME”
“1”,”1023.5″,”243.56″,”675″,”FOOBAR”

Now write a foobar.csvt file like this one:

“Integer”,”Real”,”Real”,”Integer”,”String”

The driver will then use the types you specified for the csv columns. The types recognized are Integer, Real and String, DateTime, and Date.

Convert KML to CSV (WKT)

First find layers:

ogrinfo -so myfile.kml

Then convert KML to CSV:

ogr2ogr -f CSV out.csv myfile.kml -sql “select *,OGR_GEOM_WKT from myfilelayer”
cat out.csv

Or use the cool online converter: http://geoconverter.hsr.ch

GDAL raster data tips and tricks


Reading GRASS data through GDAL/OGR support

Example 1: We write out a GRASS raster map to GeoTIFF — this format
includes the coordinates within the file’s metadata:

gdal_translate -of Gtiff /usr/local/share/grassdata/spearfish/PERMANENT/cellhd/soils soilmap.tif

ogr2ogr roadsmap.shp /usr/local/share/grassdata/spearfish/PERMANENT/vector/roads/head

Fast image display with tiling
If you want fast access you might want to try converting e.g. a BIL files to a tiled TIFF, and build overviews. You can build overviews for BIL too, but it can’t be directly tiled:

# add -co “PROFILE=BASELINE” for TIF/TFW
gdal_translate source_bil global30.tif -co “TILED=YES” -co “TFW=YES” -co “PROFILE=BASELINE”
gdaladdo global30.tif 2 4 8 16

GDAL performance problem?
GDAL_CACHEMAX is normally a number of megabytes (default is 10 or so). So something like:
gdal_translate -of GTIFF -co TILED=YES –config GDAL_CACHEMAX 120 madison_1f_01.jpg madison_1f_01.tif
would use a 120MB cache.

GDAL and 1 bit maps
With a trick you can get those:
gdal_merge.py -co NBITS=1 -o dst.tif src.tif

Generate 8 bit maps for Mapserver
gdal_translate -scale in.tif out.tif
Note: As of MapServer 4.4 support has been added for classifying non-8bit raster inputs

Greyscale conversion
A “proper” conversion would involve a colorspace transformation on the RGB image into IHS or something like that, and then taking the intensity. GRASS can do things like that.

Generate an OGC WKT (SRS)
In WKT the ellipsoid is described by two parameters: the semi-major axis and the inverse flattening. For a sphere the flattening is 0 and so the inverse flattening is infinity.

# in the GDAL source code:
cd apps
make testepsg

./testepsg ‘+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m’
Validate Succeeds.
WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
PROJCS[“unnamed”,
GEOGCS[“GRS 1980(IUGG, 1980)”,
DATUM[“unknown”,
SPHEROID[“GRS80”,6378137,298.257222101]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]],
PROJECTION[“Lambert_Conformal_Conic_2SP”],
PARAMETER[“standard_parallel_1”,35],
PARAMETER[“standard_parallel_2”,65],
PARAMETER[“latitude_of_origin”,52],
PARAMETER[“central_meridian”,10],
PARAMETER[“false_easting”,4000000],
PARAMETER[“false_northing”,2800000],
UNIT[“Meter”,1]]

Simplified WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
PROJCS[“unnamed”,
GEOGCS[“GRS 1980(IUGG, 1980)”,
DATUM[“unknown”,
[..]

Extracting spatial subset (subregion)
W N E S
gdal_translate -of GTiff -projwin 636861 5152686 745617 5054047.5 \
p192r28_5t19920809_nn1.tif test1_utm.tif

Fixing broken projection/datum info for raster data
gdal_translate -of HFA -a_srs epsg:32735 /cdrom/173072lsat.img \
173072lsat_fixed.img

# or, using a WKT file
gdal_translate -of HFA -a_srs file.prj /cdrom/173072lsat.img \
173072lsat_fixed.img

Merge various import maps, re-project on the fly and extract spatial subset according to current GRASS region
eval `g.region -g`
gdalwarp -te $w $s $e $n *.TIF \
srtm_cgiar3_italy_north_LL.tif

Export to (limited) TIFF readers such as ArcView, or ImageMagick
Many tools have trouble reading multi-band TIFFs with “band interleaving”, the GDAL output default. Best is to use the INTERLEAVE=PIXEL creation option. Just add to the gdal_translate command line:
-co INTERLEAVE=PIXEL

Inserting metadata (metadata tags)
gdal_translate -outsize 37.5% 37.5% \
-mo TIFFTAG_XRESOLUTION=300 -mo TIFFTAG_YRESOLUTION=300 \
in.tif out.tif

Raster map reprojection (warping)
gdalwarp -t_srs ‘+init=epsg:26591 +towgs84=-225,-65,9’ test1.tif \
test1_gb.tif

Raster map reprojection (warping) maintaining NULL values (sea etc):

gdalwarp -r bilinear -tr 1000 1000 \
-srcnodata “-32768” -dstnodata “-32768” \
-wo “INIT_DEST=-32768” \

-t_srs epsg:32632 italy_LL.tif italy_UTM32.tif

Reprojecting external map to current GRASS location externally
gdalwarp -t_srs “`g.proj -wf`” aster.tif aster_tmerc.tif

Cut out region of interest with gdalwarp (in target coords)
Add to command line (insert values instead of letters of course:
#damn order, differs from -projwin!!
-te W S E N

Merging many small adjacent DEMs into one big map (A)
This needs GDAL compiled with Python and numpy installed:
# if not installed in standard site-packages directory
export PYTHONPATH=/usr/local/lib/python2.5/site-packages
gdal_merge.py -v -o spearfishdem.tif -n “-32768” d*.tif

Merging many small adjacent DEMs into one big map (B)
Even easier, just use gdalwarp:
gdalwarp C_1mX1m/dtm*.tif big.tif
Or just a few tiles:
gdalwarp C_1mX1m/dtm0010[4-5]* big_selection.tif

Merge various map/bands into a RGB composite
gdal_merge.py -of HFA -separate band1.img band2.img band3.img -o out.img

GDAL gdalwarp interpolation comments
Which method -rn, rb, -rc or -rcs should one use for DEM and which for data like e.g. Landsat TM reprojecting?

-tps: Enable use of thin plate spline transformer based on available GCPs.
-rn: Use nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
-rb: Use bilinear resampling.
-rc: Use cubic resampling.
-rcs: Use cubic spline resampling (slowest algorithm).

FrankW suggests:
I would suggest -rb for DEMs, and one of the cubic kernels for landsat data. Of course, there are various factors that you should take into account. Using -rb (bilinear) for the DEM will perform local averaging of the nearby pixel values in the source. This give reasonable results without introducing any risky “overshoot” effects you might see with cubic that could be disturbing for analysis or visualization in a DEM. The cubic should in theory do better at preserving edges and general visual crispness than using bilinar or nearest neighbour. However, if you are wanting to do analysis with the landsat (such as multispectral classification) I would suggest just using -rn (nearest neighbour) so as to avoid causing odd effects to the spectral values.
Nobody can’t tell you what method should be used in your case. Generally speaking, in the case of upsampling spline and cubic interpolators are more suitable (-rcs and -rc). In the case of downsampling and the same resolution it is completely up to you what method looks better. Just try them all and select the one which is most appropriate for you.

Geocoding with ‘gdal_translate’
FrankW suggests:
As far as I know there is not on-screen method for doing this, but it certainly isn’t too difficult with a little bit of semi-manual work. Open two OpenEV views, one with the unreferenced image, one with the geo-reference base you want to use. Move your cursor over the non-referenced one (let’s call it image1), record (read: write down!) the pixel x,y values. Then look at the same location in image2. Write down the geocoordinate for the pixel. You should have four numbers for each location you want to pin the image to. And so on and so on. Then use gdal_translate to translate image1.tif to image1_georefd.tif but adding the -GCP parameter for each set of coordinates. Like so…

gdal_translate -gcp 1 1 500000 5000000 \
-gcp 200 400 550000 5250000 image1.tif \
image1_geo.tif

Reading HDF ASTER
gdalinfo pg-PR1B0000-2002031402_100_001

To select a channel and warp to UTM (or whatever is inside):

gdalwarp HDF4_SDS:ASTER_L1B:”pg-PR1B0000-2002031402_100_001″:2 aster_2.tif
gdalinfo aster_2.tif