1

I have a number of rasters (rgb images) that I want to re-project and crop. I have established a workflow for doing this in QGIS, but I'd like to do it using gdal modules, because the output is being imported into a GRASS location anyway and I could more easily automate the workflow if I were using gdal for the whole process.

Here's my QGIS workflow: I have a QGIS project that with the CRS (EPSG:32142) set to the target CRS and the desired extent bookmarked (xmin=-985631.39, ymin=7464178.33, xmax=2296529.53, ymax=9628782.84).

I import the raster that I want to convert, zoom to my extent, choose save layer as . . . In the 'save as' dialog box, I choose:

output mode=raw data extent=map view extent resolution=layer resolution

I understand how to use gdalwarp to reproject and set the extent of the output like this:

gdalwarp -t_srs EPSG:32142 \
  inputfile.tif \
  -te xmin ymin xmax ymax \
  targetfile.tif

I also understand (thanks to MappaGnosis) that I can set the resolution for the target file using the -tr argument to gdalwarp. I don't know how to get the source file's native resolution programatically as MappaGnosis suggested in my previous post, and what I'm imagining is a messy process, like using awk or Python to parse the output of gdalinfo to create the required -tr argument for gdalwarp.

  • If you don't specify the target resolution, gdalwarp does the following: "A resolution is computed with the intent that the length of the distance from the top left corner of the output imagery to the bottom right corner would represent the same number of pixels as in the source image." Wouldn't this built-in algorithm be good enough? – Hermann Oct 2 '13 at 11:52
  • That is not the behavior I observed when I did not specify a target resolution. Rather, the target resolution was significantly lower than the source. – Gregory Oct 2 '13 at 16:51
  • Can you provide an example? – Hermann Oct 3 '13 at 5:34
0

Reprojecting and Cropping

The approach to reproduce this in GRASS is as follows (instructions concern GRASS GIS' command line):

  1. Create & enter a Location based on EPSG:32142 (grass64 -c "epsg:32142" /GRASS_DB_Directory/Location_Name or grass7 -c "epsg:32647" ... for GRASS GIS ver. 7)
  2. Import the imagery to be processed (r.in.gdal)
  3. Define the desired computational region extent & resolution (g.region)
  4. "Copy" data (r.mapcalc). The result will be limited inside the defined computational region extent.

If the "re-projection" concerns re-projecting the source data from their "original" reference system which is another than the (only one) mentioned in the original question (EPSG:32142), you need to add 1-2 extra steps. Specifically,

  1. Create & enter a Location based on your data's CRS (given you have GeoTIFFs beforehand, the instruction could be as follows: grass64 -c GeoTIFF_Raster.tif /GRASS_DB_Directory/Location_Name)
  2. Import the imagery to be processed (r.in.gdal)
  3. Exit from the current location (or open another shell in parallel) and create & enter into a Location based on the target CRS, i.e. if the target CRS is EPSG:32142, do as in Step 1 in the first set of example instructions
  4. Define the desired computational region extent (g.region)
  5. Reproject the data in to the target CRS by restricting the operation within the defined region extent (r.proj along with the -r flag)

MASK-ing

The step of defining the desired computational region extent (ATTENTION, and Resolution!) can be replaced by using another Raster or Vector Map as a MASK (r.mask).

Using the location option with r.in.gdal

From the r.in.gdal and v.in.ogr manuals:

it is possible to have v.in.ogr/ r.in.gdal automatically create a new location based on the projection and extents of the file being read. This is accomplished by passing the name to be used for the new location via the location parameter. Upon completion of the command, a new location will have been created (with only a PERMANENT mapset), and the raster will have been imported with the indicated output name into the PERMANENT mapset.

Therefore, as noted in the GRASS-Wiki page Map Reprojection:

Another (similar) strategy frequently needed when maps need to be imported using, e.g. v.in.ogr or r.in.gdal, is to tell these programs to create a new location (argument "location"). If the imported data contains proper PROJ_INFO files, they can be reprojected into the current mapset, using the the projection in the current location/mapset.

Finally, no matter which of the above mini-workflows, export your data by utilising the r.out.gdal module. There are certainly even more "smart" details to share if more details are provided in the original question (how many images, potential RGB composites, et.c.)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.