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

〇、目錄

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

一、前言

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

---

二、基本了解

• 什么是NC格式

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

  ---

• 什么是GeoTiff格式

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

  ---

• Python處理nc文件需要用到的庫

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

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

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

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

三、功能實現(xiàn)

查看NC文件信息

import netCDF4 as nc import numpy as npitem = r'***.nc'#選取一個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”(經度)的基本信息 print('{:-<100}'.format('----'))TEMP = NC_DS.variables['temp'][666,123,:] # 這是一個三維數(shù)據(jù)的NC文件，["緯度","經度","溫度"]，需要按實際修改 print(TEMP) # 查看特定維度數(shù)據(jù)

---

NC文件轉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],即生成 某天 全緯度、經度范圍 圖像PREC = np.asarray(PREC) #數(shù)據(jù)類型轉換PREC[np.where(PREC == -33321)] = np.nan #異常值處理#影像的部分信息LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] #構建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)#設置影像的顯示范圍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)#這里通常選取的是左上角點的坐標(lonmin,latmax)，但如此我的圖像會上下顛倒，不知為何#圖像分辨率xsize ysize通常在NC屬性中給出，也可以 xsize = (lonmax-lonmin)/len(lon)#設置地理坐標系和投影坐標系，這里設置為WGS_84,索引號4326srs = osr.SpatialReference()# 獲取地理坐標系統(tǒng)信息，用于選取需要的地理坐標系統(tǒng)print(type(srs))print(srs)srs.ImportFromEPSG(4326) # 定義輸出的坐標系為"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ù)寫入內存prec_ds.GetRasterBand(1).GetStatistics(0,1)# 添加統(tǒng)計信息prec_ds.FlushCache() # 將數(shù)據(jù)寫入硬盤prec_ds = None # 關閉

---

設置自定義投影

#設置地理坐標系和投影坐標系，這里設置為WGS_84_Alberssr = osr.SpatialReference() #獲取地理坐標系統(tǒng)信息，用于選取需要的地理坐標系統(tǒng)print(type(sr))print(sr)sr.SetProjCS('WGS_84_Albers')sr.SetWellKnownGeogCS('WGS84')#設置地理坐標系為WGS_84sr.SetLCC(35,50,0,90,0,0)#設置等面積圓錐投影的參數(shù)wkt = sr.ExportToWkt()#print(wkt)sr.ImportFromWkt(wkt) # 將定義好的參數(shù)導出為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)#設置nodata值為-9999，ArcGIS自動識別為nodata# 添加統(tǒng)計信息 ds.GetRasterBand(1).GetStatistics(0,1) # 關閉文件ds = None

---

四、成圖預覽

生成的 Tiff 圖像

---

ArcGIS打開屬性

---

掩膜裁剪成圖

---

五、參考

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

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

Python3.GDAL從shp文件生成mask.

---

六、總結

本文包括了NC文件讀取，NC文件轉TIFF、TIFF文件掩膜裁剪、自定投影操作。吸收了多篇文章內容整合而成。代碼多有錯誤，煩請指正。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。