MODIS Tile 拼接和重投影

柚柚柚子爱吃葡萄柚

已于 2024-06-14 23:08:40 修改

阅读量789

点赞数

分类专栏： 2. 数据分析（代码和软件）文章标签： python

于 2022-04-08 17:39:46 首次发布

本文链接：https://blog.csdn.net/Wangxr_cau/article/details/124046948

版权

2. 数据分析（代码和软件）专栏收录该内容

8 篇文章 0 订阅

订阅专栏

import arcpy, os, sys
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
from os import listdir
from os.path import isfile, join
import numpy as np
from osgeo import gdal

##Mosaic To New Raster
##Usage: MosaicToNewRaster_management inputs;inputs... output_location raster_dataset_name_with_extension 
##                                    {coordinate_system_for_the_raster} 8_BIT_UNSIGNED | 1_BIT | 2_BIT | 4_BIT 
##                                    | 8_BIT_SIGNED | 16_BIT_UNSIGNED | 16_BIT_SIGNED | 32_BIT_FLOAT | 32_BIT_UNSIGNED 
##                                    | 32_BIT_SIGNED | | 64_BIT {cellsize} number_of_bands {LAST | FIRST | BLEND  | MEAN 
##                                    | MINIMUM | MAXIMUM} {FIRST | REJECT | LAST | MATCH}    


# ----------------------
# STEP1: Read data 
# ----------------------

# Date 
year         = 2008;  
date_list = np.arange(1,362,8);

# Output folder
folder_ouput = "E:/MODIS/GPP500m/gpp_"+str(year)+"/gpp_"+str(year)+"_mosaic";


# Input file folder
folder_input_list = ['E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h08v04',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h08v05',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h08v06',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h09v04',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h09v05',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h09v06',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h10v04',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h10v05',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h10v06',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h11v04',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h11v05',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h12v04',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h12v05',\
                     'E:/MODIS/GPP500m/gpp_'+str(year)+"/gpp_"+str(year)+'/h13v04',\
                     ]
    
# Tiles which we want to mosaic
tile_list = ['.h08v04.tif',
             '.h08v05.tif',
             '.h08v06.tif',
             '.h09v04.tif',
             '.h09v05.tif',
             '.h09v06.tif',  
             '.h10v04.tif', 
             '.h10v05.tif', 
             '.h10v06.tif', 
             '.h11v04.tif', 
             '.h11v05.tif', 
             '.h12v04.tif', 
             '.h12v05.tif',
             '.h13v04.tif' ]
             
           
# ----------------------
# STEP2: Processing
# ----------------------
for date in date_list:
    
    #### Defeine output file name ####
    file_list = []
    filename1  = 'GPP.'+str(year)+'{:03d}'.format(date)+'.tif'
    filename2  = 'GPP.'+str(year)+'{:03d}'.format(date)+'_wgs84_0005D'+'.tif'
    
    #### Read input file ####
    for tile, folder in zip(tile_list,folder_input_list):
    
       file_list.append( folder+ '/'+'GPP.'+str(year)+'{:03d}'.format(date)+ tile)
    
    
    #### Merge all raster ####
    mosaic = arcpy.MosaicToNewRaster_management(file_list, folder_ouput, filename1, "",  "16_BIT_SIGNED", "", "1", "LAST","FIRST")
    print(filename1)
    


    #### eproject to wgs84 ####
    
    # open dataset
    #dataset = gdal.Open(filename1,gdal.GA_ReadOnly)
    
    # gdalwarp
    #kwargs = {'format': 'GTiff', 'dstSRS': 'EPSG:4326','xRes':0.005,'yRes':0.005}
    #ds = gdal.Warp(destNameOrDestDS=filename2,srcDSOrSrcDSTab=dataset, **kwargs)
    #del ds 



[1] MosaicToNewRaster_management 
主要参考：https://community.esri.com/t5/python-questions/merge-rasters-mosaic-to-new-raster/td-p/88308
官方文档：https://pro.arcgis.com/zh-cn/pro-app/2.7/tool-reference/data-management/mosaic-to-new-raster.htm
CSDN: https://www.twblogs.net/a/5d539b54bd9eee541c317cd9/?lang=zh-cn

[2] Reproject to wgs84
reproject to wgs84:https://gis.stackexchange.com/questions/386381/converting-
reproject to wgs84: https://www.programcreek.com/python/example/116446/gdal.Warp (!)