代码公开，用python做全球精细地形渲染图

小猿猴GISer · 公众号 · · 2025-02-18 19:19

正文

用python编程做一版全球地形渲染图。

1. 下载数据

利用python批量下载全球30m分辨率地形数据。

代码：

#!/usr/bin/python# _*_ coding: utf-8 _*_
from urllib.request import urlopenfrom urllib.request import Requestimport urllibimport reimport osimport sys
http_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 False
inpath = '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 = 0
for 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 glob
src_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 None
outfile_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 gdal
def 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, lat    def 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 = None
def 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 = True
dst_Res = 20037508.342787001 * 2 / 256 / (2 ** ge_level)psize_x = dst_Respsize_y = -dst_Res
wkt_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_ds
table_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:        None
outpath_hillshade = os.path.join(outpath_color_hillshade,'hillshade')if not os.path.exists(outpath_hillshade):    try:        os.mkdir(outpath_hillshade)    except:        None
outpath_color_relief = os.path.join(outpath_color_hillshade,'color_relief')if not os.path.exists(outpath_color_relief):

代码公开，用python做全球精细地形渲染图

正文

请到「今天看啥」查看全文