[OSM-talk-fr] Openhikingmap, OSM et la rando
Stéphane Brunner
stephane.brunner at gmail.com
Lun 19 Jan 09:11:50 UTC 2009
Hello !
Quand je parlais de tuile je pensais au tuile de 1°x1° fournie par la
NASA. et week-end j'en n'ai a nouveau générée quelque unes pour test
avec un void filling fais maison (images réduties a 12%, car
imageshack limite la taille des images, ce qui atenue les défauts de
mon void filling !):
Himalaya :
http://img296.imageshack.us/my.php?image=n28e086colorshadeim6.png
http://img294.imageshack.us/my.php?image=n27e086colorshadesv9.png
Marseille :
http://img297.imageshack.us/my.php?image=n43e005colorshadero1.png
Angleterre :
http://img177.imageshack.us/my.php?image=n52w003colorshadecw7.png
http://img181.imageshack.us/my.php?image=n52w004colorshadega0.png
Suisse :
http://img297.imageshack.us/my.php?image=n43e005colorshadero1.png
http://img299.imageshack.us/my.php?image=n46e006colorshadedd1.png
http://img403.imageshack.us/my.php?image=n46e007colorshadevv2.png
http://img243.imageshack.us/my.php?image=n47e006colorshadezf8.png
http://img172.imageshack.us/my.php?image=n47e007colorshadess3.png
http://img174.imageshack.us/my.php?image=n47e008colorshadeti2.png
Je suis en train d'essayer le script grass : r.fillnulls, ce sera
surement pour la prochaine version.
Es-ce qu'il y a un endroit ou je pourrais envoyer sans problème de
taille les "dalles" que je génère ce qui permettrait de tester plus
loin !
CU
Stéphane
Source code :
Scale.txt
8000 148 103 50
2000 247 171 96
500 242 243 159
100 134 201 116
-500 255 255 255
-32767 0 0 0
color-shade.cpp :
/****************************************************************************
* color-hillshade.cpp
* Author: Stéphane Brunner, Matthew Perry, Paul Surgeon
* License :
Copyright 2005 Matthew T. Perry, 2005 Paul Surgeon, 2008 Stéphane Brunner
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
* calculates a shaded relief image from a gdal-supported raster DEM
CHANGELOG
* updated nodata handling so that 0 is nodata and 1 - 255 are the shade values
****************************************************************************/
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <fstream>
#include <list>
#include "gdal_priv.h"
#include "stringtok.h"
using namespace std;
struct SColor
{
int Red;
int Green;
int Blue;
};
struct SColorPoint
{
int Elevation;
SColor Color;
};
vector<SColorPoint*> ColorPointList;
//=============================================================================
void ReadColorScale(const string& ScaleFileName)
{
ifstream ScaleFile;
string Buffer;
list<string> StringList;
SColorPoint* TempColorPoint;
ScaleFile.open(ScaleFileName.c_str(), ios::in);
if (!ScaleFile.is_open())
{
cout << "Error opening color scale file : " << ScaleFileName << endl;
exit (1);
}
while (!ScaleFile.eof())
{
StringList.clear();
getline(ScaleFile, Buffer);
// Strip spaces in case we have a blank line
while (Buffer[0] == ' ')
{
Buffer.erase(0);
}
// If not a blank line
if (Buffer != "")
{
TempColorPoint = new SColorPoint;
stringtok(StringList, Buffer, " ");
list<string>::iterator i = StringList.begin();
TempColorPoint->Elevation = atoi(string(*i).c_str());
i++;
TempColorPoint->Color.Red = atoi(string(*i).c_str());
i++;
TempColorPoint->Color.Green = atoi(string(*i).c_str());
i++;
TempColorPoint->Color.Blue = atoi(string(*i).c_str());
ColorPointList.push_back(TempColorPoint);
}
}
ScaleFile.close();
}
//=============================================================================
// Given an elevation calculate a color based on the color points table
// At the moment we're only doing linear color gradients
//=============================================================================
SColor GetColor(float Elevation)
{
SColor Color = {0,0,0};
SColorPoint* LowerColorPoint = NULL;
SColorPoint* UpperColorPoint = NULL;
SColorPoint* TempColorPoint;
float TempElev;
float DiffFactor;
// Find closest pair of color points that the elevation falls between
// Lower color point
TempElev = -64000;
for (unsigned int i = 0; i < ColorPointList.size(); i++)
{
TempColorPoint = ColorPointList[i];
if ((TempColorPoint->Elevation <= Elevation) && (TempElev <
TempColorPoint->Elevation))
{
TempElev = TempColorPoint->Elevation;
LowerColorPoint = TempColorPoint;
}
}
// Upper color point
TempElev = 64000;
for (unsigned int i = 0; i < ColorPointList.size(); i++)
{
TempColorPoint = ColorPointList[i];
if ((TempColorPoint->Elevation >= Elevation) && (TempElev >
TempColorPoint->Elevation))
{
TempElev = TempColorPoint->Elevation;
UpperColorPoint = TempColorPoint;
}
}
// I should change the following clipping logic to use the color scale instead
// Return blue
if (LowerColorPoint == NULL)
{
Color.Red = 150;
Color.Green = 150;
Color.Blue = 255;
return Color;
}
// Return white
if (UpperColorPoint == NULL)
{
Color.Red = 255;
Color.Green = 255;
Color.Blue = 255;
return Color;
}
// Work out the factor the elevation is between the lower and upper
color point elevations
// If the upper and lower color points point to the same color point
object then
// it means that the elevation falls exactly on a color point
if (LowerColorPoint != UpperColorPoint)
{
DiffFactor = (Elevation - LowerColorPoint->Elevation) /
(UpperColorPoint->Elevation - LowerColorPoint->Elevation);
Color.Red = (int)((UpperColorPoint->Color.Red -
LowerColorPoint->Color.Red) * DiffFactor) +
LowerColorPoint->Color.Red;
Color.Green = (int)((UpperColorPoint->Color.Green -
LowerColorPoint->Color.Green) * DiffFactor) +
LowerColorPoint->Color.Green;
Color.Blue = (int)((UpperColorPoint->Color.Blue -
LowerColorPoint->Color.Blue) * DiffFactor) +
LowerColorPoint->Color.Blue;
}
else
{
Color.Red = LowerColorPoint->Color.Red;
Color.Green = LowerColorPoint->Color.Green;
Color.Blue = LowerColorPoint->Color.Blue;
}
return Color;
}
void get(float* InPixels, int j, int i, int x, int y, int nYSize, float* out) {
int cx;
int cy;
for (cx = 0 ; cx < x ; cx++) {
for (cy = 0 ; cy < y ; cy++) {
out[cx * y + cy] = InPixels[(j+cx) * nYSize + i+cy];
}
}
}
int main(int argc, char* argv[])
{
GDALDataset* poDataset;
double adfGeoTransform[6];
float* RowRed;
float* RowGreen;
float* RowBlue;
float InPixel;
int i;
int j;
const char* Format = "GTiff";
SColor TempColor;
float* InPixels;
if (argc < 3)
{
cout << "color-shade generates a color relief map from any
GDAL-supported elevation raster and a shaded relief map from any
GDAL-supported elevation raster.";
cout << "The final result is the multiplication of the booth maps." << endl;
cout << endl << "Usage:" << endl;
cout << "color-shade <input_dem> <input_color_scale>
<output_relief_map>" << endl;
cout << "[-a addFactor (0..1, default 0)]" << endl;
cout << " [-z ZFactor (default=1)] [-s scale*
(default=1)]" << endl;
cout << " [-az Azimuth (default=315)] [-alt
Altitude (default=45)]" << endl << endl;
cout << "Notes about color relief map (input_color_scale):" << endl;
cout << "The input color scale is a file containing a set of
elevation points (in meters)" << endl;
cout << "and colors. Typically only a small number of elevation
and color sets will be needed" << endl;
cout << "and the rest will be interpolated by color-relief." << endl;
cout << "Example color scale file with 4000 meters set to white
and 0 meters set to green:" << endl;
cout << "4000 255 255 255" << endl;
cout << "0 0 255 0" << endl << endl;
cout << "Using true black (0 0 0) as your RGB values will yield
blank/null cells." << endl;
cout << "Note that to remove nodata from the output, set the DEM's
nodata value to rgb of 0 0 0:" << endl;
cout << "-32767 0 0 0" << endl << endl;
cout << "See the accompanying \"scale.txt\" file for a decent
example." << endl << endl;
cout << "Notes about shaded relief map:" << endl;
cout << " Scale for Feet:Latlong use scale=370400, for
Meters:LatLong use scale=111120" << endl << endl;;
cout << "Notes map meging:" << endl;
cout << " The add factor can be use to add some britness." << endl;
cout << " the used calcul: result = colorPixel * (addFactor +
shadePixel / 256.0)." << endl << endl;
exit(1);
}
const char* InFilename = argv[1];
const string ScaleFilename = argv[2];
const char* OutFilename = argv[3];
// Open and read color scale file
ReadColorScale(ScaleFilename);
GDALAllRegister();
// Open dataset and get raster band
poDataset = (GDALDataset*) GDALOpen(InFilename, GA_ReadOnly);
if(poDataset == NULL)
{
cout << "Couldn't open dataset " << InFilename << endl;
}
GDALRasterBand *poInBand;
poInBand = poDataset->GetRasterBand(1);
poDataset->GetGeoTransform(adfGeoTransform);
// Get variables from input dataset
const int nXSize = poInBand->GetXSize();
const int nYSize = poInBand->GetYSize();
RowRed = (float *) CPLMalloc(sizeof(float)*nXSize);
RowGreen = (float *) CPLMalloc(sizeof(float)*nXSize);
RowBlue = (float *) CPLMalloc(sizeof(float)*nXSize);
InPixels = (float *) CPLMalloc(sizeof(float)*nXSize*nYSize);
// Create the output dataset and copy over relevant metadata
GDALDriver *poDriver;
poDriver = GetGDALDriverManager()->GetDriverByName(Format);
GDALDataset *poDS;
GDALRasterBand *poBandRed;
GDALRasterBand *poBandGreen;
GDALRasterBand *poBandBlue;
char** Options = NULL;
poDS = poDriver->Create(OutFilename,nXSize,nYSize,3,GDT_Byte,Options);
poDS->SetGeoTransform(adfGeoTransform);
poDS->SetProjection(poDataset->GetProjectionRef());
poBandRed = poDS->GetRasterBand(1);
poBandRed->SetNoDataValue(0);
poBandGreen = poDS->GetRasterBand(2);
poBandGreen->SetNoDataValue(0);
poBandBlue = poDS->GetRasterBand(3);
poBandBlue->SetNoDataValue(0);
const float radiansToDegrees = 180.0 / 3.14159;
const float degreesToRadians = 3.14159 / 180.0;
float *win = (float *) CPLMalloc(sizeof(float)*9);
int winsize = 9;
float *shadeBuf;
float x;
float y;
float aspect;
float slope;
float cang;
int n;
int containsNull;
const char *pszFormat = "GTiff";
float z = 1.0;
float scale = 1.0;
float az = 315.0;
float alt = 45.0;
float add = 0;
for ( int iArg = 4; iArg < argc; iArg++ )
{
if( EQUAL(argv[iArg],"-a") ||
EQUAL(argv[iArg],"-add"))
add = atof(argv[iArg+1]);
if( EQUAL(argv[iArg],"-z") )
z = atof(argv[iArg+1]);
if( EQUAL(argv[iArg],"-s") ||
EQUAL(argv[iArg],"-scale"))
scale = atof(argv[iArg+1]);
if( EQUAL(argv[iArg],"-az") ||
EQUAL(argv[iArg],"-azimuth"))
az = atof(argv[iArg+1]);
if( EQUAL(argv[iArg],"-alt") ||
EQUAL(argv[iArg],"-altitude"))
alt = atof(argv[iArg+1]);
}
GDALAllRegister();
/* -------------------------------------
* Get variables from input dataset
*/
const double nsres = adfGeoTransform[5];
const double ewres = adfGeoTransform[1];
const float inputNullValue = (float) poInBand->GetNoDataValue( );
const float nullValue = 0.0;
shadeBuf = (float *) CPLMalloc(sizeof(float)*nXSize);
/* -----------------------------------------
* Create the output dataset and copy over relevant metadata
*/
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
char **papszOptions = NULL;
for (i = 0; i < nYSize; i++) {
cout << "\r" << i << "/" << nYSize << " ";
for (j = 0; j < nXSize; j++) {
poInBand->RasterIO(GF_Read, j, i, 1, 1, &InPixel, 1, 1,
GDT_Float32, 0, 0);
// reduce missing data impact
if (InPixel == inputNullValue) {
/* int is = i == 0 ? 0 : i - 1;
int ie = i == nYSize-1 ? i : i + 1;
int js = j == 0 ? 0 : j - 1;;
int je = j == nXSize-1 ? j : j + 1;*/
int is = i;
int ie = i;
int js = j;
int je = j;
while (is != 0 && InPixel == inputNullValue) {
is--;
poInBand->RasterIO(GF_Read, j, is, 1, 1, &InPixel, 1, 1, GDT_Float32, 0, 0);
}
InPixel = inputNullValue;
while (ie != nYSize-1 && InPixel == inputNullValue) {
ie++;
poInBand->RasterIO(GF_Read, j, ie, 1, 1, &InPixel, 1, 1, GDT_Float32, 0, 0);
}
InPixel = inputNullValue;
while (js != 0 && InPixel == inputNullValue) {
js--;
poInBand->RasterIO(GF_Read, js, i, 1, 1, &InPixel, 1, 1, GDT_Float32, 0, 0);
}
InPixel = inputNullValue;
while (je != nXSize-1 && InPixel == inputNullValue) {
je++;
poInBand->RasterIO(GF_Read, je, i, 1, 1, &InPixel, 1, 1, GDT_Float32, 0, 0);
}
for (int k = 0 ; k < 50 ; k++) {
is = is == 0 ? 0 : is - 1;
ie = ie == nYSize-1 ? ie : ie + 1;
js = js == 0 ? 0 : js - 1;
je = je == nXSize-1 ? je : je + 1;
}
int length = (ie -is + 1) * (je - js + 1);
//if (length>100000) cout << endl << length << ", " << i-is << ", " <<
ie-i << ", " << j-js << ", " << je-j << endl;
if (winsize < length) {
CPLFree(win);
win = (float *) CPLMalloc(sizeof(float)*length);
winsize = length;
}
int il = ie -is + 1;
int jl = je -js + 1;
poInBand->RasterIO(GF_Read, js, is, jl, il, win, jl, il,
GDT_Float32, 0, 0);
float fact = 0;
float sum = 0;
for (int k = 0 ; k < length ; k++) {
if (win[k] != inputNullValue) {
int ik = is+k/jl - i;
int jk = js+k%jl - j;
float length = ik*ik+jk*jk;
length = length*length;
sum += win[k]/length;
fact += 1/length;
}
}
//if (length>100000) cout << endl << sum << "/" << fact << "=" << sum
/ fact << endl;
InPixels[j * nYSize + i] = sum / fact;
// Write lines to output raster
/* int is = i;
int ie = i;
int js = j;
int je = j;
int length = 1;
int valide = 0;
float sum = 0;
float fact = 0;
while ((float)length / valide > 2) {
int nb = sqrt(length);
for (int z = 0 ; z < nb ; z++) {
is = is == 0 ? 0 : is - 1;
ie = ie == nYSize-1 ? ie : ie + 1;
js = js == 0 ? 0 : js - 1;
je = je == nXSize-1 ? je : je + 1;
}
length = (ie -is + 1) * (je - js + 1);
if (length > 100000) cout << "\r" << i << "/" << nYSize << " - " << j
<< "/" << nXSize << ", " << length << ". ";
valide = 0;
fact = 0;
sum = 0;
if (winsize < length) {
CPLFree(win);
win = (float *) CPLMalloc(sizeof(float)*length);
winsize = length;
}
int il = ie -is + 1;
int jl = je -js + 1;
poInBand->RasterIO(GF_Read, js, is, jl, il, win, jl, il,
GDT_Float32, 0, 0);
for (int k = 0 ; k < length ; k++) {
if (win[k] != inputNullValue) {
valide++;
int ik = is+k/jl - i;
int jk = js+k%jl - j;
float length = sqrt(ik*ik+jk*jk);
sum += win[k]/length;
fact += 1/length;
}
}
}
//cout << "\r" << sum << "/" << fact << " = " << sum / fact << endl;
InPixels[j * nYSize + i] = sum / fact;*/
}
else {
InPixels[j * nYSize + i] = InPixel;
}
}
}
cout << endl;
/* ------------------------------------------
* Move a 3x3 window over each cell
* (where the cell in question is #4)
*
* 0 1 2
* 3 4 5
* 6 7 8
*
*/
for ( i = 0; i < nYSize; i++) {
for ( j = 0; j < nXSize; j++) {
containsNull = 0;
// Exclude the edges
// if (i == 0 || j == 0 || i == nYSize-1 || j == nXSize-1 )
// {
// We are at the edge so write nullValue and move on
// shadeBuf[j] = nullValue;
// continue;
// }
// For the edges
if (i == 0)
{
if (j == 0)
{
get(InPixels, j, i, 2, 2, nYSize, win);
win[8] = win[3];
win[7] = win[2];
win[5] = win[1];
win[4] = win[0];
win[0] = 2*win[4] - win[8];
win[1] = 2*win[4] - win[7];
win[2] = 2*win[5] - win[8];
win[3] = 2*win[4] - win[5];
win[6] = 2*win[7] - win[8];
}
else if (j == nXSize-1)
{
get(InPixels, j-1, i, 2, 2, nYSize, win);
win[5] = win[3];
win[4] = win[2];
win[2] = win[1];
win[1] = win[0];
win[0] = 2*win[1] - win[2];
win[3] = 2*win[4] - win[5];
win[6] = 2*win[4] - win[8];
win[7] = 2*win[4] - win[1];
win[8] = 2*win[5] - win[2];
}
else
{
get(InPixels, j-1, i, 3, 2, nYSize, win);
win[8] = win[5];
win[7] = win[4];
win[5] = win[3];
win[4] = win[2];
win[2] = win[1];
win[1] = win[0];
win[0] = 2*win[1] - win[2];
win[3] = 2*win[4] - win[5];
win[6] = 2*win[7] - win[8];
}
}
else if (i == nYSize-1)
{
if (j == 0)
{
get(InPixels, j, i-1, 2, 2, nYSize, win);
win[7] = win[3];
win[6] = win[2];
win[4] = win[1];
win[3] = win[0];
win[0] = 2*win[3] - win[6];
win[1] = 2*win[4] - win[7];
win[2] = 2*win[4] - win[6];
win[5] = 2*win[4] - win[3];
win[8] = 2*win[7] - win[6];
}
else if (j == nXSize-1)
{
get(InPixels, j-1, i-1, 2, 2, nYSize, win);
win[4] = win[3];
win[3] = win[2];
win[1] = win[1];
win[0] = win[0];
win[2] = 2*win[1] - win[0];
win[5] = 2*win[4] - win[3];
win[6] = 2*win[3] - win[0];
win[7] = 2*win[4] - win[1];
win[8] = 2*win[4] - win[0];
}
else
{
get(InPixels, j-1, i-1, 3, 2, nYSize, win);
win[7] = win[5];
win[6] = win[4];
win[4] = win[3];
win[3] = win[2];
win[1] = win[1];
win[0] = win[0];
win[2] = 2*win[1] - win[0];
win[5] = 2*win[4] - win[3];
win[8] = 2*win[7] - win[6];
}
}
else
{
if (j == 0)
{
get(InPixels, j, i-1, 2, 3, nYSize, win);
win[8] = win[5];
win[7] = win[4];
win[6] = win[3];
win[5] = win[2];
win[4] = win[1];
win[3] = win[0];
win[0] = 2*win[3] - win[6];
win[1] = 2*win[4] - win[7];
win[2] = 2*win[5] - win[8];
}
else if (j == nXSize-1)
{
get(InPixels, j-1, i-1, 2, 3, nYSize, win);
win[5] = win[5];
win[4] = win[4];
win[3] = win[3];
win[2] = win[2];
win[1] = win[1];
win[0] = win[0];
win[6] = 2*win[3] - win[0];
win[7] = 2*win[4] - win[1];
win[8] = 2*win[5] - win[2];
}
else
{
// Read in 3x3 window
get(InPixels, j-1, i-1, 3, 3, nYSize, win);
}
}
// Check if window has null value
for ( n = 0; n <= 8; n++) {
if(win[n] == inputNullValue) {
containsNull = 1;
break;
}
}
if (containsNull == 1) {
// We have nulls so write nullValue and move on
shadeBuf[j] = nullValue;
continue;
} else {
// We have a valid 3x3 window.
/* ---------------------------------------
* Compute Hillshade
*/
// First Slope ...
x = ((z*win[0] + z*win[3] + z*win[3] + z*win[6]) -
(z*win[2] + z*win[5] + z*win[5] + z*win[8])) /
(8.0 * ewres * scale);
y = ((z*win[6] + z*win[7] + z*win[7] + z*win[8]) -
(z*win[0] + z*win[1] + z*win[1] + z*win[2])) /
(8.0 * nsres * scale);
slope = 90.0 - atan(sqrt(x*x + y*y))*radiansToDegrees;
// ... then aspect...
aspect = atan2(x,y);
// ... then the shade value
cang = sin(alt*degreesToRadians) * sin(slope*degreesToRadians) +
cos(alt*degreesToRadians) * cos(slope*degreesToRadians) *
cos((az-90.0)*degreesToRadians - aspect);
if (cang <= 0.0)
cang = 1.0;
else
cang = 1.0 + (254.0 * cang);
shadeBuf[j] = cang;
}
InPixel = InPixels[j * nYSize + i];
// reduce missing data impact
/* if (shadeBuf[j] == nullValue && j != 0) {
RowRed[j] = RowRed[j-1];
RowGreen[j] = RowGreen[j-1];
RowBlue[j] = RowBlue[j-1];
} else {*/
TempColor = GetColor(InPixel);
float facteur = add + shadeBuf[j] / 256.0;
RowRed[j] = TempColor.Red * facteur;
RowGreen[j] = TempColor.Green * facteur;
RowBlue[j] = TempColor.Blue * facteur;
// }
}
// Write lines to output raster
poBandRed->RasterIO(GF_Write, 0, i, nXSize, 1, RowRed, nXSize, 1,
GDT_Float32, 0, 0);
poBandGreen->RasterIO(GF_Write, 0, i, nXSize, 1, RowGreen, nXSize,
1, GDT_Float32, 0, 0);
poBandBlue->RasterIO(GF_Write, 0, i, nXSize, 1, RowBlue, nXSize,
1, GDT_Float32, 0, 0);
}
delete poDS;
return 0;
}
Plus d'informations sur la liste de diffusion Talk-fr