通過Python實(shí)現(xiàn)NC文件轉(zhuǎn)GeoTiff格式

〇、目錄

• 通過Python實(shí)現(xiàn)NC文件轉(zhuǎn)GeoTiff格式
• • 一、前言
  • 二、基本了解
  • 三、功能實(shí)現(xiàn)
  • 四、成圖預(yù)覽
  • 五、參考
  • 六、總結(jié)

一、前言

• 基于Python和GDAL庫(kù)介紹處理nc文件的基本操作
• 簡(jiǎn)單處理nc文件為GeoTiff文件的Python函數(shù)

---

二、基本了解

• 什么是NC格式

  NC文件即NetCDF (Network Common Data Form），是一種通用的數(shù)據(jù)存儲(chǔ)格式。廣泛用于存儲(chǔ)地學(xué)、大氣科學(xué)、海洋科學(xué)等一系列多維數(shù)據(jù)，可封裝時(shí)間、經(jīng)度、緯度、降水、溫度等多個(gè)維度數(shù)據(jù)，數(shù)據(jù)結(jié)構(gòu)清晰易讀，可以很方便的提取、應(yīng)用。

  ---

• 什么是GeoTiff格式

  GeoTIFF 是一種地理柵格文件，可以理解為專門存儲(chǔ)地理信息的Tiff圖像文件格式。地理坐標(biāo)系、地圖投影等要素直接嵌入到圖像文件中。

  ---

• Python處理nc文件需要用到的庫(kù)

  • NetCDF4用于讀取NC文件；numpy用于數(shù)據(jù)修改；osgeo的子庫(kù)gdal和osr，前者用于實(shí)現(xiàn)轉(zhuǎn)換和柵格編輯功能，后者可用于確定地理坐標(biāo)系和地圖投影。

    import numpy as np import netCDF4 as nc from osgeo import gdal,osr

  • 建議使用 Anaconda 配置虛擬Python環(huán)境，可以很方便地進(jìn)行包管理，我這里使用的環(huán)境版本是

    python 3.6.12(64bit) numpy 1.19.1 netCDF4 1.4.2 gdal 3.0.2

三、功能實(shí)現(xiàn)

查看NC文件信息

import netCDF4 as nc import numpy as npitem = r'***.nc'#選取一個(gè)nc文件路徑 DS = nc.Dataset(item)print(DS,type(NC_DS)) # 返回DS的數(shù)據(jù)類型，<class 'netCDF4._netCDF4.Dataset'> print('{:-<100}'.format('----'))print(DS.variables) # 了解變量的基本信息 print('{:-<100}'.format('----'))print(DS.variables['lon']) # 了解某一變量如“l(fā)on”(經(jīng)度)的基本信息 print('{:-<100}'.format('----'))TEMP = NC_DS.variables['temp'][666,123,:] # 這是一個(gè)三維數(shù)據(jù)的NC文件，["緯度","經(jīng)度","溫度"]，需要按實(shí)際修改 print(TEMP) # 查看特定維度數(shù)據(jù)

---

NC文件轉(zhuǎn)GeoTiff

#NetCDF文件處理 def Convert(item,day): DS = nc.Dataset(item)Lat = DS.variables['lat'][:]Lon = DS.variables['lon'][:]PREC = DS.variables['prec'][day,:,:] #[day,lon,lat],即生成 某天 全緯度、經(jīng)度范圍 圖像PREC = np.asarray(PREC) #數(shù)據(jù)類型轉(zhuǎn)換PREC[np.where(PREC == -33321)] = np.nan #異常值處理#影像的部分信息LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] #構(gòu)建tiff框架prec_ds = gdal.GetDriverByName('Gtiff').Create(OutTif,len(Lon),len(Lat),1,gdal.GDT_Float32) #dst_ds = driver.Create(dst_filename, xsize, ysize,bands, eType=gdal.GDT_Byte)#設(shè)置影像的顯示范圍prec_ds.SetGeoTransform(LonMin,0.1, 0, LatMin, 0, 0.1)#SetGeoTransform(double xLeft,double dX,double yProjOnX,double yTop,double x,ProjOnY,double dY)#這里通常選取的是左上角點(diǎn)的坐標(biāo)(lonmin,latmax)，但如此我的圖像會(huì)上下顛倒，不知為何#圖像分辨率xsize ysize通常在NC屬性中給出，也可以 xsize = (lonmax-lonmin)/len(lon)#設(shè)置地理坐標(biāo)系和投影坐標(biāo)系，這里設(shè)置為WGS_84,索引號(hào)4326srs = osr.SpatialReference()# 獲取地理坐標(biāo)系統(tǒng)信息，用于選取需要的地理坐標(biāo)系統(tǒng)print(type(srs))print(srs)srs.ImportFromEPSG(4326) # 定義輸出的坐標(biāo)系為"WGS 84"，AUTHORITY["EPSG","4326"]prec_ds.SetProjection(srs.ExportToWkt()) # 給新建圖層賦予投影信息# 數(shù)據(jù)寫出prec_ds.GetRasterBand(1).SetNoDataValue(-9999)prec_ds.GetRasterBand(1).WriteArray(PREC) # 將數(shù)據(jù)寫入內(nèi)存prec_ds.GetRasterBand(1).GetStatistics(0,1)# 添加統(tǒng)計(jì)信息prec_ds.FlushCache() # 將數(shù)據(jù)寫入硬盤prec_ds = None # 關(guān)閉

---

設(shè)置自定義投影

#設(shè)置地理坐標(biāo)系和投影坐標(biāo)系，這里設(shè)置為WGS_84_Alberssr = osr.SpatialReference() #獲取地理坐標(biāo)系統(tǒng)信息，用于選取需要的地理坐標(biāo)系統(tǒng)print(type(sr))print(sr)sr.SetProjCS('WGS_84_Albers')sr.SetWellKnownGeogCS('WGS84')#設(shè)置地理坐標(biāo)系為WGS_84sr.SetLCC(35,50,0,90,0,0)#設(shè)置等面積圓錐投影的參數(shù)wkt = sr.ExportToWkt()#print(wkt)sr.ImportFromWkt(wkt) # 將定義好的參數(shù)導(dǎo)出為WKT prec_ds.SetProjection(sr.ExportToWkt()) # 給新建圖層賦予投影信息

---

用掩膜裁剪

from osgeo import gdal import os import timedef clip(input_shape,input_raster,output_raster):# 柵格文件路徑，打開柵格文件input_raster=gdal.Open(input_raster)# 裁剪，gdal.warp函數(shù)，input_shape是掩膜文件ds = gdal.Warp(output_raster,input_raster,format = 'GTiff',cutlineDSName = input_shape, dstNodata = -9999)#設(shè)置nodata值為-9999，ArcGIS自動(dòng)識(shí)別為nodata# 添加統(tǒng)計(jì)信息 ds.GetRasterBand(1).GetStatistics(0,1) # 關(guān)閉文件ds = None

---

四、成圖預(yù)覽

生成的 Tiff 圖像

---

ArcGIS打開屬性

---

掩膜裁剪成圖

---

五、參考

基于Python+GDAL實(shí)現(xiàn)nc格式轉(zhuǎn)geotiff格式.

使用python的netCDF4庫(kù)讀取.nc文件 和 創(chuàng)建.nc文件.

Python3.GDAL從shp文件生成mask.

---

六、總結(jié)

本文包括了NC文件讀取，NC文件轉(zhuǎn)TIFF、TIFF文件掩膜裁剪、自定投影操作。吸收了多篇文章內(nèi)容整合而成。代碼多有錯(cuò)誤，煩請(qǐng)指正。

總結(jié)

以上是生活随笔為你收集整理的通过Python实现NC文件转GeoTiff格式的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。