用python编程做一版全球地形渲染图。
1. 下载数据
利用python批量下载全球30m分辨率地形数据。
代码：
#!/usr/bin/python# _*_ coding: utf-8 _*_from urllib.request import urlopenfrom urllib.request import Requestimport urllibimport reimport osimport syshttp_adress= 'https://step.esa.int/auxdata/dem/SRTMGL1/'req = Request(http_adress)f = urlopen(req)localDir = 'N:/Globel_DEM/SRTM_Global/'localDir_s = [f for f in localDir.split('/') if f != '']tmp_path = localDir_s[0]for tmp_dir in localDir_s[1:]: tmp_path = os.path.join(tmp_path,tmp_dir) if not os.path.exists(tmp_path): os.mkdir(tmp_path)lines = f.readlines()urlList = []for eachLine in lines: print(eachLine) eachLine = eachLine.decode('utf8') wordList = eachLine.split('\"') for word in wordList: #print(word) if re.match('.*hgt.zip$', word): urlList.append(word)n_file = 0total_number = len(urlList)for everyFile in urlList: n_file += 1 everyURL = http_adress+everyFile print('Downloading...',n_file,'--of--',total_number,'|', everyFile) localFile = localDir + everyFile if not os.path.isfile(localFile): try: urllib.request.urlretrieve(everyURL, localFile) except: continue
运行成功后，将全球地形数据下载到本地，14296个压缩文件，87G。
2. 批量解压
利用python批量解压下载的压缩文件。
代码：
import globimport osfrom osgeo import gdalimport tracebackimport zipfiledef extract_zip(zip_file, out_dir = None): try: t = zipfile.ZipFile(zip_file,'r') for fz_file in t.namelist(): t.extract(fz_file,out_dir) t.close() return True except Exception as ex: print("have exception!") traceback.print_exc() print('can not extract zip file:',os.path.basename(zip_file)) return Falseinpath = 'N:\\Globel_DEM\SRTM_Global\\'unzippath = 'N:\\Globel_DEM\SRTM_Global_unzip\\'if not os.path.exists(unzippath): os.mkdir(unzippath)files = glob.glob(inpath + '*.zip')n_file = len(files)i = 0for zip_file in files: i += 1 print(i,'--of--',n_file,zip_file) extract_zip(zip_file,out_dir = unzippath)
运行结果为解压后全球地形数据，hgt格式。每一个zip压缩文件对应一个hgt高程数据。hgt格式数据可以用QGIS或者ARCGIS打开。
单个高程文件为单波段数据，16位整型，图像越暗，表示海拔越低。打开效果为：
3. 绘制影像分布图
近1.5万个文件，直接打开的话，会花去很长很长时间，很可能会把软件卡死。
先绘制影像分布图，了解影像分布情况。
代码：
#!/usr/bin/python# *-* coding:utf-8import osimport sysfrom osgeo import gdal,ogr,osrimport globsrc_suffix = '.hgt'inpath = 'N:\\Globel_DEM\\SRTM_Global_unzip\\'outpath = inpath#[0:-2] + 'small_image//'if not os.path.exists(outpath): os.mkdir(outpath)if not os.path.exists(inpath): sys.exit(1)#return Noneoutfile_csv = outpath + '/file_list_20210110.csv'shp_file = outpath + 'img_center.shp'prj = osr.SpatialReference()prj.ImportFromEPSG(4326) # wgs-84 Geographical coordinates EPSGprj.MorphToESRI()prj_file_name = shp_file.replace('.shp','.prj')if not os.path.exists(prj_file_name): prjfile = open(prj_file_name, 'w') prjfile.write(prj.ExportToWkt()) prjfile.close()DriverName = "ESRI Shapefile" # e.g.: GeoJSON, ESRI ShapefileFileName = shp_file#'test.shp'driver = ogr.GetDriverByName(DriverName)if os.path.exists(FileName): driver.DeleteDataSource(FileName)# Create the output shapefileoutDataSource = driver.CreateDataSource(FileName)#outDataSource.SetProjection(prj)#outLayer = outDataSource.CreateLayer("center_pt", geom_type=ogr.wkbPoint) idField = ogr.FieldDefn("id", ogr.OFTInteger)outLayer.CreateField(idField)idField = ogr.FieldDefn("name", ogr.OFTString)outLayer.CreateField(idField)idField = ogr.FieldDefn("fullname", ogr.OFTString)outLayer.CreateField(idField)out_file_s = open(outfile_csv,'w')tif_files = glob.glob(inpath + '*' + src_suffix )total_number = len(tif_files)print(total_number)n_file=0#interval = [3,5,7]for tif_file in tif_files: n_file = n_file + 1 tif_file = tif_file#[:-4] file_base = os.path.basename(tif_file) if n_file % 10 == 0: print('reading...',n_file,'--of--',total_number,'___',os.path.basename(tif_file)) try: ds_raw = gdal.Open(tif_file) except: with open(tif_file + '__________BAD_file_can_not_be_read.txt','w') as f: f.close() continue if ds_raw == None: with open(tif_file + '__________BAD_file_can_not_be_read.txt','w') as f: f.close() continue in_gt = ds_raw.GetGeoTransform() x_origin = in_gt[0] y_origin = in_gt[3] x_gsd = in_gt[1] y_gsd = in_gt[5] x_size = ds_raw.RasterXSize y_size = ds_raw.RasterYSize x_center = x_origin + x_size/2 * x_gsd y_center = y_origin + y_size/2 * y_gsd x_min = x_origin x_max = x_min + x_size * x_gsd y_max = y_origin y_min = y_origin + y_size * y_gsd img_nbands = ds_raw.RasterCount point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(x_center,y_center) featureDefn = outLayer.GetLayerDefn() feature = ogr.Feature(featureDefn) feature.SetGeometry(point) feature.SetField("id", n_file) feature.SetField("name", os.path.basename(tif_file)) feature.SetField("fullname", tif_file) outLayer.CreateFeature(feature) feature = None print("{},{},{},{},{},{},{}".format(tif_file,x_center,y_center,x_min,y_min,x_max,y_max),file=out_file_s)outDataSource = Noneout_file_s.close()
运行完成后，生成了影像列表文件和shp格式影像分布图。
全球
看看中国周边的分布图：
4. 投影转换
原地形高程数据为经纬度投影，平面距离单位为度，高程单位为米，平面距离和高程高度单位不一致，会导致制作的地形山体阴影图和设色地形图错误。因此需要进行投影转换，把所有高程数据都转换为平面距离单位为米的投影方式。目标投影选定为WGS84 WebMercator。
先绘制出图边框。
代码：
import globimport osimport sysfrom osgeo import gdaldef gecoord_to_webmercator(gecoord_x, gecoord_y, ge_level): resolution = 20037508.342787001 * 2 / 256 / (2 ** ge_level) x = gecoord_x * resolution - 20037508.342787001 y = 20037508.342787001 - gecoord_y * resolution #lon, lat = webmercator_to_lonlat(x, y) return x,y#lon, latdef draw_matched_sql_table_range_shp(table_name_list,ge_level,dst_suffix = '.tif',shp_file = None): if shp_file is None: return if not os.path.exists(os.path.dirname(shp_file)): return #if sql_file_list: if len(table_name_list) == 0: return print(shp_file) prj = osr.SpatialReference() prj.ImportFromEPSG(3857) # wgs-84 Geographical coordinates EPSG prj.MorphToESRI() prj_file_name = shp_file.replace('.shp','.prj') if not os.path.exists(prj_file_name): prjfile = open(prj_file_name, 'w') prjfile.write(prj.ExportToWkt()) prjfile.close() DriverName = "ESRI Shapefile" # e.g.: GeoJSON, ESRI Shapefile FileName = shp_file#'test.shp' driver = ogr.GetDriverByName(DriverName) if os.path.exists(FileName): driver.DeleteDataSource(FileName) # Create the output shapefile outDataSource = driver.CreateDataSource(FileName) #outDataSource.SetProjection(prj)# outLayer = outDataSource.CreateLayer("img_range", geom_type=ogr.wkbPolygon) idField = ogr.FieldDefn("id", ogr.OFTInteger) outLayer.CreateField(idField) idField = ogr.FieldDefn("name", ogr.OFTString) outLayer.CreateField(idField) idField = ogr.FieldDefn("fullname", ogr.OFTString) outLayer.CreateField(idField) #table_names = [ table_name_x + '_' + table_name_y for table_name_x,table_name_y in zip(sql_file_list[:,3],sql_file_list[:,4]) ] table_names = np.unique(np.array(table_name_list)) print(len(table_names)) ge_level_str = generate_ge_level_str(ge_level) dst_Res = 20037508.342787001 * 2 / 256 / (2 ** ge_level) psize_x = dst_Res psize_y = -dst_Res n_file = 0 for table_name in table_names: n_file += 1 i,j = table_name.split('_') table_x = int(i) table_y = int(j) ulx,uly = gecoord_to_webmercator((table_x) * 64 * 256,(table_y )* 64 * 256,ge_level) lrx,lry = gecoord_to_webmercator((table_x+1) * 64 * 256,(table_y+1)* 64 * 256,ge_level) #outrange = [ulx,uly,lrx,lry] sql_x = table_x//4 sql_y = table_y//4 if sql_x < 0 or sql_y < 0: continue #print('out_range_raw:',outrange) table_filename = ge_level_str + '_' + str(sql_x) + '_' + str(sql_y) + '_' + str(table_x) + '_' + str(table_y) + dst_suffix ulx = int(ulx / psize_x) * psize_x uly = int(uly / psize_y) * psize_y lrx = int(lrx / psize_x) * psize_x lry = int(lry / psize_y) * psize_y out_range = [ulx,uly,lrx,lry] ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(ulx,uly) ring.AddPoint(lrx,uly) ring.AddPoint(lrx,lry) ring.AddPoint(ulx,lry) ring.AddPoint(ulx,uly) poly = ogr.Geometry(ogr.wkbPolygon) poly.AddGeometry(ring) featureDefn = outLayer.GetLayerDefn() feature = ogr.Feature(featureDefn) feature.SetGeometry(poly) feature.SetField("id", n_file) feature.SetField("name", os.path.basename(table_filename)) feature.SetField("fullname", table_filename) outLayer.CreateFeature(feature) feature = None outDataSource = Nonedef match_ge_scope(src_file_list,ge_level,is_geographic = None): import math sql_file_list = list() fo = open(src_file_list, "r") #str =''#'' all_file_lists = fo.readlines( ) fo.close() i = 0 for file_name in all_file_lists: i = i+1 file_info_s = file_name[:-1].split(',') if len(file_info_s) != 7: continue min_x = float(file_info_s[3]) min_y = float(file_info_s[4]) max_x = float(file_info_s[5]) max_y = float(file_info_s[6]) table_x0 = 0 table_y0 = 0 if not is_geographic: if (np.abs(max_x - min_x) < 10) or (np.abs(max_y - min_y) < 10): is_geographic = True if is_geographic: min_x,min_y = lonlat_to_webmercator(min_x,min_y) max_x,max_y = lonlat_to_webmercator(max_x,max_y) gecoord_x_min,gecoord_y_min = webmercator_to_gecoord(min_x,max_y,ge_level) gecoord_x_max,gecoord_y_max = webmercator_to_gecoord(max_x,min_y,ge_level) table_x_min = int(math.floor(gecoord_x_min/256/64)) table_x_max = int(math.ceil(gecoord_x_max/256/64)) table_y_min = int(math.floor(gecoord_y_min/256/64)) table_y_max = int(math.ceil(gecoord_y_max/256/64)) for table_x in range(table_x_min,table_x_max): for table_y in range(table_y_min,table_y_max): sql_x = int(table_x//16) sql_y = int(table_y//16) if (table_x0 != table_x) or (table_y0 != table_y): sql_file_list.append([file_info_s[0],sql_x,sql_y,table_x,table_y]) table_x0 = table_x table_y0 = table_y return sql_file_list#-----------------src_paths = ['']out_path = ''src_suffix = 'hgt'ge_level = 13dst_suffix = '.tif'is_build_overview = Truedst_Res = 20037508.342787001 * 2 / 256 / (2 ** ge_level)psize_x = dst_Respsize_y = -dst_Reswkt_str = ''' PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere", GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]], PROJECTION["Mercator_Auxiliary_Sphere"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",0.0], PARAMETER["Standard_Parallel_1",0.0], PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] '''sr = osr.SpatialReference(wkt_str)projection = srsrc_file_list = 'N:\\Globel_DEM\\SRTM_Global_unzip\\file_list_20210110.csv'outpath = 'N:\\Globel_DEM\\SRTM_Global_webmercator\\'if not os.path.exists(outpath): os.mkdir(outpath)shp_file = os.path.join(os.path.dirname(src_file_list),'out_range_tk.shp')sql_list = match_ge_scope(src_file_list,13,is_geographic=True)sql_list = np.array(sql_list)print(sql_list.shape)file_names = (sql_list[:,0]).tolist()vrt_file = os.path.join(os.path.dirname(src_file_list),'total.vrt')if not os.path.isfile(vrt_file): #files = glob.glob(out_file_path + "*[0-9].tif") vrt_ds = gdal.BuildVRT(vrt_file, file_names) del vrt_dstable_names = [ table_name_x + '_' + table_name_y for table_name_x,table_name_y in zip(sql_list[:,3],sql_list[:,4]) ]table_names = np.unique(np.array(table_names))print(len(table_names))#ge_level_str = generate_ge_level_str(ge_level)if not os.path.exists(shp_file): draw_matched_sql_table_range_shp(table_names,ge_level,shp_file = shp_file)poly_list = read_shp(shp_file)try: vrt_ds = gdal.Open(vrt_file)except: print('can not open vrt file, exit...') sys.exit(1)total_number = len(poly_list)n = 0for poly in poly_list: #print(poly) n += 1 print(n,'---of---',total_number,poly[0]) out_file_name = os.path.join(outpath,poly[0]) ulx,lry,lrx,uly = poly[1:] buffer_size = 16 * dst_Res out_range = [ulx-buffer_size,lry-buffer_size,lrx+buffer_size,uly+buffer_size] print((lrx-ulx)/dst_Res,(-lry+uly)/dst_Res) #print(np.array(out_range)/dst_Res) print('out_range:',out_range) #break out_file_processing = out_file_name + '_processing.txt' out_file_processed = out_file_name + '_processed.txt' if os.path.exists(out_file_processed): continue if os.path.exists(out_file_processing): continue if os.path.isfile(out_file_name): continue with open(out_file_processing,'w') as f: f.close() if not os.path.isfile(out_file_name): print('warping...',out_file_name) out_ds = gdal.Warp(out_file_name, vrt_ds, options=gdal.WarpOptions( creationOptions=["TILED=YES"], srcSRS='EPSG:4326', dstSRS=projection, xRes=psize_x, yRes=-psize_y, outputBounds = out_range, resampleAlg = gdal.GRIORA_Bilinear, multithread = True, warpOptions = ['NUM_THREADS=ALL_CPUS'], dstNodata = -32768)) data = out_ds.GetRasterBand(1).ReadAsArray() if data.min() == data.max(): out_ds = None #is_build_overview = False try: os.remove(out_file_name) except: print('can not remove blank file...continue') continue else: if is_build_overview: #print('building overview') print('building overview...') out_ds.BuildOverviews('nearest',[2,4,8,16,32,64,128]) out_ds.FlushCache() out_ds = None #else: out_ds = None with open(out_file_processed,'w') as f: f.close() try: if os.path.exists(out_file_processing): os.remove(out_file_processing) except: None
运行完成后，生成webmercator投影的出图框。
出图框为shp文件，由多边形构成，每一个多边形规定投影转换后单幅图像的边界，其属性name对应的属性值给定了输出文件名。影像分辨率设置为19m，对应网络切片地图第13级。
投影转换后的文件为：
在QGIS打开效果为：
还是单波段黑白图。
5.渲染
制作色带。
全球陆地高程约-200-9000m，以此为高程范围，制作色带。
color_table.txt-10 48 84 1500 32 55 1001 65 130 70100 85 150 70250 100 160 75500 120 165 85750 140 175 1001000 155 180 1151250 180 190 1251500 185 180 1201750 170 160 1102000 155 140 1002225 150 120 902500 130 105 702750 115 85 603000 220 220 2204000 254 254 254
渲染代码：
from osgeo import gdal,osrimport osimport globimport numpy as npimport tracebackfrom image_dehaze import *path = 'N:\\Globel_DEM\\SRTM_Global_webmercator\\'outpath_color_hillshade = 'N:\\Globel_DEM\\SRTM_Global_webmercator_color_hillshade\\'if not os.path.exists(outpath_color_hillshade): try: os.mkdir(outpath_color_hillshade) except: Noneoutpath_hillshade = os.path.join(outpath_color_hillshade,'hillshade')if not os.path.exists(outpath_hillshade): try: os.mkdir(outpath_hillshade) except: Noneoutpath_color_relief = os.path.join(outpath_color_hillshade,'color_relief')if not os.path.exists(outpath_color_relief): try: os.mkdir(outpath_color_relief) except: Noneall_files = glob.glob(path+'*.tif')i_s = 0total_number = len(all_files)is_build_overview = Truefor tif_file in all_files: i_s += 1 try: tif_file_out = tif_file#os.path.join(outpath_color_hillshade, os.path.basename(tif_file)) tif_file_out_color_hillshade = os.path.join(outpath_color_hillshade, os.path.basename(tif_file_out)) tif_file_out_hillshade = os.path.join(outpath_hillshade, os.path.basename(tif_file_out)) tif_file_out_color_relief = os.path.join(outpath_color_relief, os.path.basename(tif_file_out)) tif_file_out_processed = tif_file_out_color_hillshade + '_____processed.txt' tif_file_out_processing = tif_file_out_color_hillshade + '_____processing.txt' print(i_s,'---of---',total_number,tif_file_out) if os.path.exists(tif_file_out_processed): #print(tif_file_out_processed) continue if os.path.exists(tif_file_out_processing): #print(tif_file_out_processing) continue if os.path.exists(tif_file_out_color_hillshade): if not os.path.exists(tif_file_out_processed): with open(tif_file_out_processed,'w') as f: f.close() continue with open(tif_file_out_processing,'w') as f: f.close() in_ds = gdal.Open(tif_file_out) img_data = in_ds.GetRasterBand(int(1)).ReadAsArray() if (img_data.min() == img_data.max()): in_ds = None with open(tif_file_out_processing + '___________is_an_unvalid_file.txt','w') as f: f.close() continue if (0 == img_data.max()): in_ds = None with open(tif_file_out_processing + '___________is_an_unvalid_file.txt','w') as f: f.close() continue if not os.path.exists(tif_file_out_hillshade): gdal_dem_hillshade_code = 'gdaldem.exe hillshade ' + tif_file_out +' ' +tif_file_out_hillshade#D:\N07W012.tif D:\N07W012_shade.tif os.system(gdal_dem_hillshade_code) if not os.path.exists(tif_file_out_color_relief): gdal_dem_hillshade_code = 'gdaldem.exe color-relief ' + tif_file_out + ' ' + os.path.join(path,'global_dem_color_table.txt') + ' ' + tif_file_out_color_relief#D:\N07W012.tif D:\N07W012_shade.tif os.system(gdal_dem_hillshade_code) print(img_data.min(),img_data.max()) in_ds_hillshade = gdal.Open(tif_file_out_hillshade) in_ds_colorrelief = gdal.Open(tif_file_out_color_relief) xsize = in_ds.RasterXSize ysize = in_ds.RasterYSize in_geotrans = in_ds.GetGeoTransform() #print(in_geotrans) in_projection = in_ds.GetProjection() #print(in_projection) if len(in_projection) == 0: in_projection = None if in_geotrans[0] == 0 or in_geotrans[3]== 0: in_projection = None gtiff_driver = gdal.GetDriverByName('GTiff') create_options=['BIGTIFF=NO','TILED=YES','COMPRESS=LZW'] if xsize * ysize < 4000 * 4000: create_options=['BIGTIFF=NO'] img_bands = 3 out_ds = gtiff_driver.Create(tif_file_out_color_hillshade,xsize,ysize,img_bands,gdal.GDT_Byte,options = create_options) if in_geotrans is not None: try: #print('set geotrans:',tif_file_out_color_hillshade) out_ds.SetGeoTransform(in_geotrans) except: print('can not set geotrans:',tif_file_out_color_hillshade) if in_projection is not None: try: #print('set projection:',tif_file_out_color_hillshade) out_ds.SetProjection(in_projection) except: print('can not set projection:',tif_file_out_color_hillshade) img_data_hillshade = in_ds_hillshade.GetRasterBand(int(1)).ReadAsArray() #img_data = in_ds.GetRasterBand(int(1)).ReadAsArray() #print(img_data.min(),img_data.max()) nodata_value = -32768 nodata_area = GetMaskArea((img_data != nodata_value).astype(np.uint8),isExpand=False) water_area = GetMaskArea((img_data > 0).astype(np.uint8),isExpand=False) if water_area is not None: img_data_hillshade[water_area] = 255 if nodata_area is not None : img_data_hillshade[nodata_area] = 0 img_data_hillshade = img_data_hillshade.astype(np.float32)/255.0 for img_band in np.arange(img_bands): img_data_color_relief = in_ds_colorrelief.GetRasterBand(int(img_band+1)).ReadAsArray() #img_data = (img_data_hillshade * 0.6 + img_data_color_relief * 0.4).astype(np.uint8) img_data = (img_data_hillshade * img_data_color_relief) #img_data[img_data<1] = 1 img_data = img_data.astype(np.uint8) print(img_data.min(),img_data.max()) out_ds.GetRasterBand(int(img_band+1)).WriteArray(img_data) if is_build_overview or True: print('building overview...') #in_ds = gdal.Open(tif_file_out)#merge( names = None,format = 'GTiff',out_file = 'out.tif',outrange=None,psize_x = None,psize_y = None,nodata = None,a_nodata = None ): out_ds.BuildOverviews('average',[2,4,8,16,32,64,128]) out_ds = None in_ds = None in_ds_colorrelief = None in_ds_hillshade = None with open(tif_file_out_processed,'w') as f: f.close() try: os.remove(tif_file_out_processing) except: None except Exception as ex: print("have exception!") traceback.print_exc() break
运行结果为彩色渲染地形图，单幅打开是这样：
放大了看：
电脑性能好的话可以把全部影像拖进去，效果是这样：
亚洲腹部，一片巨大的高原，青藏高原。
珠穆朗玛峰。
以及非洲，规模宏大的乞力马扎罗山，一座巨大的火山。
性能再好的电脑，打开这么多影像，也会非常非常慢，最好的办法还是做切片。
10个月宝宝每天需要喝多少奶粉？
