处理VOC数据

一、先提取出来VOC格网数据

1、voc提取代码是将提取出来的数据给保存成xlsx格式。

# -*- coding: utf-8 -*-
import os
from netCDF4 import Dataset
from openpyxl.workbook import Workbook


def sbdx(source_file,type):
    nc_obj1 = Dataset(r"F:\VOC_clear\cmaq.lat_lon") #读取lat_lon文件 （这个格式是nc文件）
    nc_obj2 = Dataset(source_file) #读取conc文件（这个格式也是nc文件）
    lat = (nc_obj1.variables['LAT'][:]) #读取lat_lon文件中的lat
    lon = (nc_obj1.variables['LON'][:]) #读取lat_lon文件中的lon
    day = (nc_obj2.variables['TFLAG'][:]) #读取conc文件中存储日期的属性
    tp = (nc_obj2.variables[type][:]) #读取conc文件中存储变量中的属性

    wb = Workbook() #首先获取excel文件生成对象wb 使用openpyxl.load_workbook方法
    sh = wb.active #获取当前工作簿
    sh.cell(row=1,column=1).value = "lat" #第一行第一列为标题lat
    sh.cell(row=1,column=2).value = "lon" #第一行第二列为标题lon
    year = int(str(source_file[-19:-15])) #从conc文件中识别年份 这个数字要根据给的路径具体修改,小到大

    #从各个conc文件中查看可知不同年份有数据的天份不同
    if year == 2013 or year == 2014 or year == 2015 or year == 2018 or year == 2019:
        d = 363
    elif year == 2016:
        d = 364
    elif year == 2017:
        d = 362
	#经纬度范围根据具体文件修改
    for i in range(127): 
        for j in range(197):
        	#读取127*197行经纬度数据
            sh.cell(row=197 * i + j + 2, column=1).value = lat[0, 0, i, j] 
            sh.cell(row=197 * i + j + 2, column=2).value = lon[0, 0, i, j]
            #循环时间
            for k in range(d):
                sh.cell(row=1,column=3+k).value = day[k,0,0] #给第一行写入时间
                sh.cell(row = 197 * i + j + 2,column=3+k).value = tp[k,0,i,j]
               
    #保存文件
    wb.save("F:/VOC_clear/" + str(type) + "_" + str(year) + "_.xlsx")


if "__name__ = __main__":

    path = "F:/VOC_clear/yuanshiData/"
    path_list = os.listdir(path)
    path_list.sort(key=lambda x:str(x[:-5]))
    for filename in path_list:
        filename = path+filename  #上面获取的只是文件名字，要传递路径，还要重新加一下路径
        #成分变量
        type1 = ["ALK", "HCHO", "ACIDS", "VOC", "ACETONE", "ACETYLENE", "ACROLEIN", "BACL", "BUTADIENE13", "CCHO",
                 "CRES", "ETHENE", "ETOH", "GLY", "MEK", "MEOH", "MGLY", "MVK", "OLE1", "OLE2", "PRD2", "PROPENE",
                 "RCHO", "RNO3", "AORGAJ", "AORGBJ"]
        #循环输入成分变量
        for type in type1:
            sbdx(filename,type)
    print("successful!")  

2、然后在用xlsx格式转成csv格式。

import pandas as pd

#修改一下输入和输出路径就可以
data = pd.read_excel(r"输入xlsx文件路径",index_col=0,engine="openpyxl")   #xlrd旧版本不能打开xlsx格式了，使用openyxl替代打开xlsx  index_col=0的作用是去掉第一列索引 0，1，2.。。（不加的话会添加一列索引）,
data.to_csv('输出csv文件路径',encoding='utf-8')

注： 提取VOC代码在师兄指导下完成（我太水了）

二、清理至区县上面

1、首先将voc数据提取出来的经纬度在ArcGIS中展开（该数据时兰伯特投影），做出点shp文件，然后和标准地图做一个空间连接，如下图

2、查找空的区县

   这样操作肯定会有点没有落到标准区县上面，所以此时我们要做的是找出空的区县（如何找看以前的博客文章），然后使用 -要素转点 -功能，找出空的区县的中心点，然后使用-近邻分析-功能（邻近要素选择的是空间连接的索引），找出空的区县中心点最近点的索引坐标，将这个和前面做出来不缺的区县的索引放在一块，然后得到了2901个区县的全部格网文件。

3、使用python的groupby分组功能将数据清理至区县上面

说明： 这里匹配的时候使用的是索引来匹配的（缺的区县的NEAR_FID和空间连接的FID（验证一下经纬度是否一样）），使用经纬度来匹配的话会因为小数点问题，匹配不到数据（困扰了好长时间）

代码如下：

import pandas  as pd
import os

def clear_to_county(voc_lujing):
    county_grid = pd.read_csv(r"F:\VOC_clear\完整的区县和经纬度坐标.csv",encoding="GBK") #区县经纬度格网
    FID_grid = pd.read_csv(r"F:\VOC_clear\FID用于连接voc格网进行匹配合并.csv") #FID格网，用于匹配voc格网
    voc_grid = pd.read_csv(voc_lujing) #voc格网

    VOC_FID_lianjie = pd.concat([FID_grid,voc_grid],axis=1) #用于列合并

    lianjie = pd.merge(county_grid,VOC_FID_lianjie,left_on='FID1',right_on="FID",how='inner')

    #分组代码
    county_clear = lianjie.groupby(['PAC','NAME','省','市代码','市']).mean()
    #删除不要的列
    county_clear.drop(['FID','FID1','lat_x','lon_x','lat_y','lon_y'], axis=1, inplace=True) # inplace=True, 直接从内部删除 del一次只能删除一个，具体的看csdn
    #截取voc变量名字
    var_voc = voc_lujing[15:-4] #具体路径修改数字
    county_clear.to_csv("F:/VOC_clear/" + str(var_voc) + ".csv", encoding="utf-8")  #输出结果

if "__name__ = __main__":
    #批量读取voc格网数据
    path = "F:/VOC_clear/1/"
    path_list = os.listdir(path)
    path_list.sort(key=lambda x:str(x[:-5]))
    for filename in path_list:
        filename = path+filename  #上面获取的只是文件名字，要传递路径，还要重新加一下路径
        clear_to_county(filename)

感想： 一开始用R语言来做的，结果经纬度问题阻挡了我很长时间，后来发现，R中的summarize（）函数还不能把365天的数据各自汇总求均，8号这天突然想到了python，然后就用python来匹配经纬度，结果还是小数点问题，后来突然想到了索引，然后翻阅之前做的格网文件，发现空区县有索引，正好VOC数据和标准地图空间连接的索引是相对应的（手动查找对照），然后就这样出来了啊，期间也出现了好多问题（如python如何合并列），度娘了多少个博客才找到解决办法，总之，真的是一份耕耘一分收获，收获很大。
PS： 代码要是很强势的话也不会这么难受了，难搞啊