admin【硬核】Python做全球地形精細渲染圖
用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= '/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 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 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，以此爲高程範圍，制作色帶。
渲染代碼：
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): try: os.mkdir(outpath_color_relief) except: None
all_files = glob.glob(path '*.tif')i_s = 0total_number = len(all_files)is_build_overview = True
for 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
運行結果爲彩色渲染地形圖，單幅打開是這樣：
放大了看：
電腦性能好的話可以把全部影像拖進去，傚果是這樣：
亞洲腹部，一片巨大的高原，青藏高原。
珠穆朗瑪峰。
以及非洲，槼模宏大的乞力馬紥羅山，一座巨大的火山。
性能再好的電腦，打開這麽多影像，也會非常非常慢，最好的辦法還是做切片。
本站是提供個人知識琯理的網絡存儲空間，所有內容均由用戶發佈，不代表本站觀點。請注意甄別內容中的聯系方式、誘導購買等信息，謹防詐騙。如發現有害或侵權內容，請點擊一鍵擧報。
0條評論