Tuesday, 16 February 2016

Overlaying Raster and Vector datasets with different projection parameters in GRASS GIS


GRASS GIS is an extremely powerful application for manipulating and modelling large raster datasets such as Landsat 8 tiff bands. As GRASS users we are aware the application does not implement 'on-the-fly-projections' to ease the pain when adding data in a different projection in the way QGIS does.

The grass developers insist doing so would introduce artifacts and some distortions when taking measurements off the overlayed data, so we have to deal with a 'one location = one projection/zone/datum/ellipsoid' combination. 

How do we get our raster images to underlay or overlay our vectors? The golden rule when using GRASS is if in doubt create a new 'location' any time you want to import or use data in a different projection. 

I mainly work with just two projections, EPSG:27700 (OSGB36) and EPSG:4326 (WGS 84) and In my case I would define two locations and associated mapsets, one for my everyday working projects in (OSGB36 here in the UK) and one for my working landsat data in (WGS84). 

 can easily change between these two locations any time and re-project (r.proj) data from one the other. For example I would simply open the (OSGB36) location and use (r.proj) module to overlay data from the (WGS84) location. Lets look at how this is done. 

1. Let's open a GRASS GIS session and create our first new 'location' and 'mapset' for our (OSGB36) location using the 'Grass Location Wizard'. This will be my everyday working project location here in the UK.

When the 'Splash' screen opens click on the left 'New' location button
Select first (default) option then Click 'Next'
Use the magnifier to find your code if you know it.

Select the option 'Used in Great Britain'
That's it... the 'Summary' confirms success.
Ignore the fact this is showing 'OSGB362' yours will show 'OSGB36'
2. Import your vector data into it with 'v.in.ogr' in the normal way.

I'm only going to import one layer from my 'sqlite' database stored inside my 'PERMANENT' directory.


Zoom to the full extent of your layer simply to check it is viewable in the projection.
3. Create a second 'location' and 'mapset' (projection code = 'utm' datum code = 'wgs84') the same as in part 1. above.
4. Import all your landsat bands (into WGS84) using 'v.in.ogr' specifying 'override projection use locations projection' this will be your working landsat project location. 

Again Zoom to the full extent of your layer to check it is viewable in this projection.
5. Change back to the 'OSGB36' location and mapset and display your vector data in a display window. Go into 'settings' > 'GRASS working environment' > 'Change Location & mapset' you will then see this:-

Change the 'location' to OSGB36 and 'mapset' to 'PERMANENT'
  6. 'Zoom out' the display so your data appears no bigger than a small blob (see below) in the middle of the display (this is to allow for the much larger extent of the Landsat tile to fit within the current location parameters of the vector mapset). 

Click on the '-' magnifier about 4/5 times.
7. Set 'computational region extent from display' option (under the sixth magnifier from the right) button in the display entitled 'various zoom options') this sets the new extent created in 6. above. 

8. 'r.proj' your Landsat image/s you wish to overlap enter text as below and then hit 'RUN'.  
Source Tab is 'WGS84' and 'PERMANENT' and the 'Input' is the raster you wish to load across into the 'OSGB36' projection.
Target Tab - leave output raster blank, leave interpolation as 'nearest' and add the resolution to 15m (this may be 30m) depending on what yours is.
9. Both layers will now overlap.

Correctly overlapping vector on raster image.
10. Set 'g.region' to the imported Landsat 8 image. 

Important note:-

If we import (instead of project) any new raw landsat images through 'v.in.ogr' whilst inside the 'OSGB36' location the new images will misalign vertically because these images are still hardcoded to their 'WGS84' projection/datum parameters even though we are importing them into a 'OSGB36' projection.

Notice the vertical misaligned vector data in the bottom of the window because i have 'imported' the landsat image (WGS84) (top) into 'OSGB36' location.
This can cause a lot of frustration for new users who may be used to 'on-the-fly-projections'.

THE GOLDEN RULE: Small 'extents' fit into bigger 'extents', Bigger 'extents' will not fit into smaller 'extents'.

Further reading:-

https://grasswiki.osgeo.org/wiki/Location_and_Mapsets

https://grasswiki.osgeo.org/wiki/GRASS_Location_Wizard

https://grasswiki.osgeo.org/wiki/Map_Reprojection

------------------------------------------------------------------------

If you wish to check your grass environment files (those which control your Grass session in the 'PERMANENT' directory) you can use the following as a guide for checking. In other words the contents of your 'WIND', 'DEFAULT_WIND', 'PROJ_EPSG', 'PROJ_INFO and 'PROJ_UNITS' files if using Landsat images.

Checking your settings (optional) for 'OSGB36' projection.

type 'g.gisenv -n' in the command console.

MAPSET=PERMANENT
GISDBASE=/home/paul/Dropbox
LOCATION_NAME=OSGB36
GUI=wxpython

type 'g.region -p' in the command console.

projection: 99 (Transverse Mercator)
zone: 0
datum: osgb36
ellipsoid: airy
north: 324535
south: 77919
west: 128407
east: 371768
nsres: 29.99829704
ewres: 30.00012327
rows: 8221
cols: 8112
cells: 66688752

WIND  (file)

