zoukankan      html  css  js  c++  java
  • 计算山体阴影

    #!/usr/bin/env python
    # -*- coding: utf-8 -*-
    from osgeo import gdal
    from numpy import gradient
    from numpy import pi
    from numpy import arctan
    from numpy import arctan2
    from numpy import sin
    from numpy import cos
    from numpy import sqrt
    import matplotlib.pyplot as plt
    import subprocess
    
    
    def hillshade(array, azimuth, angle_altitude):
        """
        :param array: input USGS ASCII DEM / CDED .dem
        :param azimuth: sun position
        :param angle_altitude: sun angle
        :return: numpy array
        """
        #计算梯度
         x, y = gradient(array)
        #计算坡度
         slope = pi/2. - arctan(sqrt(x*x + y*y))
        #计算坡向
         aspect = arctan2(-x, y)
        #计算方位角
         azimuthrad = azimuth * pi / 180.
        altituderad = angle_altitude * pi / 180.
    
    
        shaded = sin(altituderad) * sin(slope)
         + cos(altituderad) * cos(slope)
         * cos(azimuthrad - aspect)
        # return 255*(shaded + 1)/2
        return aspect
    
    ds = gdal.Open('../geodata/092j02_0200_demw.dem')
    arr = ds.ReadAsArray()
    
    hs_array = hillshade(arr, 90, 45)
    plt.imshow(hs_array,cmap='Greys')
    plt.savefig('../geodata/hillshade_whistler.png')
    plt.show()
    
    # gdal command line tool called gdaldem
    # link  http://www.gdal.org/gdaldem.html
    # usage:
    # gdaldem hillshade input_dem output_hillshade
    # [-z ZFactor (default=1)] [-s scale* (default=1)]"
    # [-az Azimuth (default=315)] [-alt Altitude (default=45)]
    # [-alg ZevenbergenThorne] [-combined]
    # [-compute_edges] [-b Band (default=1)] [-of format] [-co "NAME=VALUE"]* [-q]
    
    create_hillshade = '''gdaldem hillshade -az 315 -alt 45 ../geodata/092j02_0200_demw.dem ../geodata/hillshade.tif'''
    
    #通过标准库中的subprocess包来fork一个子进程,并运行一个外部的程序
    subprocess.call(create_hillshade)
  • 相关阅读:
    MarkDown学习记录
    AngularJS+Ionic开发-2.项目结构介绍
    【二分+贪心+倍增】【NOIP2012】疫情控制
    【考试总结】NOIP模拟 test11-1
    【考试总结】NOIP模拟 test10-27
    【考试总结】AHSOFNU&QZQZ test10-25
    【差分+二分答案】【NOIP2012】借教室
    【动态规划】【NOIP2015】子串
    【高精度+贪心】【NOIP2012】国王游戏
    ICPC Central Russia Regional Contest (CRRC 19)题解
  • 原文地址:https://www.cnblogs.com/gispathfinder/p/5790130.html
Copyright © 2011-2022 走看看