nc2TFF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "start_time": "2023-08-11T10:20:43.184118Z",
     "end_time": "2023-08-11T10:20:43.184250Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as nc\n",
    "from osgeo import gdal,osr\n",
    "import numpy as np\n",
    "from utils import getfiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "data": {
      "text/plain": "{'lon': <class 'netCDF4._netCDF4.Dimension'>: name = 'lon', size = 576,\n 'lat': <class 'netCDF4._netCDF4.Dimension'>: name = 'lat', size = 361,\n 'time': <class 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'time', size = 24}"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filePath = '/Users/bytedance/Downloads/MERRA2_400.tavgU_2d_slv_Nx.202306.nc4'\n",
    "target = 'TS'\n",
    "outtif = f'./out/{target}.tiff'\n",
    "ncdata = nc.Dataset(filePath, mode='r')\n",
    "ncdata.dimensions"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "start_time": "2023-08-11T10:41:30.259149Z",
     "end_time": "2023-08-11T10:41:30.259436Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/b1/0fd1b6hs7lz0fm_mh346lybm0000gn/T/ipykernel_76654/615470046.py:1: RuntimeWarning: invalid value encountered in cast\n",
      "  time = ncdata.variables['time'][:]\n",
      "/var/folders/b1/0fd1b6hs7lz0fm_mh346lybm0000gn/T/ipykernel_76654/615470046.py:1: UserWarning: WARNING: valid_range not used since it\n",
      "cannot be safely cast to variable data type\n",
      "  time = ncdata.variables['time'][:]\n"
     ]
    },
    {
     "data": {
      "text/plain": "masked_array(data=[0],\n             mask=False,\n       fill_value=999999,\n            dtype=int32)"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "time = ncdata.variables['time'][:]\n",
    "time"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "start_time": "2023-08-11T10:23:55.111497Z",
     "end_time": "2023-08-11T10:23:55.111573Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'TS'",
     "output_type": "error",
     "traceback": [
      "\u001B[0;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[0;31mKeyError\u001B[0m                                  Traceback (most recent call last)",
      "Cell \u001B[0;32mIn[3], line 15\u001B[0m\n\u001B[1;32m     13\u001B[0m GeoTransform \u001B[38;5;241m=\u001B[39m (LonMin,Lon_Res, \u001B[38;5;241m0\u001B[39m, LatMax,\u001B[38;5;241m0\u001B[39m,\u001B[38;5;241m-\u001B[39mLat_Res)\n\u001B[1;32m     14\u001B[0m \u001B[38;5;66;03m# 读取数据，处理异常值\u001B[39;00m\n\u001B[0;32m---> 15\u001B[0m var \u001B[38;5;241m=\u001B[39m \u001B[43mncdata\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mvariables\u001B[49m\u001B[43m[\u001B[49m\u001B[43mtarget\u001B[49m\u001B[43m]\u001B[49m\n\u001B[1;32m     16\u001B[0m var_nans \u001B[38;5;241m=\u001B[39m var[:]\n\u001B[1;32m     17\u001B[0m _FillValue \u001B[38;5;241m=\u001B[39m var\u001B[38;5;241m.\u001B[39m_FillValue\n",
      "\u001B[0;31mKeyError\u001B[0m: 'TS'"
     ]
    }
   ],
   "source": [
    "# longitude and latitude\n",
    "lon = ncdata.variables['lon'][:]\n",
    "lat = ncdata.variables['lat'][:]\n",
    "# time = data.variables['time'][:]\n",
    "#获取四个脚点坐标\n",
    "LonMin,LatMax,LonMax,LatMin = [np.min(lon),np.max(lat),np.max(lon),np.min(lat)]\n",
    "# 设置图像分辨率\n",
    "row = len(lat)\n",
    "col = len(lon)\n",
    "Lon_Res=(LonMax - LonMin) / (float(col) - 1)\n",
    "Lat_Res=(LatMax- LatMin) / (float(row) - 1)\n",
    "\n",
    "GeoTransform = (LonMin,Lon_Res, 0, LatMax,0,-Lat_Res)\n",
    "# 读取数据，处理异常值\n",
    "var = ncdata.variables[target]\n",
    "var_nans = var[:]\n",
    "_FillValue = var._FillValue\n",
    "var_nans[var_nans == _FillValue] = np.nan\n",
    "tmp_data = var_nans[0,:,:]\n",
    "# 灰度值转换\n",
    "# tmp_data = 256 * (tmp_data - np.min(var))/(np.max(var)-np.min(var))\n",
    "tmp_data = tmp_data[::-1]#这里是需要倒置一下的 这个很重要！！！！！\n",
    "\n",
    "# 判断栅格数据的数据类型\n",
    "if 'int8' in tmp_data.dtype.name:\n",
    "    datatype = gdal.GDT_Byte\n",
    "elif 'int16' in tmp_data.dtype.name:\n",
    "    datatype = gdal.GDT_UInt16\n",
    "else:\n",
    "    datatype = gdal.GDT_Float32\n",
    "\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "outRaster = driver.Create(outtif, col, row, 1, datatype)\n",
    "\n",
    "outRaster.SetGeoTransform(tuple(GeoTransform))\n",
    "# 获取数据集第一个波段，是从1开始，不是从0开始\n",
    "outband = outRaster.GetRasterBand(1)\n",
    "outband.WriteArray(ncdata)\n",
    "\n",
    "# 数据投影信息\n",
    "outRasterSRS = osr.SpatialReference()\n",
    "# 代码4326表示WGS84坐标\n",
    "outRasterSRS.ImportFromEPSG(4326)\n",
    "outRaster.SetProjection(outRasterSRS.ExportToWkt())\n",
    "\n",
    "outband.FlushCache()\n",
    "out_ds=outband= None # 关闭spei_ds指针，注意必须关闭"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}