Error in multiple netCDF file processing followed by raster projection

1644
8
Jump to solution
07-26-2017 08:23 AM
DuminduJayasekera
New Contributor III

Hi,

I have multiple netCDF files I need to make a "netCDF raster" layer and then "project raster" and save as a ".tif" file. I have the following script but this gives me an error message.

Please see the code and error message below:

Thanks in advance.

# Import system modules  
import arcpy, sys 
from arcpy import env  
from arcpy.sa import *  
  
# Input data source  
arcpy.env.workspace = r"C:/Users/jayaskeradl/Desktop/weekly_avg_project"  
arcpy.env.overwriteOutput = True  
  
# Set output folder  
OutputFolder = r"C:/Users/jayaskeradl/Desktop/weekly_avg_project/output"  
  
# Loop through a list of files in the workspace  
NCfiles = arcpy.ListFiles("*.nc")  
  
for filename in NCfiles:  
     print("Processing: " + filename)  
     inNCfiles = arcpy.env.workspace + "//" + filename  
     fileroot = filename[0:(len(filename)-3)]  
     TempLayerFile = "SM_amount"  
     outRaster = OutputFolder + "//" + fileroot  
    # Process: Make NetCDF Raster Layer  
     arcpy.MakeNetCDFRasterLayer_md(filename, "var250", "lon", "lat", TempLayerFile, "", "", "BY_VALUE")  
     arcpy.ProjectRaster_management(filename, TempLayerFile,"GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "NEAREST", "0.125 0.125", "NAD_1983_To_WGS_1984_1", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")
    # Process: Copy Raster  
     arcpy.CopyRaster_management(TempLayerFile, outRaster + ".tif", "", "", "", "NONE", "NONE", "")  
  
  
print "***DONE!!!"  
print arcpy.GetMessages() 


‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍ExecuteError: Failed to execute. Parameters are not valid.
Undefined coordinate system for input dataset.
Failed to execute (ProjectRaster).
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
1 Solution

Accepted Solutions
IanMurray
Frequent Contributor

You have the same input for Make NetCDF Layer and for Project Raster(filename or the .nc file), while you should take the output from Make NetCDF Layer (an actual raster layer, which should have a coordinate system) and use that for the input of Project Raster.  The output of Project Raster would need a different name than you are currently using for the input.

View solution in original post

8 Replies
IanMurray
Frequent Contributor

First off, when posting code to GeoNet please use the code formatting to make your script more legible(See https://community.esri.com/people/curtvprice/blog/2014/09/25/posting-code-blocks-in-the-new-geonet?s...‌ for more on that).

Currently you are getting an indent error, so either something is indented where it shouldn't be, you are mixing spaces and tabs for your indenting.

I do not believe that will be the only issue since you code has some fairly serious problems.  since the foreslash character is an escape character, the file paths you used in your workspace and output folder will not be read properly.  This is easily fixed by placing an r character prior to the string to make it a raw string and treat the foreslash as a raw string not an escape character.

arcpy.env.workspace = r"C:/Users/Desktop/weekly_avg_project"

OutputFolder = r"C:/Users/Desktop/weekly_avg_project/output"

You have a similar problem with escape characters in the following lines:

inNCfiles = arcpy.env.workspace + "/" + filename

outRaster = OutputFolder + "/" + fileroot

they would need to be double foreslashes to end up with a single foreslash character in python(or even better, use os.path.join to concatenate filepaths instead of just string concatentation. See: https://docs.python.org/2/library/os.path.html)

inNCfiles = arcpy.env.workspace + "//" + filename

outRaster = OutputFolder + "//" + fileroot

Hopefully this gets you closer to the script working, if not repost your updated code with the code formatting and new error message and we can go from there.

0 Kudos
DuminduJayasekera
New Contributor III

Thanks Ian.

I appreciate your reply. I have considered your comments and modified in the script. (Please see the script and error above). My code still does not work. I get a different error now.  I really appreciate if you can help me the script working. 

Thanks a lot again. 

0 Kudos
IanMurray
Frequent Contributor

Code looks much better now thanks.

Issue now is how your output coordinate system parameter is formatted, as it will not accept the string name for the coordinate system nor the WKID for the coordinate system

"The coordinate system to which the input raster will be projected. The default value is set based on the Output Coordinate System environment setting.

Valid values for this parameter are

  • A file with the ".prj" extension (the prj files which ship with ArcGIS can be found in "C:\Program Files\ArcGIS\Coordinate Systems").
  • An existing feature class, feature dataset, raster catalog (basically anything with a coordinate system).
  • The string representation of a coordinate system. These lengthy strings can be generated by adding a coordinate system variable to ModelBuilder, setting the variable's value as desired, then exporting the model to a Python script."

You have a few options, one is to have a .prj file somewhere with the output spatial reference you want and use that as the parameter(http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/project-raster.htm).  Similarly if you had a raster you manually projected into that coordinate system somewhere in your workspace, you could reference that as an input.

Another option I believe is to set the Output Coordinate System environment setting in your script to as it should then become the default value for the output coordinate system variable(you can do this by placing arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4269) with the rest of your environment variables(Actually I'm curious if this works, someone else had a similar error here and I wonder if you left the 3rd parameter blank the environment setting would entirely take care of things.)

0 Kudos
DuminduJayasekera
New Contributor III

Thanks again Ian. I have tried what you mentioned in the above comment. (Please see the code and errors above). I would like to know how I can copy the string representation of a coordinate system using the model builder as you noted in 3rd bullet. 

I have opened the properties of the raster of same coordinate system but it does not provide string representation of a coordinate system. See the image below.

Again your help is highly appreciated to make the script working and moving forward. 

0 Kudos
IanMurray
Frequent Contributor

It looks like the environmental setting isn't working at all so ignore that.  I feel at this point it would just be easiest to write the file path to a dataset or .prj on your computer that has the coordinate system you are wanting to use instead of that long drawn string for the coordinate system.

If you do want the coordinate system string, I saw you had another thread similar to this using Model Builder, you can export your model to a python script from Model Builder and as long as you set the coordinate system properly on the Project Raster Tool in Modelbuilder, it should have the correct string for you to use in python for that coordinate system.

0 Kudos
DuminduJayasekera
New Contributor III

OK. Thanks again for taking your time and reply.

I was able to export to the model builder to obtain the projection string as you can see in the code. I ran the model and it still gives me an error saying that "Undefined coordinate system for input dataset". The input file is a netCDF file which does not have projection information. How can i solve this problem and make the script running.?

Please see the attached sample of netCDF files I am using. 

Thanks and appreciate your help again. 

0 Kudos
IanMurray
Frequent Contributor

You have the same input for Make NetCDF Layer and for Project Raster(filename or the .nc file), while you should take the output from Make NetCDF Layer (an actual raster layer, which should have a coordinate system) and use that for the input of Project Raster.  The output of Project Raster would need a different name than you are currently using for the input.

DuminduJayasekera
New Contributor III

Thanks a bunch again. I made it work. What you suggested was correct. I made the changes and it's correctly projecting the files. Thanks a lot again. 

0 Kudos