站内推广方式有哪些,wordpress 插件安装在哪个文件夹,有哪些做短租的网站好,杭州专门做网站这个是之前写的代码了#xff0c;正好今天有空所以就和大家分享一下。我们在处理项目时#xff0c;有时候需要高清底图作为辅助数据源去对比数据#xff0c;所以可能会需要卫星数据。所以今天就和大家分享一下如何使用Python基于矢量范围批量下载高清遥感瓦片数据。 1 读取矢… 这个是之前写的代码了正好今天有空所以就和大家分享一下。我们在处理项目时有时候需要高清底图作为辅助数据源去对比数据所以可能会需要卫星数据。所以今天就和大家分享一下如何使用Python基于矢量范围批量下载高清遥感瓦片数据。 1 读取矢量边界 这里我们使用osgeo中的osr、ogr库去读取矢量的地理范围。之前也分享过感兴趣的可以去PythonGIS专栏里面看一看。注意我这里只算了四至范围而不是整个矢量的边界范围因为边界范围计算行列数不太好计算。正因如此我们下载的瓦片数据也不具备坐标系所以也裁剪不了。
def Open_Vector(path_shp)::param path_shp: 输入84坐标矢量:return: 返回四至范围ds ogr.Open(path_shp, True)# True表示以读写方式打开layer ds.GetLayer(0)# 获取图层feature layer.GetFeature(0)geom feature.GetGeometryRef()# 获取该要素的地理空间范围left, right, down, up geom.GetEnvelope()# 获取图层的地理范围return left, right, down, up
2 通过经纬度计算航带数 这里没什么好说的就是基础的公式直接计算行列数即可。这个函数在整个函数作为辅助函数主程序会自己调用它。
def calculation_tile(lat, lon, zoom)::param lat: 84坐标纬度:param lon: 84坐标经度:param zoom: 缩放级别:return: 瓦片的行列号# 将经纬度从WGS84坐标系转换为GCJ02坐标系# lon_deg,lat_deg WGS84_To_GCJ02(lon_deg,lat_deg)# 根据缩放级别计算格网数量n 2.0 ** zoom# 将纬度从度转换为弧度lat_radio math.radians(lat)# 计算瓦片中的x坐标tile_x int((lon 180.0) / 360.0 * n)# 计算瓦片中的y坐标tile_y int((1.0 - math.log(math.tan(lat_radio) (1 / math.cos(lat_radio))) / math.pi) / 2.0 * n)# 返回计算得到的瓦片坐标行和列return tile_x, tile_y
3 获取瓦片下载链接 这里使用了基础的反爬虫方法随机调用请求头也可以自己添加或减少请求头。即使如此有时还会爬取失败感兴趣的可以自己改进一下反爬方法如果不会就多运行几次。
def Get_image(url, x, y):agents [Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.1 (KHTML, like Gecko) Chrome/13.0.782.24 Safari/535.1,Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/534.27 (KHTML, like Gecko) Chrome/12.0.712.0 Safari/534.27,Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/60.0.3112.101 Safari/537.36,Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/532.5 (KHTML, like Gecko) Chrome/4.0.249.0 Safari/532.5,Mozilla/5.0 (Windows; U; Windows NT 5.2; en-US) AppleWebKit/532.9 (KHTML, like Gecko) Chrome/5.0.310.0 Safari/532.9,Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US) AppleWebKit/534.7 (KHTML, like Gecko) Chrome/7.0.514.0 Safari/534.7,Mozilla/5.0 (Windows; U; Windows NT 6.0; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/9.0.601.0 Safari/534.14,Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.14 (KHTML, like Gecko) Chrome/10.0.601.0 Safari/534.14,Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US) AppleWebKit/534.20 (KHTML, like Gecko) Chrome/11.0.672.2 Safari/534.20,]try:# 打印下载成功的消息显示瓦片的位置和下载状态print(瓦片 str(x) _ str(y) 下载成功)# 创建一个请求对象使用指定的URLrequests urllib.request.Request(url)# 为请求添加一个随机的User-Agent头以模拟不同的浏览器或客户端requests.add_header(User-Agent, random.choice(agents)) # 换用随机请求头# 使用指定的请求打开URL并设置超时时间为60秒image urllib.request.urlopen(requests, timeout60)# 读取返回的图像数据image_io image.read()# 使用BytesIO将图像数据转换为可处理的字节流对象image_bytes io.BytesIO(image_io)# 使用PIL库打开图像image Image.open(image_bytes)# 将图像从RGB格式转换为BGR格式OpenCV需要的格式image cv2.cvtColor(np.asarray(image), cv2.COLOR_RGB2BGR)except EOFError:# 如果发生EOFError可能由于网络问题、超时等打印下载失败的消息并尝试重试print(瓦片 str(x) _ str(y) 下载失败,正在重试......)Get_image(url, x, y) # 递归调用Get_image函数进行重试# 返回处理后的图像数据return image
4 主程序 这里就不过多解释了我的代码注释非常完善如果有什么不懂的直接留言即可。整体就是筛选出最大的缩放级别再调用行列数函数计算出瓦片的编号再调运瓦片的下载链接最后拼接起来即可。
# -*- coding: utf-8 -*-Time 2023/4/9 14:37
Auth RS迷途小书童
File Vector Download Remote Sensing Tile Data.py
IDE PyCharm
Purpose根据矢量范围下载三方地图瓦片
Web博客地址:https://blog.csdn.net/m0_56729804import io
import cv2
import math
import random
import numpy as np
from osgeo import ogr
import urllib.request
from PIL import Imagedef Write_image(lat1, lon1, lat2, lon2, out_path)::param lat1: 左上角纬度:param lon1: 左上角经度:param lat2: 右下角纬度:param lon2: 右下角经度:return: 返回瓦片影像zooms list()# 创建一个空列表zooms用于存储所有的缩放级别for i in range(1, 19):# 循环缩放级别col calculation_tile(lat1, lon1, i)# 将经纬度转换为对应的地图瓦片编号结果存储在col中row calculation_tile(lat2, lon2, i)if col[0] - row[0] 0 or col[1] - row[1] 0:continueelse:zooms.append(i)# 如果差值不为0将当前的缩放级别i添加到zooms列表中zoom zooms[-1]# 获取zooms列表中的最后一个元素即最大的缩放级别并存储在zoom变量中left_up calculation_tile(lat1, lon1, zoom)# 使用最大的缩放级别和第一个经纬度范围调用函数获取左上角的地图瓦片编号存储在left_up中right_down calculation_tile(lat2, lon2, zoom)# 使用最大的缩放级别和第二个经纬度范围调用函数获取右下角的地图瓦片编号存储在right_down中images_columns list()# 创建一个空列表images_columns用于存储所有的地图瓦片图像列print(当前瓦片行数, right_down[0]-left_up[0])print(当前瓦片列数, right_down[1] - left_up[1])print(--------------------------------------数据获取--------------------------------------)for x in range(left_up[0], right_down[0]):# 循环行images_rows list()# 创建一个空列表images_rows用于存储所有的地图瓦片图像行for y in range(left_up[1], right_down[1]):# 循环列tile_url http://t4.tianditu.com/DataServer?Timg_wxstr(x)ystr(y)lstr(zoom) \tk45c78b2bc2ecfa2b35a3e4e454ada5ceimage Get_image(tile_url, x, y)cv2.imwrite(out_path /%s.jpg % (str(x)_str(y)), image)images_rows.append(image)# 将获取到的瓦片图像添加到images_rows列表中用于后续的图像合成img_column_new np.vstack(images_rows)# 使用NumPy的v stack函数将images_rows列表中的所有图像竖直堆叠起来形成一个新的图像列images_columns.append(img_column_new)# 将这个新的图像列添加到images_columns列表中用于后续的图像合成print(正在拼接瓦片数据......)result np.hstack(images_columns)# 使用NumPy的h stack函数将images_columns列表中的所有图像水平堆叠起来形成一个最终的大图像print(正在保存瓦片数据......)cv2.imwrite(out_path /result.jpg, result)return result
5 总结
下列有一些三方底图链接也可以查看【PythonGIS】第三方地图服务汇总。 tile_url http://www.google.cn/maps/vt/pb!1m4!1m3!1istr(zoom)!2istr(x)!3istr(y)!2m3!1e0!2sm!3i345013117!3m8!2szh-CN!3scn!5e1105!12m4!1e68!2m2!1sset!2sRoadmap!4e0
# Google地图瓦片
tile_url http://mt3.google.cn/vt/lyrss110hlzh-CNglcnsrcappxstr(x)ystr(y)zstr(zoom)sG
# Google影像瓦片
tile_url http://t4.tianditu.com/DataServer?Timg_wxstr(x)ystr(y)lstr(zoom)tk45c78b2bc2ecfa2b35a3e4e454ada5ce
# 天地图卫星数据vec_w电子地图2000坐标系
http://wprd01.is.autonavi.com/appmaptile?langzh_cnsize1scl1style6x str(x) y str(y) z str(zoom) ltype3
# 高德底图偏移火星坐标系这里输入的矢量需要是WGS84坐标系的经纬度不能是投影坐标系哦。此外如果使用高德、百度等底图可能会有一定的偏移因为我国需要加密成火星坐标系但是还是可以用的略微有偏移而已作为对比图够用了。天地图就无所谓它的坐标是准的。
