精华内容
下载资源
问答
  • 好吧,我洋洋洒洒,1分钟不到,就把切片服务调出来了。附录:代码。抄代码前,点个赞,能活999! 注意下面代码是离线环境的(所以你要在open layer官网搞个ol.css和ol.js)要是不想找,也找不到,可以下载我准备好的...

    啊,最近遇到一个扯皮的单位,一些事情,

    坑一  :一个奇葩搞项目开发开发商,非要说俺们发的服务有问题。……(其实可以用arcgis api 调用相关服务看,arcgis api 调用没问题)

    坑二:非要俺用openlayer4调用地图服务给他们看,其实他们不会openlayer。找arcgis 售后,售后说openlayer非arcgis产品,概不支持。

    好吧,

    古人云——工作就是填坑的过程。

    古人又云——老板请你来,就是请你用你的才华,把坑填了。

    吃完一顿饭,再外面逛了一圈,那个开发商还没用openlayer把服务调出来。

    好吧,我洋洋洒洒,1分钟不到,就把切片服务调出来了。附录:代码。抄代码前,点个赞,能活999!

    注意下面代码是离线环境的(所以你要在open layer官网搞个ol.css和ol.js)要是不想找,也找不到,可以下载我准备好的文件,代码!!哇,送代码啦!https://download.csdn.net/download/ucs426/13122436

    <!DOCTYPE html>
    <html>
      <head>
        <title>XYZ Esri</title>
        <link rel="stylesheet" href="./ol.css" type="text/css">
        <!-- The line below is only needed for old environments like Internet Explorer and Android 4.x -->
        <!-- <script src="https://cdn.polyfill.io/v2/polyfill.min.js?features=requestAnimationFrame,Element.prototype.classList,URL"></script> -->
        <script src="./ol.js"></script>
      </head>
      <body>
        <div id="map" class="map"></div>
        <script>
          var map = new ol.Map({
            target: 'map',
            layers: [
              new ol.layer.Tile({
                source: new ol.source.XYZ({
                  attributions: 'Tiles © <a href="https://xxxxxxx:6443/arcgis/rest/services/NC/MapServer">ArcGIS</a>',
                  url: 'https://xxxxxxx:6443/arcgis/rest/services/NC/MapServer/tile/{z}/{y}/{x}'
                })
              })
            ],
            view: new ol.View({
              center: ol.proj.fromLonLat([115.85129,28.703333]),
              zoom: 11
            })
          });
        </script>
      </body>
    </html>
    欢迎关注我的微信公众号!

     

     

    展开全文
  • 现在想在web前端用JavaScript 加载本地的map imagelayer.下图是在portal里的情况。 通过代码 <script> require([ "esri/Map", "esri/views/SceneView... ], function(Map, SceneView, MapImageLayer) {

    现在想在web前端用JavaScript 加载本地的map imagelayer.下图是在portal里的情况。

    通过代码

     <script>
            require([
              "esri/Map",
              "esri/views/SceneView",
              "esri/layers/MapImageLayer"
            ], function(Map, SceneView, MapImageLayer) {
              /*****************************************************************
              
               *****************************************************************/
              var permitsLayer = new MapImageLayer({
                portalItem: {
                  // autocasts as new PortalItem()
                  id: "27221c4efe7e417093e0110a75cb71b7"
                }
              });
      
              /*****************************************************************
               * Add the layer to a map
               *****************************************************************/
              var map = new Map({
                basemap: "",
                layers: [permitsLayer]
              });
      
              var view = new SceneView({
                container: "viewDiv",
                map: map
              });
      
      
            });
          </script>

    效果如下所示

    arcgis for javascript系列博客

    arcgis for javascript 示例(一)如何调用本地portal 发布的webmap

    arcgis for javascript 示例(二)如何调用本地portal 里的 webscene

     

    展开全文
  • 遥感影像切分切片

    千次阅读 2019-11-19 13:57:29
    文章目录遥感影像切片切片拼接 遥感影像切片 生活当中,我们可能经常遇到处理一个很大的遥感影像情况,例如做一些逐像元运算,不便于并行处理。因此制作了将影像切分成多个小片的程序,切分后保持原有的数值、数据...

    遥感影像切片

    生活当中,我们可能经常遇到处理一个很大的遥感影像情况,例如做一些逐像元运算,不便于并行处理。因此制作了将影像切分成多个小片的程序,切分后保持原有的数值、数据类型、波段数、投影。

    切片

    用法为
    python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数
    例如:
    python ./splitImage.py aa.tif ./output 8 1024

    splitImage.py内容如下

    import numpy as np
    import gdal
    import os
    import sys
    from  multiprocessing import Pool
    
    NP2GDAL_CONVERSION = {
      "uint8": 1,
      "int8": 1,
      "uint16": 2,
      "int16": 3,
      "uint32": 4,
      "int32": 5,
      "float32": 6,
      "float64": 7,
      "complex64": 10,
      "complex128": 11,
    }
    def split(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata):
        dataset = gdal.Open(rasterpath)
        if ((i+1)*size > cols) | ((j+1)*size>rows):
            #x向越界
            if ((i+1)*size > cols) & ((j+1)*size<=rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,size)
            #y向越界
            elif ((i+1)*size <= cols) & ((j+1)*size>rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,rows-j*size)
            #xy方向均越界
            else:
                print(cols-i*size,rows-j*size)
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,rows-j*size)
        else:
            data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,size)
        #如果第一个波段的最大值等于最小值,认为是无效值,不对其创建分片
        if data0.max() == data0.min():
            return
    
    
        format = "GTiff"
        driver = gdal.GetDriverByName(format)
        bandsNum = dataset.RasterCount
        basename = os.path.basename(rasterpath)
        newbasename = basename[:basename.rfind(".")]
        dirpath = os.path.dirname(rasterpath)
        dst_ds = driver.Create(outpath + newbasename+'_%s%s.tif' % (dex[i], dex[j]), size, size, bandsNum, gdaltype)
        gtnew = (gt[0] + i * size * gt[1], gt[1], gt[2], gt[3] + j * size * gt[5], gt[4], gt[5])
        dst_ds.SetGeoTransform(gtnew)
        dst_ds.SetProjection(proj)
        for k in range(bandsNum):
            if ((i + 1) * size > cols) | ((j + 1) * size > rows):
                # x向越界
                if ((i + 1) * size > cols) & ((j + 1) * size <= rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, size)
                # y向越界
                elif ((i + 1) * size <= cols) & ((j + 1) * size > rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, size, rows - j * size)
                # xy方向均越界
                else:
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, rows - j * size)
                smally,smallx=data2.shape
                if nodata==None:
                    data1=np.zeros((size,size),dtype=datatype)
                else:
                    data1 = np.ones((size, size), dtype=datatype)*nodata
                data1[0:smally,0:smallx]=data2
    
            else:
                data1=dataset.GetRasterBand(k+1).ReadAsArray(i*size,j*size,size,size)
    
            dst_ds.GetRasterBand(k+1).WriteArray(data1)
            if nodata != None:
                dst_ds.GetRasterBand(k + 1).SetNoDataValue(nodata)
        dataset=None
        dst_ds=None
    
    if __name__=="__main__":
        args=sys.argv
        try:
            rasterpath = args[1]
            outpath = args[2]
            parelle = int(args[3])
            # 分片的像元行列数,输出为正方形的
            size = int(args[4])
        except:
            print("""Usage:python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数\nUsage:python ./splitImage.py  aa.tif ./output 8 1024""")
            raise
    
        # rasterpath = "H:\\testsplit\\myimage.tif"
        # outpath = 'H:\\testsplit\\test4'
        # parelle=6
        # # 分片的像元行列数,输出为正方形的
        # size = 2048
        # 区分每块的位置,000到999
        dex = ["%.3d" % i for i in range(1000)]
    
        if not os.path.exists(outpath):
            os.makedirs(outpath)
        if not outpath.endswith(os.path.sep):
            outpath = outpath + os.path.sep
        dataset=gdal.Open(rasterpath)
        proj=dataset.GetProjection()
        gt=dataset.GetGeoTransform()
        datatype=dataset.GetRasterBand(1).ReadAsArray(0,0,1,1).dtype.name
        #获取无效值,认为每个波段的无效值是相同的
        nodata=dataset.GetRasterBand(1).GetNoDataValue()
        #numpy的数据类型,转换为gdal的数据类型
        gdaltype=NP2GDAL_CONVERSION[datatype]
        #总行列数
        cols,rows=dataset.RasterXSize,dataset.RasterYSize
        dataset=None
        #获取分片的行数和列数
        numx,numy=int(np.ceil(cols/size)),int(np.ceil(rows/size))
        pool = Pool(parelle)
        for i in range(numx):
            for j in range(numy):
                pool.apply_async(split, args=(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata))
        pool.close()
        pool.join()
    

    切分后的文件列表如下

    .
    ├── myimage_000000.tif
    ├── myimage_000001.tif
    ├── myimage_000002.tif
    ├── myimage_000003.tif
    ├── myimage_000004.tif
    ├── myimage_000005.tif
    ├── myimage_001000.tif
    ├── myimage_001001.tif
    ├── myimage_001002.tif
    ├── myimage_001003.tif
    ├── myimage_001004.tif
    ├── myimage_001005.tif
    ├── myimage_002000.tif
    ├── myimage_002001.tif
    ├── myimage_002002.tif
    ├── myimage_002003.tif
    ├── myimage_002004.tif
    ├── myimage_002005.tif
    ├── myimage_003000.tif
    ├── myimage_003001.tif
    └── myimage_003002.tif

    多进程版

    import numpy as np
    import gdal
    import os
    import sys
    from  multiprocessing import Pool
    
    NP2GDAL_CONVERSION = {
      "uint8": 1,
      "int8": 1,
      "uint16": 2,
      "int16": 3,
      "uint32": 4,
      "int32": 5,
      "float32": 6,
      "float64": 7,
      "complex64": 10,
      "complex128": 11,
    }
    def split(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata):
        dataset = gdal.Open(rasterpath)
        if ((i+1)*size > cols) | ((j+1)*size>rows):
            #x向越界
            if ((i+1)*size > cols) & ((j+1)*size<=rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,size)
            #y向越界
            elif ((i+1)*size <= cols) & ((j+1)*size>rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,rows-j*size)
            #xy方向均越界
            else:
                print(cols-i*size,rows-j*size)
                data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,rows-j*size)
        else:
            data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,size)
        #如果第一个波段的最大值等于最小值,认为是无效值,不对其创建分片
        if data0.max() == data0.min():
            return
        #S表示分片的意思
    
        format = "GTiff"
        driver = gdal.GetDriverByName(format)
        bandsNum = dataset.RasterCount
        basename = os.path.basename(rasterpath)
        newbasename = basename[:basename.rfind(".")]
        dirpath = os.path.dirname(rasterpath)
        dst_ds = driver.Create(outpath + newbasename+'_%s%s.tif' % (dex[i], dex[j]), size, size, bandsNum, gdaltype)
        gtnew = (gt[0] + i * size * gt[1], gt[1], gt[2], gt[3] + j * size * gt[5], gt[4], gt[5])
        dst_ds.SetGeoTransform(gtnew)
        dst_ds.SetProjection(proj)
        for k in range(bandsNum):
            if ((i + 1) * size > cols) | ((j + 1) * size > rows):
                # x向越界
                if ((i + 1) * size > cols) & ((j + 1) * size <= rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, size)
                # y向越界
                elif ((i + 1) * size <= cols) & ((j + 1) * size > rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, size, rows - j * size)
                # xy方向均越界
                else:
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, rows - j * size)
                smally,smallx=data2.shape
                if nodata==None:
                    data1=np.zeros((size,size),dtype=datatype)
                else:
                    data1 = np.ones((size, size), dtype=datatype)*nodata
                data1[0:smally,0:smallx]=data2
    
            else:
                data1=dataset.GetRasterBand(k+1).ReadAsArray(i*size,j*size,size,size)
    
            dst_ds.GetRasterBand(k+1).WriteArray(data1)
            if nodata != None:
                dst_ds.GetRasterBand(k + 1).SetNoDataValue(nodata)
        dataset=None
        dst_ds=None
    
    if __name__=="__main__":
        args=sys.argv
        try:
            rasterpath = args[1]
            outpath = args[2]
            parelle = int(args[3])
            # 分片的像元行列数,输出为正方形的
            size = int(args[4])
        except:
            print("""Usage:python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数\nUsage:python ./splitImage.py  aa.tif ./output 8 1024""")
            raise
    
        # rasterpath = "H:\\testsplit\\myimage.tif"
        # outpath = 'H:\\testsplit\\test4'
        # parelle=6
        # # 分片的像元行列数,输出为正方形的
        # size = 2048
        # 区分每块的位置,000到999
        dex = ["%.3d" % i for i in range(1000)]
    
        if not os.path.exists(outpath):
            os.makedirs(outpath)
        if not outpath.endswith(os.path.sep):
            outpath = outpath + os.path.sep
        dataset=gdal.Open(rasterpath)
        proj=dataset.GetProjection()
        gt=dataset.GetGeoTransform()
        datatype=dataset.GetRasterBand(1).ReadAsArray(0,0,1,1).dtype.name
        #获取无效值,认为每个波段的无效值是相同的
        nodata=dataset.GetRasterBand(1).GetNoDataValue()
        #numpy的数据类型,转换为gdal的数据类型
        gdaltype=NP2GDAL_CONVERSION[datatype]
        #总行列数
        cols,rows=dataset.RasterXSize,dataset.RasterYSize
        dataset=None
        #获取分片的行数和列数
        numx,numy=int(np.ceil(cols/size)),int(np.ceil(rows/size))
        pool = Pool(parelle)
        for i in range(numx):
            for j in range(numy):
                pool.apply_async(split, args=(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata))
        pool.close()
        pool.join()
    

    带一定的重合度(多进程版)

    • 相邻的切图之间,包括上下和左右,重合一定数量的像元,防止拼接时候出现一些裂缝等现象
    • 此方法更加通用
    import numpy as np
    import gdal
    import os
    import sys
    from  multiprocessing import Pool
    
    NP2GDAL_CONVERSION = {
      "uint8": 1,
      "int8": 1,
      "uint16": 2,
      "int16": 3,
      "uint32": 4,
      "int32": 5,
      "float32": 6,
      "float64": 7,
      "complex64": 10,
      "complex128": 11,
    }
    def split(rasterpath,outpath,dex,i, j,size,lap,cols,rows,gdaltype,datatype,gt,proj,nodata):
        dataset = gdal.Open(rasterpath)
        if (i+size > cols) | (j+size>rows):
            #x向越界
            if (i+size > cols) & (j+size<=rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i,j, cols-i,size)
            #y向越界
            elif (i+size <= cols) & (j+size>rows):
                data0 = dataset.GetRasterBand(1).ReadAsArray(i , j, size,rows-j)
            #xy方向均越界
            else:
                #print(cols-i*size,rows-j*size)
                data0 = dataset.GetRasterBand(1).ReadAsArray(i , j , cols-i,rows-j)
        else:
            data0 = dataset.GetRasterBand(1).ReadAsArray(i, j , size,size)
        #如果第一个波段的最大值等于最小值,认为是无效值,不对其创建分片
        if data0.max() == data0.min():
            return
    
        format = "GTiff"
        driver = gdal.GetDriverByName(format)
        bandsNum = dataset.RasterCount
        basename = os.path.basename(rasterpath)
        newbasename = basename[:basename.rfind(".")]
        dirpath = os.path.dirname(rasterpath)
        dst_ds = driver.Create(outpath + newbasename+'_%s%s.tif' % (dex[int(i/(size-lap))], dex[int(j/(size-lap))]), size, size, bandsNum, gdaltype)
        gtnew = (gt[0] + i  * gt[1], gt[1], gt[2], gt[3] + j * gt[5], gt[4], gt[5])
        dst_ds.SetGeoTransform(gtnew)
        dst_ds.SetProjection(proj)
        for k in range(bandsNum):
            if (i+size > cols) | (j+size>rows):
                # x向越界
                if (i+size > cols) & (j+size<=rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i , j , cols - i , size)
                # y向越界
                elif (i+size <= cols) & (j+size>rows):
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i , j , size, rows - j )
                # xy方向均越界
                else:
                    data2 = dataset.GetRasterBand(k+1).ReadAsArray(i , j , cols - i , rows - j )
                smally,smallx=data2.shape
                if nodata==None:
                    data1=np.zeros((size,size),dtype=datatype)
                else:
                    data1 = np.ones((size, size), dtype=datatype)*nodata
                data1[0:smally,0:smallx]=data2
    
            else:
                data1=dataset.GetRasterBand(k+1).ReadAsArray(i,j,size,size)
    
            dst_ds.GetRasterBand(k+1).WriteArray(data1)
            if nodata != None:
                dst_ds.GetRasterBand(k + 1).SetNoDataValue(nodata)
        dataset=None
        dst_ds=None
    
    if __name__=="__main__":
        args=sys.argv
        try:
            rasterpath = args[1]
            outpath = args[2]
            parelle = int(args[3])
            # 分片的像元行列数,输出为正方形的
            size = int(args[4])
            lap=int(args[5])
        except:
            print("""Usage:python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数 重叠像元数\nUsage:python ./splitImage.py  aa.tif ./output 8 1024 10""")
            raise
    
    
        # 区分每块的位置,000到999
        dex = ["%.3d" % i for i in range(1000)]
    
        if not os.path.exists(outpath):
            os.makedirs(outpath)
        if not outpath.endswith(os.path.sep):
            outpath = outpath + os.path.sep
        dataset=gdal.Open(rasterpath)
        proj=dataset.GetProjection()
        gt=dataset.GetGeoTransform()
        datatype=dataset.GetRasterBand(1).ReadAsArray(0,0,1,1).dtype.name
        #获取无效值,认为每个波段的无效值是相同的
        nodata=dataset.GetRasterBand(1).GetNoDataValue()
        #numpy的数据类型,转换为gdal的数据类型
        gdaltype=NP2GDAL_CONVERSION[datatype]
        #总行列数
        cols,rows=dataset.RasterXSize,dataset.RasterYSize
        dataset=None
        #获取分片的行数和列数
        numx=int(np.ceil((cols-size)/(size-lap)+1))
        numy=int(np.ceil((rows-size)/(size-lap)+1))
    
        #numx,numy=int(np.ceil(cols/size)),int(np.ceil(rows/size))
    
        pool = Pool(parelle)
        numxs=[(size-lap)*(i-1) for i in range(1,numx+1)]
        numys = [(size - lap) * (i - 1) for i in range(1, numy + 1)]
        #直接传入左上角的索引
        for i in numxs:
            for j in numys:
                pool.apply_async(split, args=(rasterpath,outpath,dex,i, j,size,lap,cols,rows,gdaltype,datatype,gt,proj,nodata))
        pool.close()
        pool.join()
    
    

    拼接

    gdal_merge.py拼接

    上一步拆分完了,做一些处理,最后可能还需拼接成一个大的影像。这里比较简单,直接调用gdal自带的gdal_merge.py拼接程序。如果电脑安装了python的gdal库的话就会有这个文件,可以通过搜索查找。
    使用前首先将gdal_merge.py复制到当前目录下。
    这里将一个文件夹下的所有tif文件拼接成一个文件。 具体可以修改tif变量。
    用法为
    python 此文件的路径 输入文件夹 输出tif nodata值
    例如:
    python ./mergeImage.py ./output /out.tif 0
    mergeImage.py文件内容

    import os
    import sys
    args=sys.argv
    if len(args)!=4:
        print("usage:python thispy folder_path outtif nodata")
        sys.exit(1)
    path=args[1]
    outtif=args[2]
    nodata=eval(args[3])
    tifs=[os.path.join(path,i) for i in os.listdir(path) if i.endswith(".tif")]
    #print(tifs)
    #将 gdal_merge.py复制到当前目录下
    os.system("python gdal_merge.py -init %s -n %s -a_nodata %s -o %s %s"%(nodata,nodata,nodata,outtif," ".join(tifs)))
    
    

    文件数量过多导致的命令太长的问题

    gdal拼接的代码修改版,原有的自带的gdal_merge.py如果输入文件名称加起来的长度太长的话容易导致报错,因为cmd的命令长度有限制,修改后解决这个问题
    gdal_merge_my.py的地址

    python gdal_merge_my.py -init 0 -n 0 -a_nodata 0 -o 输出影像名称  输出路径\*.tif
    
    展开全文
  • 不管是商业软件还是开源软件发布的影像瓦片服务,本质都是将遥感影像按照不同的级别(不同比例尺)切片成一张张小的瓦片地图,在浏览器上查看的时候根据比例尺和范围查询相应的瓦片地图进行加载。但是发布在互联网上...

    0 前言

    随着传感器的精度不断提高,遥感影像所占的空间也越来越大,为了方便分享浏览,通常会发布成影像瓦片服务。不管是商业软件还是开源软件发布的影像瓦片服务,本质都是将遥感影像按照不同的级别(不同比例尺)切片成一张张小的瓦片地图,在浏览器上查看的时候根据比例尺和范围查询相应的瓦片地图进行加载。但是发布在互联网上的影像瓦片地图无法进行本地化的分析处理,这个时候就需要将瓦片地图抓取到本地进行拼接处理。

    1 抓取瓦片图的原理

    瓦片地图服务分享在互联网或局域网上,遵循网络传输协议和特定的加载规则。平时我们在地图浏览页面浏览影像的操作,比如放大缩小和平移,都是以地图窗口为视角作为查询范围,结合当前地图的缩放级别,来查询特定级别特定范围下的瓦片数据。然后瓦片地图就像拼图一样加载在浏览器页面上,形成我们看到的整幅地图的效果。要想抓取瓦片地图,首先要了解瓦片地图服务的构造。本文以智慧广州公众服务“天地图·广州 智慧广州时空信息云平台”为例,进行讲解。广州天地图主页展示如下:
    天地图1
    因为需要使用到浏览器的开发者工具,在这里推荐使用谷歌浏览器。在谷歌浏览器的界面下按F12调出开发者工具,切换到Network选项卡,然后在地图窗口内放大一下地图来请求新的瓦片,请求的瓦片信息会在Network选项卡中显示出来,如下:
    天地图2

    我们选中其中一个瓦片图,查看Request URL的内容,复制出来观察一下。

    http://gzmap.gov.cn:12345/ServiceAccess/MapService-T/%E5%BD%B1%E5%83%8F/1e85e146b0728d2fb3a5312c75089400/tile/2/377/284

    可以看到URL后面的“tile/2/377/284”表示当前瓦片的级别和行列号,该瓦片处在影像级别2,行号377,列号284的位置中。

    打开前面的URL地图服务部分,可以看出是用ArcServer平台发布的服务。我们可以查看这个影像服务的坐标系,瓦片的大小和精度,影像有多少个级别,影像的分辨率等信息。这里的影像分了12个级别,从 Level ID 0 到 Level ID 11,其中每个级别的Start TileEnd Tile存储着起始瓦片和终止瓦片的信息。服务页面的部分信息展示如下:
    server_info
    通过右击页面查看网页源代码,可以看到每个级别影像暴露出来的Start TileEnd Tile的URL地址,这也是我们自动化下载影像的入口。
    leve_info
    通过Start TileEnd Tile显示出来的行列号,我们可以构造出该级别下所有瓦片的行列号,进而构造出该级别所有瓦片的地址。一个地址对应的是一张瓦片PNG,我们只需要全部下载后拼接就能得到一整幅地图。

    比如本例中Level 0级别的 Start Tile 和 End Tile 分别为“93/70”和“95/71”。然后我们需要构造改级别下所有的行列号,如下:
    (93,70),(93,71)
    (94,70),(94,71)
    (95,70),(95,71)
    可以看出在级别0的影像瓦片只有6张,非常的少,但是随着级别的提高,瓦片的数量呈几何上升。通过行列号构造出所有的瓦片URL后就能下载所有的瓦片png了,进而再将所有的png拼接起来。

    2 如何使用FME抓取瓦片图

    我们在前面讲了抓取瓦片图的原理,现在就开始讲如何用FME抓取影像并拼接。因为我做的FME工作流较长,我就挑重点转换器的原理来讲解,最后会提供FME模板供大家下载参考

    在新建的FME工程中,现新建2个公共参数,如下:

    • Tile_URL :存储影像服务地址
    • Level_id:存储要抓取地图的影像级别

    2.1 获取影像页面信息

    我们在得到一幅影像的服务地址后(比如天地图的影像服务),可以使用HttpCaller转换器来获取整个页面的信息,再用StringSearcher转换器获取该页面上12个级别的Start TileEnd Tile信息。

    • HttpCaller:Http请求转换器,通过URL可以发起Http\Https请求来获取站点信息。

    • StringSearcher:字符查找转换器,可以进行简单字符查找,可以使用正则表达式进行高级查找。

    获取了12个级别的Start TileEnd Tile信息后,我们需要指定级别来获取瓦片,否则12个级别的瓦片数量会非常庞大,普通电脑的存储空间根本不够。博主这边测试了Level0到Level5的瓦片下载和拼接,Level8级别的瓦片因数量过多(186110个),FME在进行瓦片镶嵌时崩溃了。

    2.2 生成矩阵数组

    在指定要获取的影像界别后,需要基于Start TileEnd Tile生成所有瓦片的地址信息。在这里博主使用自定义转换器结合Python代码来实现。博主自定义了一个名叫“MatrixArray2”的转换器,接收起始位置和终止位置的XY信息,可以生成步长1的矩阵列表,输出的矩阵列表表示为array_x和array_y,另外通过index_x和index_y保留从0开始计数的矩阵位置信息。
    MatrixArray

    2.3 获取瓦片图

    在获取整个影像或者局部影像的时候,我们首先是需要获取到起始和终止影像的信息,然后来构造区域内所有瓦片图的地址信息。我们在上一个步骤中获取到的矩阵数组,需要进一步构造出每个瓦片图的地址信息。这里我们使用AttributeCreater转换器,然后输入以下表达式:

    $(Tile_URL)/tile/@Value(level_id)/@Value(array_x)/@Value(array_y)
    

    在构造出每个瓦片图的地址信息后,使用ImageFetcher转换器来下载瓦片图。在Image URL参数中输入构造的瓦片图地址,在Image Type参数中输入要下载图片类别。图片类型需要根据影像服务对应的图片类型来设置,否则会下载不了,大部分是使用PNG格式,也存在PNG和JPG混合使用的模式,该模式需要分开下载。
    ImagFetcher

    2.4 瓦片图定位

    下载好的每个瓦片图是没有空间信息的,需要我们通过需要来为每个瓦片图设置相对的空间位置。要为瓦片图设置空间位置,有2个转换器可以实现,分别是Offsetter和RasterGeoreferencer转换器。

    • Offsetter:设置空间偏移
    • RasterGeoreferencer:为栅格设置空间位置

    两者设置的空间偏移量都是一致的:
    Offsetter

    2.5保存瓦片并镶嵌瓦片

    保存瓦片直接添加写模块然后设置好变量文件名就可以了,瓦片镶嵌使用RasterMosaicker转换器,然后添加写模块保存到本地。这个步骤简单,就不多介绍了。

    3 工程下载

    FME工程下载地址:getTilePNG.fmw

    4 总结

    影像抓取也属于爬虫的一类,暴力的信息抓取会给服务器造成巨大的负担,影像或者坐标信息的泄露也有可能会承担相应的责任,希望大家能够合理依法的利用技术,维护国家安全,拥有保密意识。

    展开全文
  • 文章目录背景 背景 我们知道,地图分为栅格和矢量两种。...现在网络上看到的百度、高德、腾讯地图等,基本都是基于矢量切片来显示的,而遥感卫星影像,自然还是栅格切片。 百度地图个性在线编辑器旧版 ...
  • 在利用高分影像半自动化构建单木树种遥感影像样本集过程中,采用了影像冠层切片(CSI)圈定、人工标注、数据增强等方法;同时为了训练单木树种遥感影像样本集,对5个CNN模型进行针对性改写。通过对比分析发现:LeNet5_relu...
  • 本博文介绍的脚本,能够较为方便在指定区域批量地将遥感影像裁剪成固定大小的切片。 2. 样本准备 影像以及对应的点矢量 3. 基于gdal的裁剪代码 # -*- coding: utf-8 -*- from osgeo import ogr import os, sys ...
  • 针对近海、内河场景中船只检测准确性低的问题,提出了一种基于短波红外遥感影像实现水体分割和船只自动检测的方法。利用水体在短波红外波段反射率低的特点,采取阈值分割和形态学处理的方法,从影像中快速准确地提取...
  • 解决影像切片化的快速服务发布,支持国内外不同类型的影像,可用于支撑影像生产和质检、成果影像管理和发布、历史影像归档和管理等应用。 平台特点 • 支持TIFF和IMG格式原始影像发布,坐标系为WGS84或CGCS2000...
  • 本文主要介绍基于C#+GDAL-Python实用工具开发的遥感影像批量处理工具,该工具目前主要包括影像批量切片生成KML文件和影像批量切片生成Tiles文件。该工具.Net框架版本为4.0,GDAL版本为1.11_32位版,Python版本为...
  • python 遥感图像分类

    2021-06-04 15:36:54
    本文采用阈值切片法对遥感影像进行分类,使用数据为landsat8遥感影像,下载地址: http://github.com/GeospatialPython/Learn/blob/master/thermal.zip?raw=true 用envida'k
  •  最近做遥感影像切片,用到了pillow包来处理图片,今天就简单梳理一下相关的内容,分享给大家也备日后查阅。  导入相关包,打开影像并打印相关信息,查看图片: from PIL import Image from PIL import ...
  • arcmap影像变色问题记录

    千次阅读 2019-06-29 19:19:27
    之前用GDAL做了遥感影像切片,但是在ArcMap中显示的时候会变色,刚开始以为是保存的影像的位数问题,即影像是UInt16,然后java中没有UInt类型。最后发现是因为arcmap应该会对输入的影像进行拉伸处理,也就是说在...
  • 影像黑边裁剪的方法

    2020-03-24 14:13:00
    使用影像下载工具下载的影像一般都会有黑边的情况出现,当我们进行切片操作之前,要将其黑边处理为透明色再进行切片,否则加载切片后的影像在场景中也会显示黑边,严重影响视觉效果,基于此梳理遥感影像黑边处理流程...
  • 遥感技术的飞速发展为人们提供了大量的地球空间信息,然而制作地图缓存即切图过程本身是一项非常耗时的工作,按传统的切缓存做法,存在着影像数据不能及时上线使用,即时效性差;切片制作过程频繁、难以管理,需要...
  • 当前主流GIS软件以及...实验采用不同级别大小的遥感影像进行测试,结果表明ParaTile在面对不同规模的数据时,无论从速度还是算法稳定性上都较现有算法和工具具有显著优势,特别是当数据量越大时,这种优势愈加明显。
  • 三种方法我都试过,踩了不少坑,最后是通过第二种方法完成遥感影像的分割。 第一种-栅格分割(10.5的arcgis失败了) 分割栅格工具在Data Management Tools-栅格-栅格处理下,分割栅格有三种类型,分别是 SIZE_...
  • WebGIS开创影像应用新模式,如今,我们不再仅仅借助桌面端实现遥感影像应用,通过构建在线应用模式可以充分利用数据中心的软硬件资源及带宽,同时可以实现多人共享协同的交互操作,高效的挖掘影像价值。 目前,我们...
  • 1、遥感影像发布为切片服务,这是最常规的一种地图服务发布方式; 2、移动端使用的tpk地图包发布,通常也可以通过第一种导出为tpk地图包; 3、数据量较小的shp要素文件发布为动态服务; 4、发布3维场景(dem+遥感...
  • GIS与大数据

    2016-12-13 10:12:35
    说明:文章转自网上多篇文章的部分内容,自己网上看...GIS界最头疼的地图缓存切片生成和存储的问题、海量遥感影像的计算、大规模批空间插值、海量空间数据聚合、空间数据处理等等都有望通过大数据得到一定程度的解决。
  • 遥感影像tif数据通过Geoserver进行发布,为了加快服务器响应速度,我们使用Geoserver对遥感数据做了缓存处理,如下图,使用离线瓦片,我这里出现了中文不识别的问题,因此没法直接Submit去执行切片任务,F12查看...
  • Cesium - 离线使用方法

    千次阅读 2018-04-27 20:31:11
    http://www.cnblogs.com/defineconst/p/6107609.html使用Cesium可以直观的看基于DEM切片产生的Terrain地形数据,有种身临其境的感觉,但缺点是Cesium默认缺省加载了微软Bing提供的地形以及遥感影像数据,可以跟踪...
  • geotools原生java工具包,适用于空间数据的解析和转换,比如点线面矢量数据的增删改查,坐标转换等操作,遥感影像的镶嵌、裁剪等操作。 官网文档地址: http://docs.geotools.org/latest/userguide/ 2.2.Image ...
  • GIS中JPG、PNG、TIF等几种图片格式的区别JPEGPNGTIFMrSID在切片服务中几种图片格式的取舍PNGPNG8PNG24PNG32JPEGMIXEDLERC ...不要使用JPEG或者其他任何有损的压缩格式去对遥感影像进行压缩,分辨率

空空如也

空空如也

1 2
收藏数 37
精华内容 14
关键字:

遥感影像切片