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:
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, sysimport Tkinter, tkFileDialogimport globimport arcpyfrom arcpy import envfrom arcpy.sa import *# Check out the ArcGIS Spatial Analyst extension licensearcpy.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 pointsroot = Tkinter.Tk()root.withdraw()CSV_file = tkFileDialog.askopenfilename(parent=root,initialdir="/",title='Please select CSV file')# Set the arcgis workspacearcpy.env.workspace = working_dirprint "Your working directory is:"print working_dir##################################################### Make a XY layer from the CSVtry: # 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 .shpshapefile = 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 extensionfiles_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, isn't that strange?
I'm not using the h5py nor the netCDF4 packages for python in the code above. Is that a problem?
It's my first post on this forum so I would appreciate any comments regarding the best way to make these questions and or/python code best practice.
أكثر...
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, sysimport Tkinter, tkFileDialogimport globimport arcpyfrom arcpy import envfrom arcpy.sa import *# Check out the ArcGIS Spatial Analyst extension licensearcpy.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 pointsroot = Tkinter.Tk()root.withdraw()CSV_file = tkFileDialog.askopenfilename(parent=root,initialdir="/",title='Please select CSV file')# Set the arcgis workspacearcpy.env.workspace = working_dirprint "Your working directory is:"print working_dir##################################################### Make a XY layer from the CSVtry: # 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 .shpshapefile = 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 extensionfiles_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, isn't that strange?
I'm not using the h5py nor the netCDF4 packages for python in the code above. Is that a problem?
It's my first post on this forum so I would appreciate any comments regarding the best way to make these questions and or/python code best practice.
أكثر...