Dear all,
About a month ago I emailed this group to ask for advice on projecting a lot
a netCDF files. Dave Allured advised me to use NCL (thank you for that).
Since I am fairly new to the netCDF community and have little experience
with its tools (including Java), I decided to give it a shot with ArcGIS 9.2
that has some netCDF capabilities. One of the advantages of ArcGIS is that
it knows very well the various coordinate systems and has automated routines
for projecting. Especially, ArcGIS can help deal with Sphere/Spheroid
projections. ArcGIS can resample data as well. Unfortunately it also loses
the metadata.
The following code is a Python script that imports multiple netCDF files,
saves them, project them, resamples them and exports them as new netCDF
files. It uses the Spatial Analyst tool box (requires special license) of
ArcGIS 9.2 as well as others (don?t require special licenses).
I just thought I should share,
Cedric
#
---------------------------------------------------------------------------
# projectStIV.py
# Created on: Fri Aug 03 2007 01:29:06 PM
# (generated by ArcGIS/ModelBuilder)
#
---------------------------------------------------------------------------
# Import system modules
import sys, string, os, arcgisscripting
# Create the Geoprocessor object
gp = arcgisscripting.create()
# Check out any necessary licenses
gp.CheckOutExtension("spatial")
# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management
Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Multidimension
Tools.tbx")
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst
Tools.tbx")
gp.overwriteoutput = 1
for year in ["2006"]:
for month in ["05"]:
#for day in
["01","02","03","04","05","06","07","08","09","10","11","12","13","14","15",
"16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","
31"]:
for day in
["18"]:#,"19","20","21","22","23","24","25","26","27","28","29","30","31"]:
#for hour in
["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14",
"15","16","17","18","19","20","21","22","23"]:
for hour in ["23"]:
# Local variables...
netcdf_1 = "G:\\Research\\GIS\\Toolboxes\\StIVfiles\\ST4." +
year + month + day + hour +".01h.nc"
precip = "precip"
precip_layer_lyr =
"G:\\Research\\GIS\\Toolboxes\\StIVfiles\\precip_layer.lyr"
precip_Raster_img =
"G:\\Research\\GIS\\Toolboxes\\StIVfiles\\precip_Raster.img"
projected_img =
"G:\\Research\\GIS\\Toolboxes\\StIVfiles\\projected.img"
Resampled_img =
"G:\\Research\\GIS\\Toolboxes\\StIVfiles\\Resampled.img"
Mask =
"G:\\Research\\GIS\\forNoah\\Processed_fromCedric\\rasters\\hgt_900m"
Extracted_img =
"G:\\Research\\GIS\\Toolboxes\\StIVfiles\\Extracted.img"
netcdf_2 = "G:\\Research\\GIS\\Toolboxes\\StIVfiles\\" +
year + month + day + hour + "_project.nc"
print "Starting the projection of " + "ST4." + year + month
+ day + hour +".01h.nc"
# Process: Make NetCDF Raster Layer...
gp.MakeNetCDFRasterLayer_md(netcdf_1, "Total_precipitation",
"x", "y", precip, "", "", "BY_VALUE")
# Process: Save To Layer File...
gp.SaveToLayerFile_management(precip, precip_layer_lyr)
# Process: Copy Raster...
gp.CopyRaster_management(precip_layer_lyr,
precip_Raster_img, "", "", "", "NONE", "NONE", "")
# Process: Define Projection...
gp.DefineProjection_management(precip_Raster_img,
"PROJCS['stereographic_wgs84',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHER
OID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degre
e',0.0174532925199433]],PROJECTION['Stereographic'],PARAMETER['False_Easting
',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-105.0],
PARAMETER['Scale_Factor',0.933],PARAMETER['Latitude_Of_Origin',90.0],UNIT['K
ilometer',1000.0]]")
# Process: Project...
gp.ProjectRaster_management(precip_Raster_img,
projected_img,
"PROJCS['PCS_LCC_Guad(GCS_NAD83)',GEOGCS['GCS_North_American_1983',DATUM['D_
North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['G
reenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Confor
mal_Conic'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],P
ARAMETER['Central_Meridian',-98.19],PARAMETER['Standard_Parallel_1',28.4],PA
RAMETER['Standard_Parallel_2',30.3],PARAMETER['Latitude_Of_Origin',29.34],UN
IT['Meter',1.0]];IsHighPrecision", "NEAREST", "4763",
"NAD_1983_To_WGS_1984_1", "",
"PROJCS['stereographic_wgs84',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHER
OID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degre
e',0.0174532925199433]],PROJECTION['Stereographic'],PARAMETER['False_Easting
',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-105.0],
PARAMETER['Scale_Factor',0.933],PARAMETER['Latitude_Of_Origin',90.0],UNIT['K
ilometer',1000.0]];IsHighPrecision")
# Process: Resample...
gp.Resample(projected_img, Resampled_img, "900", "NEAREST")
# Process: Extract by Mask...
gp.ExtractByMask(Resampled_img, Mask, Extracted_img)
# Process: Raster to NetCDF...
gp.RasterToNetCDF_md(Extracted_img, netcdf_2,
"Total_precipitation", "", "x", "y", "", "")
print "Projection completed, file created: " + year + month
+ day + hour + "_project.nc"
