2

I'm getting some strange results in a multivalues to points extraction from a NETCDF4 layer from the MERRA reanalysis product. In particular, it seems like the values being assigned are way off in some areas, as you can see from the different colours of the points and the raster, which should be the same since they share the same scale: Extracted values and their corresponding raster values It seems like its not a reprojection problem, since the points seem to be where they should be, judging from the contour map of Peru also included.

The code I'm using to extract the multivalues to points (and I'm looping over a long bunch of rasters) is the following:

import os, sys
import Tkinter, tkFileDialog
import glob
import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Ask user to select a path to Working directory (i.e. where the raw .nc4 files are)
root = Tkinter.Tk()
root.withdraw()
working_dir =    tkFileDialog.askdirectory(parent=root,initialdir="/",title='Please select a directory')

# Ask user to select CSV file with ENAHO X/Y points
root = Tkinter.Tk()
root.withdraw()
CSV_file = tkFileDialog.askopenfilename(parent=root,initialdir="/",title='Please select CSV file')

# Set the arcgis workspace
arcpy.env.workspace = working_dir
print "Your working directory is:"
print working_dir

####################################################
# Make a XY layer from the CSV
try:
        # Set variables
        x_coords = "longitud"
        y_coords = "latitud"
        out_layer = "enahoXY"
        # 4326 is the code for GCS_WGS_1984
        sr = arcpy.SpatialReference(4326)

        # Make the XY layer
        arcpy.MakeXYEventLayer_management(CSV_file, x_coords, y_coords, out_layer, "", "")

except:
        # If an error occurs print message to screen
        print arcpy.GetMessages()

# Save as .shp
shapefile = os.path.join(working_dir,'/results.shp')
arcpy.CopyFeatures_management(out_layer, shapefile)

####################################################

# Loop around files in directory, but only if they have a .nc extension
files_in_dir = os.listdir(working_dir)
print "Current file processing:"
for f in files_in_dir:
        if f.endswith('.nc4'):

                # Get a part of file that is unique (eg. A2007033)
                parts = f.split('.',2)
                header = parts[2]
                data1 = "tmax"
                extension = ".tif"

                # Program extraction
                temp = "T2MMAX" # maximum daily temperature at 2m (in K)
                x_dimension = "lon"
                y_dimension = "lat"
                band_dimension = "#"
                dimension = "#"
                valueSelectionMethod = "BY_VALUE"

                # Extract temperature
                subset = header + data1 + extension
                print subset
                arcpy.MakeNetCDFRasterLayer_md(f, temp, x_dimension, y_dimension, subset, band_dimension, dimension, valueSelectionMethod)
                # Extract values from Rasters - the column takes the name from the subset
                ExtractMultiValuesToPoints(shapefile, subset, "NONE")

Any ideas what could be the problem?

I'm pretty sure I've entered the correct CRS for the shapefile in the code, but ArcMap tells me the "results.shp" file has no geographical references, which seems strange.

I'm not using the h5py nor the netCDF4 packages for python in the code above.

closed as off-topic by PolyGeo Sep 10 '16 at 5:01

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • I think you have presented too much code which is probably unrelated to the question that you are asking and that you should review the advice on writing code snippets to try and create one which can be used to make your question more answerable. – PolyGeo Sep 10 '16 at 5:03