proj: 99
zone: 0 
north: 324535 
south: 77919 
east: 371768 
west: 128407 
cols: 8112 
rows: 8221 
e-w resol: 30.00012327 
n-s resol: 29.99829704 
top: 1.000000000000000 
bottom: 0.000000000000000 
cols3: 8112 
rows3: 8221 
depths: 1 
e-w resol3: 30.00012327 
 n-s resol3: 29.99829704 
t-b resol: 1 

DEFAULT_WIND (file)

proj: 99 
zone: 0 
north: 1 
south: 0 
east: 1 
west: 0 
cols: 1 
rows: 1 
e-w resol: 1 
n-s resol: 1 
top: 1.000000000000000 
bottom: 0.000000000000000 
cols3: 1 
rows3: 1 
depths: 1 
e-w resol3: 1 
n-s resol3: 1 
t-b resol: 1 

PROJ_EPSG 
epsg: 27700 

PROJ_INFO (file)
name: Transverse Mercator proj: tmerc datum: osgb36 ellps: airy lat_0: 49 lon_0: -2 k: 0.9996012717 x_0: 400000 y_0: -100000 no_defs: defined towgs84: 368.000,-120.000,425.000 

PROJ_UNITS 
unit: meter units: meters meters: 1

Checking your settings (optional) for 'WGS84' projection.

Type 'g.gisenv -n' in the command console 

MAPSET=PERMANENT 
GISDBASE=/home/paul/Dropbox
LOCATION_NAME=WGS84 
GUI=wxpython 

Type 'g.region -p' in the command console.

projection: 1 (UTM) 
zone: 30 
datum: wgs84 
ellipsoid: wgs84 
north: 5848815 
south: 5605485 
west: 299085 
east: 539115 
nsres: 30 
ewres: 30 
rows: 8111 
cols: 8001 
cells: 64896111 

WIND (file)

proj: 1 
zone: 30 
north: 5848815 
south: 5605485 
east: 539115 
west: 299085 
cols: 8001 
rows: 8111 
e-w resol: 30 
n-s resol: 30 
top: 1.000000000000000 
bottom: 0.000000000000000 
cols3: 8001 
rows3: 8111 
depths: 1 
e-w resol3: 30 
n-s resol3: 30 
t-b resol: 1 

DEFAULT_WIND (file) 

proj: 1 
zone: 30 
north: 1 
south: 0 
east: 1 
west: 0 
cols: 1 
rows: 1 
e-w resol: 1 
n-s resol: 1 
top: 1.000000000000000 
bottom: 0.000000000000000 
cols3: 1 
rows3: 1 
depths: 1 
e-w resol3: 1 
n-s resol3: 1 
t-b resol: 1 

PROJ_INFO (file)

name: Universal Transverse Mercator proj: utm zone: 30 no_defs: defined datum: wgs84 ellps: wgs84 towgs84: 0.000,0.000,0.000 

PROJ_UNITS  (file)
unit: Meter units: Meters meters: 1

19 comments:

  1. Thanks for sharing your post. You have nicely presented the concept on this post. I think, you know about Raster To Vector Service.

    ReplyDelete
  2. Effective post you created here . I read your blog and appreciate it . Thanks for shared . Remvoe White Background | Clipping Path | Product Photo Editing

    ReplyDelete
  3. Your post is very much informative and the method shows a great talent . Thanks for sharing.


    Image Manipulation Service

    Retouching Image

    ReplyDelete
  4. Always creative writing from your side. That's why I always recommend your blogs to my friends. Thanks again.
    clipping path service

    ReplyDelete
  5. Much informative article you have written.We expect more in future.

    ReplyDelete
  6. Really enjoyed your article to read.so much helpful.Thanks.keep updating more and more.

    ReplyDelete
  7. It's a great informative post for me.Thousands of thanks for this post.Love the blog also.

    Raster to vector convertion

    ReplyDelete
  8. Pick monkey is a very good app for photo editing.Thanks for sharing
    Thanks for this awesome post. It's really a great app for Bloggers! For the Photo Restoration Services you must have clipping path service to get perfect result.

    ReplyDelete
  9. It's time to work with new project and make more clients. Make photography business more effective by taking our editing service. In fact, if you are really serious about your financial success, you can deal with Clipping Path Unit. We will give you the best editing service within your budget and it's on time.

    ReplyDelete
  10. This comment has been removed by the author.

    ReplyDelete
  11. Beautiful post. I learn a lot of things from you. Thanks.
    fusionmug.com

    ReplyDelete
  12. We offer qualitative image editing services, image masking, clipping path, image restoration, image manipulation, real estate
    image editing services & Raster to vector service .For more information please contact us. clipping path service provider

    ReplyDelete
  13. Very effective article about Raster to vector service .I appreciate that your post is so much informative.Thanks for sharing .
    Clipping Path service .

    ReplyDelete
  14. Stunning post,Just simple and simplicity is the best.like it

    ReplyDelete
  15. The article was up to the point and described the information very effectively. Thanks to blog author for wonderful and informative post.
    infinite logo design
    logo design app

    ReplyDelete
  16. I am very glad to comment on this site with the excellent content of your site
    clipping path service
    Raster To vector

    ReplyDelete

Featured post

Qgis-server...Installing the QGIS Lizmap Plugin & Lizmap Web Client

This post follows on from my previous three (most recent first in list) linked below. There is no doubt that in just a short time from no...