如何用python画彩色（用python画彩虹）？如果你对这个不了解，来看看！
故宫下雪了！我用Python给它画了一组手绘图，仅用45秒（附代码），下面是数据派THU给大家的分享，一起来看看。
如何用python画彩色
来源：恋习Python
本文约1300字，建议阅读7分钟。
本文介绍图像处理的基本概念与原理，并手把手教你图片的手绘效果处理。
这几天，许多城市，迎来了2019年的第一场雪。
13日早晨，当北京市民拉开窗帘时发现，窗外雪花纷纷扬扬在空中飘落而且越下越大，树上、草地、屋顶、道路上...都落满雪花京城银装素裹，这是今冬以来北京迎来的第二场降雪一下雪，北京就变成了北平，故宫就变成了紫禁城八万张门票在雪花飘下来之前，便早已预订一空。
看着朋友圈、微博好友都在纷纷晒图，小编只能羡慕不已。
不过，我们突然想到，可以通过Python将故宫的建筑物图片，转化为手绘图（素描效果）。效果图如下：
一、概念与原理
我们都知道手绘图效果的特征主要有：
黑白灰色；边界线条较重；相同或相近色彩趋于白色；略有光源效果。核心原理：利用像素之间的梯度值和虚拟深度值对图像进行重构，根据灰度变化来模拟人类视觉的模拟程度。
把图像看成二维离散函数，灰度梯度其实就是这个二维离散函数的求导，用差分代替微分，求取图像的灰度梯度。常用的一些灰度梯度模板有：Roberts 梯度、Sobel 梯度、Prewitt 梯度、Laplacian 梯度。
以Sobel梯度计算来解释：
首先计算出
、
，然后计算梯度角
梯度方向及图像灰度增大的方向，其中梯度方向的梯度夹角大于平坦区域的梯度夹角。如下图所示，灰度值增加的方向梯度夹角大，此时梯度夹角大的方向为梯度方向。对应在图像中寻找某一点的梯度方向即通过计算该点与其8邻域点的梯度角，梯度角最大即为梯度方向。
二、图像的数组形式与变换
其中，需要用到的方法：
Image.open( ): 打开图片np.array( ) : 将图像转化为数组convert("L"): 将图片转换成二维灰度图片Image.fromarray( ): 将数组还原成图像uint8格式代码如下：
from PIL import Imageimport numpy as npim = Image.open(r"C:\Users\Administrator\Desktop\gugong\微信图片_20190216152248.jpg").convert('L')a=np.asarray(im).astype('float')print(a.shape,a.dtype)(1080, 608) float64#(1080, 608)分别表示高度，宽度三、图像的手绘效果处理
实现思路步骤：
梯度的重构构造guan光源效果梯度归一化图像生成1. 梯度的重构
numpy的梯度函数的介绍：
np.gradient(a) ： 计算数组a中元素的梯度，f为多维时，返回每个维度的梯度 离散梯度： xy坐标轴连续三个x轴坐标对应的y轴值：a, b, c 其中b的梯度是（c-a）/2 ，而c的梯度是： (c-b)/1当为二维数组时，np.gradient(a) 得出两个数组，第一个数组对应最外层维度的梯度，第二个数组对应第二层维度的梯度。
代码如下：
grad=np.gradient(a)grad_x,grad_y=gradgrad_x = grad_x * depth / 100.#对grad_x值进行归一化grad_y = grad_y * depth / 100.#对grad_y值进行归一化2. 构造guan光源效果
设计一个位于图像斜上方的虚拟光源光源相对于图像的视角为Elevation,方位角为Azimuth建立光源对各点梯度值的影响函数运算出各点的新像素值其中：
np.cos(evc.el) : 单位光线在地平面上的投射长度dx,dy,dz：光源对x，y，z三方向的影响程度3. 梯度归一化
构造x和y轴梯度的三维归一化单位坐标系；梯度与光源相互作用，将梯度转化为灰度。4. 图像生成
具体详情代码如下：
from PIL import Imageimport numpy as npimport osimport joinimport timedef image(sta,end,depths=10): a = np.asarray(Image.open(sta).convert('L')).astype('float') depth = depths # 深度的取值范围(0-100)，标准取10 grad = np.gradient(a) # 取图像灰度的梯度值 grad_x, grad_y = grad # 分别取横纵图像梯度值 grad_x = grad_x * depth / 100.#对grad_x值进行归一化 grad_y = grad_y * depth / 100.#对grad_y值进行归一化 A = np.sqrt(grad_x ** 2 + grad_y ** 2 + 1.) uni_x = grad_x / A uni_y = grad_y / A uni_z = 1. / A vec_el = np.pi / 2.2 # 光源的俯视角度，弧度值 vec_az = np.pi / 4. # 光源的方位角度，弧度值 dx = np.cos(vec_el) * np.cos(vec_az) # 光源对x 轴的影响 dy = np.cos(vec_el) * np.sin(vec_az) # 光源对y 轴的影响 dz = np.sin(vec_el) # 光源对z 轴的影响 b = 255 * (dx * uni_x + dy * uni_y + dz * uni_z) # 光源归一化 b = b.clip(0, 255) im = Image.fromarray(b.astype('uint8')) # 重构图像 im.save(end)def main(): xs=10 start_time = time.clock() startss = os.listdir(r"C:\Users\Administrator\Desktop\gugong") time.sleep(2) for starts in startss: start = ''.join(starts) sta = 'C:/Users/Administrator/Desktop/gugong/' + start end = 'C:/Users/Administrator/Desktop/gugong/' + 'HD_' + start image(sta=sta,end=end,depths=xs) end_time = time.clock() print('程序运行了 ----' + str(end_time - start_time) + ' 秒') time.sleep(3)main()程序运行了 ----43.01828205879955 秒 #一共35张图片最终效果图对比：
其他图片就不一一列举，若需要获取更多图片素材，后台回复“故宫雪景”即可获得；你也可以通过此代码为自己画一张手绘图；也可以为自己的家乡或母校画。
最后温馨提示：
人生苦短，我用Python；
除了生娃，啥都能干！！
欢迎大家在留言处，留言自己曾经用Python做过哪些有意思的事！
参考资料
北京理工大学的嵩天老师的网络课程http://www.icourse163.org/learn/BIT-1001870002?tid=1001963001#/learn/announce编辑：黄继彦
校对：陈瑞清
— 完 —
关注清华-青岛数据科学研究院官方微信公众平台“THU数据派”及姊妹号“数据派THU”获取更多讲座福利及优质内容。
用python画彩虹
你认识GDAL吗？你了解GDAL吗？
GDAL是栅格和矢量地理空间数据格式的转换器库，由开源地理空间基金会以X / MIT样式的开源许可证发布。作为一个库，它为调用的应用程序提供了所有支持的格式的单个栅格抽象数据模型和单个矢量抽象数据模型。 它还带有用于数据转换和处理的各种有用的命令行实用程序。反正我不是很熟，我只知道它很酷、很快、而且很厉害，不过好像也有很多缺陷。很酷很快很厉害我就不多说了，网上的彩虹屁多的和牛毛一样。我主要想说说几个问题，虽然这些问题可能只有我会遇到（意思是可能我安装的姿势不太对）。
你要是为我花过一块钱（上一篇），你可能就会发现一点点问题。为什么我得单独加载一个WGS坐标的tif图呢？这就是我想说的问题所在了，你要是看不下去就跳到下一段吧。python的gdal库可以读取很多的遥感影像格式，但是读进来并不意味着认出来，就好比上回的OMI，读进来后并不能准确的认出它的ds.GetGeoTransform()中的6个参数（分别是左下角x，左上角y，x和y的长宽以及两个0值---如果有错，请无情的批判我），好像也认不到它的坐标系？具体我忘了。总之，是有些问题的。
但是，今天的内容并不集中在python的gdal库，而是正！经！ GDAL 的C语言程序上，通过python的subprocess包调用cmd函数。起码这个方式上，我没有遇到那么多问题~
好了~废话不要那么多，不如直接上干货。。。
等等，上干货之前，你总要了解一下GDAL的安装吧，我觉得安装GDAL也是满满的血坑，但是。。。我就是不告诉你我摔了个怎么样的大趴。。。毕竟：
跳过GDAL安装，愉快的lu python代码吧
import subprocess#首先是拼接subprocess.call('D:/Anaconda3/python.exe'+' '+'gdal_merge.py '+'-of GTiff '+'-o '+str(oup)+' '+str(inp1))#这里，out是输出的路径+文件名，inp1是输入的subprocess.call("gdalwarp -cutline " + str(shp) + " -crop_to_cutline " + str(inp) + " -tr " + str(res) + " " + str(res) + " " + str(out), shell = True)#和上面一样的，out是输出的路径+文件名，inp1是输入的，两个res分别是空间分辨率的长和宽，shp是裁剪用的矢量哦你们要是问我为什么写的这么简单，解释的这么不酷？
还能为什么？我懒啊！
还在等什么？还不赶紧跨过GDAL安装的门槛，愉快的用GDAL快速的裁剪、拼接成千上万的遥感影像？
放两个结果给大家：
第一个是裁剪：
（猜猜这是哪？）
下一个是拼接：
（这个你们可能看不清，它是两个MAIAC的影像拼接而成，分为上下两部分）
佛系关注公众号：一个有趣的灵魂W
带给你不一样的惊喜
最后，希望我能和大家一样的幸福快乐~