风云2号云图Mosaic Dataset处理

# ---------------------------------------------------------------------------
# -*- coding: utf-8 -*-
# QiXiangAnalyze.py
# Created on: 2012-07-04 08:15:21.00000
# Creater: GISPathfinder
# Usage: isolineAnalyze <station> <reclassificationString> <isoResult>
# Description:
# ---------------------------------------------------------------------------

import time
import datetime
from struct import *
from bz2 import BZ2File
import sys, os, arcpy

arcpy.env.overwriteOutput = True

start = time.clock()

file = open(r"F:/geodata/shanxi/MapData/GP/QiXiang/FL2/1607080900.000", "r")

#
zonName,flag,dataName = unpack("8s3s40s", file.read(51))
zonName = zonName.decode("gbk").rstrip('x00')
dataName = dataName.decode("gbk").rstrip('x00')

#
print(zonName)
print(dataName)
#
#file.read(5+3+3+3)
year,month,day,hour = unpack("5s3s3s3s", file.read(14)) # 时间说明
#print(year + " " + month+ " " + day+ " " + hour)
dateTimeData = datetime.datetime(int(year),int(month), int(day), int(hour), 0, 0, 0)
fileNameParseResult = dateTimeData.strftime("%Y%m%d%H%M%S")
asciiGrid = "F:\geodata\shanxi\MapData\GP\QiXiang\temp\"+fileNameParseResult+".txt"
print(fileNameParseResult)

XNumGrids,YNumGrids,= unpack("5s5s", file.read(10))
XNumGrids = int(XNumGrids.decode("gbk").rstrip('x00'))
YNumGrids = int(YNumGrids.decode("gbk").rstrip('x00'))
print("X: " + str(XNumGrids)+" Y: "+str(YNumGrids))

StartLon,StartLat,= unpack("8s8s", file.read(16)) # 开始经度 开始纬度
StartLon = float(StartLon.decode("gbk").rstrip('x00'))
StartLat = float(StartLat.decode("gbk").rstrip('x00'))
print("StartLon "+str(StartLon)+" StartLat"+str(StartLat))

project,reso,imageType,= unpack("2s5s2s", file.read(9))
project = int(project.decode("gbk").rstrip('x00')) # 投影方式
reso = float(reso.decode("gbk").rstrip('x00')) # 分辨率
imageType = int(imageType.decode("gbk").rstrip('x00')) # 图象类型
print("project " + str(project))
print("reso " + str(reso))
print("imageType " + str(imageType))

refTable, = unpack("12s", file.read(12))
refTable = refTable.decode("gbk").rstrip('x00')

CenterLon,CenterLat,= unpack("8s8s", file.read(16)) # 中心经度 中心纬度
CenterLon = float(CenterLon.decode("gbk").rstrip('x00'))
CenterLat = float(CenterLat.decode("gbk").rstrip('x00'))
print("CenterLon "+str(CenterLon)+" CenterLat "+str(CenterLat))

textYX = []
for j in range(0, YNumGrids):
textX = []
for k in range(0, XNumGrids):
value, = unpack("B", file.read(1))
textX.append(str(value))
textYX.append(' '.join(textX)+' ')
file.close()

#
textYX.reverse()
#
#-------------------------------------------------------------------------------
#file_object = open('ASCIIDataQPF.txt', 'w')
#
file_object = open(asciiGrid, 'w')
file_object.write("NCOLS " + str(XNumGrids) + " ")
file_object.write("NROWS " + str(YNumGrids) + " ")
file_object.write("XLLCENTER " + str(StartLon) + " ")
file_object.write("YLLCENTER " + str(StartLat) + " ") # round(YReso, 3) *
file_object.write("CELLSIZE " + str(reso) + " ")
file_object.write("NODATA_VALUE " + str(-9999) + " ")
#
#
file_object.writelines(textYX)
file_object.close()
end = time.clock()
dateSpanTransfer = end - start
print("read: %f s" % dateSpanTransfer)
dateSpanTransfer = end - start

#file.read(1)
#-------------------------------------------------------------------------------

projXa80Proj = 'PROJCS["Xian80_Projected_Coordinate_System",GEOGCS["GCS_Xian_1980",DATUM["D_Xian_1980",SPHEROID["Xian_1980",6378140.0,298.257]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",110.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",0.0],UNIT["Meters",1.0]]'

projXa80 = 'GEOGCS["GCS_Xian_1980",DATUM["D_Xian_1980",SPHEROID["Xian_1980",6378140.0,298.257]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'

projWGS84 = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"

rasterData = "F:\geodata\shanxi\MapData\GP\QiXiang\temp\" + fileNameParseResult+".img"

#rasterProject = "F:\geodata\shanxi\MapData\GP\QiXiang\Result\FL2Raster.gdb\FY" + fileNameParseResult
rasterProject = "F:\geodata\shanxi\MapData\GP\QiXiang\Result\" + fileNameParseResult+".img"

inputLinkFile = "F:\geodata\shanxi\MapData\GP\QiXiang\Result\DBZTic.txt"

FL2Dataset = "F:\geodata\shanxi\MapData\GP\QiXiang\Result\FL2Raster.gdb\FL2Dataset"

colorMapFile = "F:\geodata\shanxi\MapData\GP\QiXiang\V-02.clr"

# Process: ASCII to Raster
arcpy.ASCIIToRaster_conversion(asciiGrid, rasterData, "INTEGER")

# Process: Define Projection
arcpy.DefineProjection_management(rasterData, projXa80Proj)

# Process: Register Raster
arcpy.RegisterRaster_management(rasterData, "REGISTER", "", inputLinkFile, "POLYORDER1", "")

# Process: Project Raster
arcpy.ProjectRaster_management(rasterData, rasterProject, projXa80, "NEAREST", 0.15, "", "", projXa80Proj)

# Process: Define Projection
arcpy.DefineProjection_management(rasterProject, projWGS84)

# Process: Add Colormap
arcpy.AddColormap_management(rasterProject, "", colorMapFile)

# Process: Add Rasters To Mosaic Dataset
arcpy.AddRastersToMosaicDataset_management(FL2Dataset, "Raster Dataset", rasterProject, "UPDATE_CELL_SIZES", "UPDATE_BOUNDARY", "NO_OVERVIEWS", "", "0", "1500", "", "", "SUBFOLDERS", "ALLOW_DUPLICATES", "BUILD_PYRAMIDS", "CALCULATE_STATISTICS", "NO_THUMBNAILS", "", "NO_FORCE_SPATIAL_REFERENCE")

#arcpy.AddRastersToMosaicDataset_management(FL2Dataset, "Raster Dataset", rasterProject, "NO_CELL_SIZES", "NO_BOUNDARY", "NO_OVERVIEWS", "", "0", "1500", "", "", "SUBFOLDERS", "ALLOW_DUPLICATES", "BUILD_PYRAMIDS", "CALCULATE_STATISTICS", "NO_THUMBNAILS", "", "NO_FORCE_SPATIAL_REFERENCE")

print("OK")