Or you could do it without GRASS, using a script like this.
Roger
--
On Tue, Apr 21, 2009 at 7:41 AM, Dylan Beaudette
<[email protected]>wrote:
> Hi,
>
> Here is a start:
> http://grass.osgeo.org/wiki/Main_Page
>
> the main idea:
>
> # make a new mapset/location
> # import dem:
> r.in.gdal
>
> # convert to classes
> r.reclass
>
> # fix colors:
> r.colors
>
> # export image:
> r.out.tiff or d.out.file
>
> Cheers,
> Dylan
>
> On Tue, Apr 21, 2009 at 6:33 AM, leonidas <[email protected]>
> wrote:
> >
> > Any tutorial how to implement that with grass gis?
> >
> > Dylan Beaudette-2 wrote:
> >>
> >> You could always apply the elevation classification in a standard GIS
> >> like GRASS, and then save the color-mapped image and use that instead
> >> of the original DEM.
> >>
> >> Cheers,
> >>
> >> Dylan
> >>
> >>
> >
> > --
> > View this message in context:
> http://n2.nabble.com/Dem-classification-tp2662869p2670014.html
> > Sent from the Mapserver - User mailing list archive at Nabble.com.
> >
> > _______________________________________________
> > mapserver-users mailing list
> > [email protected]
> > http://lists.osgeo.org/mailman/listinfo/mapserver-users
> >
> _______________________________________________
> mapserver-users mailing list
> [email protected]
> http://lists.osgeo.org/mailman/listinfo/mapserver-users
>
#! /usr/bin/env python
# Create a colored image based on LUT in MakeColor function
# Author: Roger Andre, October 2008
# Usage: discreet_gray2color.py <infile> <outfile>
from osgeo import gdal
import sys
import numpy
import os.path
gdal.TermProgress = gdal.TermProgress_nocb
src_file = sys.argv[1]
dst_file = sys.argv[2]
out_bands = 3
def MakeColor(z_value):
'''LUT for color ramp. Keys are pixel values, hash values are RGB triplets.'''
color_dict = {
0:[102,0,255],
25:[20,82,255],
50:[0,194,224],
75:[0,255,122],
100:[41,255,0],
125:[204,255,0],
150:[245,194,0],
175:[224,118,0],
200:[168,46,0],
225:[105,0,0],
250:[64,0,0],
255:[0,0,0]}
key_list = color_dict.keys()
key_list.sort()
while len(key_list) > 0:
last_val = key_list[-1]
if z_value >= last_val:
return color_dict[last_val]
else:
key_list.remove(last_val)
# Print some info
print "Creating %s" % (dst_file)
# Open source file
src_ds = gdal.Open( src_file )
src_band = src_ds.GetRasterBand(1)
# create destination file
dst_driver = gdal.GetDriverByName('GTiff')
dst_ds = dst_driver.Create(dst_file, src_ds.RasterXSize, src_ds.RasterYSize, out_bands, gdal.GDT_Byte)
# create output bands
band1 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])
band2 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])
band3 = numpy.zeros([src_ds.RasterYSize, src_ds.RasterXSize])
# set the projection and georeferencing info
dst_ds.SetProjection( src_ds.GetProjection() )
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
# read the source file
gdal.TermProgress( 0.0 )
for iY in range(src_ds.RasterYSize):
src_data = src_band.ReadAsArray(0,iY,src_ds.RasterXSize,1)
col_values = src_data[0] # array of z_values, one per row in source data
for iX in range(src_ds.RasterXSize):
z_value = col_values[iX]
# print z_value # uncomment to see what value breaks color ramp
[R,G,B] = MakeColor(z_value)
band1[iY][iX] = R
band2[iY][iX] = G
band3[iY][iX] = B
gdal.TermProgress( (iY+1.0) / src_ds.RasterYSize )
# write each band out
dst_ds.GetRasterBand(1).WriteArray(band1)
dst_ds.GetRasterBand(2).WriteArray(band2)
dst_ds.GetRasterBand(3).WriteArray(band3)
dst_ds = None
_______________________________________________
mapserver-users mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/mapserver-users
