2017-04-12 1 views
0

J'ai tracé une carte avec le fichier de forme dans le fond de carte, sur la base de cela, je veux ajouter les données DEM (.tif) dans le carte. J'ai téléchargé un morceau de données .tif (la longitude et la latitude étaient vraiment à l'intérieur de la plage du fichier de formes) sur le site SRTM, mais finalement quand j'ai lancé le programme, il n'affichait pas les données raster dans la loction spécifique de la carteComment ajouter des données DEM (.tif) dans un fichier de forme dans gdal avec python

from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt 
import numpy as np 
from matplotlib.patches import Polygon 
from osgeo import gdal 
from numpy import linspace 
from numpy import meshgrid 

fig = plt.figure(figsize=(8,6), dpi=80) 
ax1 = fig.add_axes([0.1,0.1,1.0,1.0]) 
map = Basemap(llcrnrlon=114.7, 
      llcrnrlat=29.3, 
      urcrnrlon=120, 
      urcrnrlat=34.6, 
     resolution='h', projection='tmerc', lat_0 =31.5,lon_0=116.5,ax=ax1) 
shp_info = map.readshapefile("bou2_4l",'state',color='k',linewidth='1',drawbounds=True) 

ds = gdal.Open("srtm_60_06.tif") 
data = ds.ReadAsArray() 
loc3=[30,35] 
lat3=[115,120] 
x1,y1 = map(loc3,lat3) 
x = linspace(x1[0],x1[1], data.shape[1]) 
y = linspace(y1[0],y1[1], data.shape[0]) 
xx, yy = meshgrid(x, y) 
map.pcolormesh(xx, yy, data) 
plt.show() 

Répondre

0

Vous pouvez avoir une réponse plus concise si vous avez demandé cela à gis.stackexchange.

De toute façon, cela me semble une erreur de projection. Je suis censé être corrigé, mais les tiffs STRM sont des GeoTiffs donc ils ont l'information de projection dans l'en-tête.

Il semble que vous utilisiez Transverse Mercator alors que les données SRTM utilisent un système de coordonnées géographiques.

Essayez votre reprojeter TIFF, je pense que gdal.Warp devrait le faire:

gdal.Warp(output_raster,input_raster,dstSRS='EPSG:your_basemap_projection')