Select to view content in your preferred language

NetCDF to csv script

7472
1
09-24-2012 01:35 PM
ChristopherStrother
Deactivated User
I have 27 .nc files and want to extract the precipitation "prcp" values above 25mm and write them along with their x,y coordinates to a .csv file for further analysis in ArcGIS. I have 10.1 installed as well as Python 2.7 with numpy and scipy modules. Any ideas on why this code isn't working?

Error:
Traceback (most recent call last):
  File "E:\UGA\COWEETA\Daymet\netcdf2csv.py", line 69, in <module>
    lat = var.getValue()
  File "C:\Python27\ArcGIS10.1\lib\site-packages\scipy\io\netcdf.py", line 784, in getValue
    return self.data.item()
ValueError: can only convert an array  of size 1 to a Python scalar

import sys
import os
import csv
import numpy
import scipy.io
from scipy.io import netcdf

path_netcdf='C:\\PRCP11028\\'

THRESHOLD = 25

x_dimension = "x"
y_dimension = "y"
band_dimension = ""
dimension = "time"
dimension_values = ""
valueSelectionMethod = "BY_VALUE"

extension = ".nc"

files=os.listdir(path_netcdf)

fcsv = open(path_netcdf+"extremes.txt","wb")
csvFile = csv.writer(fcsv)

for fname in files:
    
    if fname.lower().endswith(extension):

        print fname
        
        variable, extension = fname.split(".")

        netcdfFile = netcdf.netcdf_file(os.path.join(path_netcdf,fname), 'r')

        allDimNames = netcdfFile.dimensions.keys()

        for dimName in allDimNames:
            dimValue = netcdfFile.dimensions[dimName]
##            print dimName, ' : ', dimValue

            tileID = getattr(netcdfFile, 'tileid')
            strTile = str(tileID)
            #print tileID

            year = getattr(netcdfFile, 'start_year')
            strYear = str(year)
            #print strYear

        
        variableNames = netcdfFile.variables.keys()
            #print variableNames

        var = netcdfFile.variables["lat"]
            #print var
        lat = var.getValue()
            #print lat
                 
        var = netcdfFile.variables["lon"]   
        lon = var.getValue() 
##        print lon

        var = netcdfFile.variables["x"]   
        x = var.getValue() 
##        print x

        var = netcdfFile.variables["y"]  
        y = var.getValue() 
##        print y

        var = netcdfFile.variables["prcp"]
        prcp = var.getValue() 

##        myshape = var.shape
##        print myshape
       
        for day in range(0,364):
            strDay = str(day+1)
            aboveThresIndices = prcp[day] >= THRESHOLD

            lats = lat[aboveThresIndices]
            lons = lon[aboveThresIndices]
            prcps = prcp[day][aboveThresIndices]

            for iLat in range(0,len(lats)):
##                print lats[iLat], " ", lons[iLat]
                csvFile.writerow([tileID,strDay,lats[iLat],lons[iLat],prcps[iLat]]) 
                

        netcdfFile.close()

fcsv.close()
        
print "-------- finished"

sys.exit() 
Tags (2)
0 Kudos
1 Reply
PhilMorefield
Frequent Contributor
The data you are trying extract by using getValue() is an array-like object (i.e., it's two dimensional). As you can see here, the getValue() method only works on scalars (e.g., an integer, a float, etc.). Consider doing this instead:
# retrieve first time step
var[0]

# retrieve second time step
var[1]

Is your goal to convert these data to rasters? If so, I would skip writing them to csv and use the NumpyArrayToRaster() function.

If what you want is a just a list of the points and values written to a text file, you'll need to look at reshaping the array-like object your data are currently in.
0 Kudos