Данные SRTM были построены с помощью следующего скрипта Python 2 (требуется Python Imaging Library и NumPy ), а затем масштабированы в GIMP для корректировки соотношения сторон исходных данных на этой широте (около 92 м x 70 м). Затенение ландшафта и гипсометрические цвета были объединены в GIMP в режиме умножения слоев.
# Прочитать файл SRTM3 и создать затемненный рельеф# 2010-04-05из структуры импорта распаковать , вычислитьfrom numpy import * import numpy as np from PIL import Imagerow_length = 1201 # row_length составляет 1201 для SRTM3 или 3601 для SRTM1 file_name = "N40E014.hgt" # из http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/ hlim = 800 # предел высоты для карты [м]ref_lat = 40.55 # эталонная широта earth_eq = 6371. * 1000. * 2. * pi x_scale = 1. / 360. * earth_eq * cos ( ref_lat / 180. * pi ) / row_length y_scale = 1. / 360. * earth_eq / row_lengthprint "1 пиксель = % u * % u m" % ( x_scale , y_scale ) print "factor" , y_scale / x_scaleh = нули (( row_length , row_length )) f = open ( file_name , 'r' ) li = []для j в диапазоне ( row_length ): для i в диапазоне ( row_length ): d = f . read ( 2 ) ( height ,) = unpack ( '> h' , d ) h [ i , j ] = height, если height < - 1000 : li . добавить (( я , j ))hmax = h . max () h3 = zeros_like ( h ) h3 [:,:] = h [:,:] print len ( li ), «недостающие точки данных»def get_nei ( z ): h2 = h [ z [ 0 ] - 1 : z [ 0 ] + 2 , z [ 1 ] - 1 : z [ 1 ] + 2 ] nn = sum ( где ( h2 < - 1000 , 0 , 1 )) av = sum ( где ( h2 > - 1000 , h2 , 0 )) / float ( nn ) return nn , av# заполняем недостающие точки методом усреднения ближайшего соседа: loop = len ( li ) lim = 7 while loop > 0 : sd = False for q in range ( len ( li )): if h [ li [ q ]] > - 1000. : продолжить п , = get_nei ( Li [ д ]) , если п > = Пт : печать Li [ д ], петли , п , , Пт h3 [ Li [ д ]] = цикл - = 1 с.о. = Истинный если не sd : lim - = 1 h [:,:] = h3 [:,:] выведите "недостающие точки выполнены" def hext ( a ): «Шестнадцатеричный цвет в триплет». r , g , b = a [ 0 : 2 ], a [ 2 : 4 ], a [ 4 : 6 ] return int ( r , 16 ), int ( g , 16 ), int ( b , 16 )# из http://en.wikipedia.org/wiki/Wikipedia:WikiProject_Maps/Conventions/Topographic_maps: col_sea = hext ( "0978ab" ) cols = "" " {{Mapcolor | r = 245 | v = 244 | b = 242 | hex = # F5F4F2 | col = black}} {{Mapcolor | r = 224 | v = 222 | b = 216 | hex = # E0DED8 | col = black}} {{Mapcolor | r = 202 | v = 195 | b = 184 | hex = # CAC3B8 | col = black}} {{Mapcolor | r = 186 | v = 174 | b = 154 | hex = # BAAE9A | col = black}} {{Mapcolor | r = 172 | v = 154 | b = 124 | hex = # AC9A7C | col = black}} {{Mapcolor | r = 170 | v = 135 | b = 83 | hex = # AA8753 | col = black}} {{Mapcolor | r = 185 | v = 152 | b = 90 | hex = # B9985A | col = black}} {{Mapcolor | r = 195 | v = 167 | b = 107 | hex = # C3A76B | col = black}} {{Mapcolor | r = 202 | v = 185 | b = 130 | hex = # CAB982 | col = black}} {{Mapcolor | r = 211 | v = 202 | b = 157 | hex = # D3CA9D | col = black}} {{Mapcolor | r = 222 | v = 214 | b = 163 | hex = # DED6A3 | col = black}} {{Mapcolor | r = 232 | v = 225 | b = 182 | hex = # E8E1B6 | col = black}} {{Mapcolor | r = 239 | v = 235 | b = 192 | hex = # EFEBC0 | col = black}} {{Mapcolor | r = 225 | v = 228 | b = 181 | hex = # E1E4B5 | col = black}} { {Mapcolor | r = 209 | v = 215 | b = 171 | hex = # D1D7AB | col = black}} {{Mapcolor | r = 189 | v = 204 | b = 150 | hex = # BDCC96 | col = black} } {{Mapcolor | r = 168 | v = 198 | b = 143 | hex = # A8C68F | col = black}} {{Mapcolor | r = 148 | v = 191 | b = 139 | hex = # 94BF8B | col = black}} {{Mapcolor | r = 172 | v = 208 | b = 165 | hex = # ACD0A5 | col = black}} "" " col = []для l в кол . splitlines (): если len ( l ) < 10 : продолжить i = l . find ( '#' ), если i > - 1 : col . append ( hext ( l [ i + 1 : i + 7 ]))col . reverse () # -> снизу вверхo = Изображение . новый ( 'RGB' , форма h . )def interp ( c , f ): «Интерполировать в таблицу цветов». r = int (( 1. - f ) * col [ c ] [ 0 ] + f * col [ c + 1 ] [ 0 ]) g = int (( 1. - f ) * col [ c ] [ 1 ] + f * col [ c + 1 ] [ 1 ]) b = int (( 1. - f ) * col [ c ] [ 2 ] + f * col [ c + 1 ] [ 2 ]) return r , g , bfor j in range ( row_length ): for i in range ( row_length ): c , f = divmod ( h [ j , i ] / hmax * ( len ( col ) - 1 ), 1 ) if 0 < h [ j , i ] < hmax : o . putpixel (( j , i ), interp ( int ( c ), f )) elif h [ i , j ] == hmax : o . putpixel (( j , i ), col [ - 1 ]) else : o . putpixel (( j , i ), col_sea )о . save ( "map_height.png" ) # сохраняем карту высот o2 = o . урожай (( 0 , 0 , 942 , 603 )) o2 . сохранить ( "map_height_cropped.png" )# взято из hillshade.py: #def illumination (idata, azdeg = 315.0, altdeg = 45.): def illumination ( idata , azdeg = 225.0 , altdeg = 45. ): # преобразовать alt, az в радианы az = azdeg * np . pi / 180.0 alt = altdeg * np . pi / 180.0 # градиент в направлениях x и y dx , dy = np . градиент ( idata ) slope = 0,5 * np . пи - нп . arctan ( np . hypot ( dx , dy )) аспект = np . arctan2 ( dx , dy ) odata = np . грех ( альт ) * нп . sin ( наклон ) + np . cos ( alt ) * np . cos ( наклон ) * np . cos ( - az - \ aspect - 0.5 * np . pi ) # масштабировать до интервала -1,1 # 1 означает максимальное пребывание на солнце, а 0 означает полную тень. odata = ( odata - odata . min ()) / ( odata . max () - odata . min ()) вернуть odatail = 255 * освещенность ( ч )o4 = Изображение . new ( 'RGBA' , il . shape ) для j в диапазоне ( row_length - 1 ): для i в диапазоне ( row_length - 1 ): v = int ( il [ j , i ]), если 0 <= v < 128 : alpha = ( 255 - 2 * v ) o4 . putpixel (( j , i ), ( 0 , 0 , 0 , alpha )) elif v == 128 : o4 . putpixel (( j , i ), ( 0 , 0 , 0 , 0 )) elif 128 < v < 256 : альфа = 2 * ( v - 128 ) o4 . putpixel (( j , i ), ( 255 , 255 , 255 , alpha )) иначе : o4 . putpixel (( j , i ), ( 255 , 255 , 255 , 0 )) o4 . save ( "il_NW_alpha.png" ) # NW- lighting (альфа-прозрачность для использования с Inkscape)