Hello ! Ce weekend je me suis posé des questions pour avoir le meilleurs rendu possible des données SRTM. Je trouve que l'image disponible sur la page http://wiki.openstreetmap.org/wiki/HikingBikingMaps donne un bon rendu de ces donnée mais n'est viable que pour de petites échelles ! Les que l'on s'éloigne trop les courbe de niveau ne donne pas quelque chose de satisfaisant ! Il faux donc ce rapprocher de rendu comme celui-ci : http://www.openstreetmap.org/?lat=43.469&lon=5.096&zoom=9&layers=00B0FTF Qui es déjà très bien mais est loin d'être parfais surtout en altitudes : - Les Alpes : http://www.openstreetmap.org/?lat=45.846&lon=7.203&zoom=9&layers=00B0FTF - L'Himalaya (on ne voie plus qu'un énorme plafond blanc) : http://www.openstreetmap.org/?lat=-71.5&lon=31.8&zoom=3&layers=00B0FTF
A cela je voie 2 problèmes : - La dernière couleurs est affectée a une altitude de 3500 ce qui es trop bas pour des très hautes altitude. - Je pesse que mettre une couleur blanche en altitude ne sois pas une bonne idée, cela fais panser à la neige mais ça à l'inconvénient de ne pas s'adapter a toutes les chaîne de montagnes. Je pensse que c'est mieus d'utiliser une coulleur rouge comme sur l'image su site de MapNik (http://www.mapnik.org/). J'ai fais quelque essais pour arriver à en résultat qui me satisfais sans toutefois être parfais : http://img158.imageshack.us/my.php?image=testyf7.png Pour cela j'ai utiliser une fichier scale.txt comme celui-ci : 8000 246 153 96 2000 247 199 152 600 242 243 159 200 134 201 116 -500 255 255 255 -10000 17 110 176 -32767 0 0 0 et j'ai en quelque sorte mergé les 2 programme hillshade et color-relief pour avoir une image avec des couleurs pour le relief et des ombres (voire attachement en espérant que cela passe !). que j'ai exécuter avec la commande : color-shade "${X%%.hgt.zip}"-warped.tif scale.txt "${X%%.hgt.zip}"-color-shade.tif -z 2 -s 1.8 -a 0.2 Je n'ai pas réussi à résoudre correctement le problème des données manquantes (j'ai fais en hack pour diminuer le problème) et pourtant il est corrigé sur OpenCycleMap. Es-ce que quelqu'un sais comment ? Que-ce que vois en penser ? Es-ce que vous avez des remarques ? Que faire pour aller plus loin ? CU Stéphane 2009/1/9 sylvain letuffe <sylv...@letuffe.org>: > >> Oui mais justement sans gestion de canal alpha, on risque d'avoir des >> très moches effets de crènelages, nop ? > > Y'a pas que les problèmes de crénelage si IE le supporte mal, avec mes essais, > si on ne place pas de canal alpha (param fill-opacity) sur les forêts, les > prairies (et même certaines icones sans canal alpha qui sont répétées en > motif), bref, toutes les surfaces, ça se pose par dessus les courbes de > niveau et c'est sacrément moche. > > j'espère avoir un peu de temps pour rapidement faire l'essai et ré-écrouler ma > bécanne ;-) > > Je vais pouvoir mettre 1Go de cache à la poubelle qui avait été > douloureusement accouchés au fûr et à mesure des navigations > > -- > sly > > _______________________________________________ > Talk-fr mailing list > Talk-fr@openstreetmap.org > http://lists.openstreetmap.org/listinfo/talk-fr > -- Stéphane Brunner mail : stephane.brun...@gmail.com messageries instantanées : stephane.brun...@gmail.com (http://talk.google.com) -------------------------------------- http://www.ubuntu-fr.org - Distribution Linux http://fr.wikipedia.org - Encyclopédie communautaire http://mozilla-europe.org - Navigateur internet / Client de messagerie http://framasoft.net - Annuaire de logiciel libre (gratuit) http://jeuxlibres.net - Jeux Libres (gratuit) http://openstreetmap.org - Cartographie libre (en développement) -------------------------------------- Il existe 10 sortes de personnes : celles qui connaissent le binaire, et les autres. -------------------------------------- Si un jour on te reproche que ton travail n'est pas un travail de professionnel, dis toi que : Des amateurs ont construit l'arche de Noé, et des professionnels le Titanic.
/**************************************************************************** * 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; } 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; 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 *poBand; poBand = poDataset->GetRasterBand(1); poDataset->GetGeoTransform(adfGeoTransform); // Get variables from input dataset const int nXSize = poBand->GetXSize(); const int nYSize = poBand->GetYSize(); RowRed = (float *) CPLMalloc(sizeof(float)*nXSize); RowGreen = (float *) CPLMalloc(sizeof(float)*nXSize); RowBlue = (float *) CPLMalloc(sizeof(float)*nXSize); // 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 *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) poBand->GetNoDataValue( ); const float nullValue = 0.0; shadeBuf = (float *) CPLMalloc(sizeof(float)*nXSize); win = (float *) CPLMalloc(sizeof(float)*9); /* ----------------------------------------- * Create the output dataset and copy over relevant metadata */ poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); char **papszOptions = NULL; /* ------------------------------------------ * 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) { poBand->RasterIO( GF_Read, j, i, 2, 2, win, 2, 2, GDT_Float32, 0, 0 ); 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) { poBand->RasterIO( GF_Read, j-1, i, 2, 2, win, 2, 2, GDT_Float32, 0, 0 ); 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 { poBand->RasterIO( GF_Read, j-1, i, 3, 2, win, 3, 2, GDT_Float32, 0, 0 ); 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) { poBand->RasterIO( GF_Read, j, i-1, 2, 2, win, 2, 2, GDT_Float32, 0, 0 ); 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) { poBand->RasterIO( GF_Read, j-1, i-1, 2, 2, win, 2, 2, GDT_Float32, 0, 0 ); 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 { poBand->RasterIO( GF_Read, j-1, i-1, 3, 2, win, 3, 2, GDT_Float32, 0, 0 ); 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) { poBand->RasterIO( GF_Read, j, i-1, 2, 3, win, 2, 3, GDT_Float32, 0, 0 ); 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) { poBand->RasterIO( GF_Read, j-1, i-1, 2, 3, win, 2, 3, GDT_Float32, 0, 0 ); 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 poBand->RasterIO( GF_Read, j-1, i-1, 3, 3, win, 3, 3, GDT_Float32, 0, 0 ); } } // 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; } poBand->RasterIO(GF_Read, j, i, 1, 1, &InPixel, 1, 1, GDT_Float32, 0, 0); // 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; }
_______________________________________________ Talk-fr mailing list Talk-fr@openstreetmap.org http://lists.openstreetmap.org/listinfo/talk-fr
