3

EDIT: So the mod asked me to narrow down my question (waiting for answers to tell me to give more info, but here we go).

How can i change my piece of python code to do the following steps:

#importing input rasters
inRas_fdr = arcpy.Raster(in_fdr)
inRas_lu = arcpy.Raster(in_lu)
inRas_ns = arcpy.Raster(in_ns)
inRas_out = arcpy.Raster(in_out)
inRas_riv = arcpy.Raster(in_riv)
inRas_slp = arcpy.Raster(in_slp)

#converting input each input raster into a numpy array 
arr_fdr = arcpy.RasterToNumPyArray(inRas_fdr)
arr_lu = arcpy.RasterToNumPyArray(inRas_lu)
arr_ns = arcpy.RasterToNumPyArray(inRas_ns)
arr_out = arcpy.RasterToNumPyArray(inRas_out)
arr_riv = arcpy.RasterToNumPyArray(inRas_riv)
arr_slp = arcpy.RasterToNumPyArray(inRas_slp)


#the slope layer is the closest to the DEM and therefore considered most likely to be correct
#x and y dimensions and cell size are taken from the slope layer
lines=(len(arr_slp[1]))
cols=(len(arr_slp[:]))
print("the dimension of your data is " + str(lines) + " x " + str(cols) + " cells")

cellsize=inRas_slp.meanCellWidth
print("the cell size of your data is: " + str(cellsize))

xextent=lines*cellsize
yextent=cols*cellsize
print("the x-extent of your data is: " + str(xextent))
print("the y-extent of your data is: " + str(yextent))

arr_x = numpy.arange(0,lines,1)
arr_y = numpy.arange(0,cols,1)

#starting editing of the netcdf dataset
#all commands concerning the dataset have to be contained inbetween
#start of the editing (netCDF4.Dataset() and close()) 
foo = netCDF4.Dataset("foo.nc", "w", format="NETCDF4",clobber=True)

#adding dimensions
foo.createDimension("t", 1)
foo.createDimension("x", lines)
foo.createDimension("y", cols)

#adding variables

var_fdr=foo.createVariable("fdr",int,("t","x","y"),zlib=False)
var_lu=foo.createVariable("lu",int,("t","x","y"),zlib=False)
var_ns=foo.createVariable("ns",int,("t","x","y"),zlib=False)
var_out=foo.createVariable("out",int,("t","x","y"),zlib=False)
var_riv=foo.createVariable("riv",int,("t","x","y"),zlib=False)
var_slp=foo.createVariable("slp","f4",("t","x","y"),zlib=False)
var_x=foo.createVariable("x",int,("x",))
var_y=foo.createVariable("y",int,("y",))
var_t=foo.createVariable("t",int,("t",))

print(var_fdr)

#populating dimensions
t = 1
var_x[:] = arr_x
var_y[:] = arr_y

#populating variables with values from the numpy arrays created earlier
var_fdr[:]=arr_fdr
var_lu[:]=arr_lu
var_ns[:]=arr_ns
var_out[:]=arr_out
var_riv[:]=arr_riv
var_slp[:]=arr_slp

#adding global attributes
foo.description = "Dataset containing all the input variables for an MMF-based erosion model"

#local attributes
var_fdr.standard_name="Flow direction"
var_lu.standard_name="Land use"
var_ns.standard_name="Precipitation"
var_out.standard_name="Catchment outlet"
var_riv.standard_name="River cells"
var_slp.standard_name="Slope"
var_x.standard_name="x-coordinate"
var_y.standard_name="y-coordinate"
var_x.axis="X"
var_y.axis="Y"
var_x.units="meter"
var_y.units="meter"
var_t.units="years"

foo.close()
  1. Loading several raster layers (.flt format) with the same extent into python

at the moment i have 6 input layers: fdr, lu, ns, out, riv, slp

  1. Storing each of the layers into a numpy array
  2. Setting up the netCDF database

having 3 dimensions: time, x, y (where I need only one timestep at the moment)

having 6 variables: fdr, lu, ns, out, riv, slp

  1. Writing the original raster value into one variable of the netCDF database

The netCDF file should have 3 dimensions: time,

  1. Conserving spatial extent of the input layers

i.e. CRS and, more importantly, cell size of 5m.

So far, the output that I get does hold the values of the input layer, yet the cells are

  1. arranged in a strange and stripewise way

link "out1.png" in the original post below (i cannot post more than 8 links..)

it should look like the link "multiband_fd.png"

  1. of cell size = 1m

I just don't understand whether it is a problem with the missing projection or with either the creation or the assignment of the dimensions/variables.


ORIGINAL POST

Using ArcGIS 10.4.1 and Python 2.7.10

Python packages used: os, sys, arcpy, numpy, ctypes, netCDF4

I am but a beginner with both python and ArcGis.


I am going crazy with the following task: I prepared input data (raster, .flt format) for an erosion model written in FORTRAN. The model is being revised at the moment and the programmer doing the job has proposed using the netCDF-format instead of multiple .flt layers both as input and output. I had never heard about this format until two days ago and am struggling with converting the .flt layers while leaving their values and spatial extent unchanged. At least until the revision, the model used to demand raster layers of the exact same extent (albeit containing loads of noData values) with the raster cells lying on top of each other. I want to keep the raster size at 5m.

My input layers are the following:

  • fdr: flow direction (output from ArcGIS flow direction tool: value=2,4,8,..,128)
  • slp: slope (in degrees, output of ArcGIS slope tool)
  • lu: land use (integers from 1 to 9, representing different land use classes)
  • ns: precipitation (value=1 for the whole area, allocating all of the catchment area to the same precipitation gauge)
  • out: catchment outlet: a single raster cell with value=1 for the outlet (the rest of the layer are NoData values)
  • riv: "river", flow paths above a threshold of 5ha have value=1 (output from ArcGIS flow accumulation tool), the rest of the layer are NoData values

  1. My first approach was trying to use the ArcGis tool "Raster to netCDF", wich failed due to the fact that is only possible to put in ONE input layer. I tried to use a raster catalog consisting of my input, yet ArcGIS refuses to take it as input. Anyway, my ArcGIS tends to crash whenever handling raster catalogs, so I gave this up.

  2. Building a multiband Raster from my input rasters with the ArcGIS tool "Composite Bands", then using "Raster to netCDF" on this multiband Raster. This works fine and also conserves the spatial extent/environments. However, it seems to be not possible at all (?) to assign labels to the bands, they will be named "Band 1" etc and you have to have a table at hand to assign them all the time.

  3. I realised the most likely way to accomplish my goal was to use python to 1) import my input rasters, 2) transform them into (numpy)-arrays, 3) create a netCDF databse and 4) write the data contained in the arrays into the created netCDF database


Following approach 3), the attached python script "script1.py" is what I came up with (having a hard time finding information on the internet). I guess it starts getting interesting only at line 135.

When the output netCDF-file "foo.nc" is loaded into ArcGIS, I only ever got this weird striped pattern (map output "out1.png"). In the current state, not even the values of the variables are right anymore, they used to be in a former version of my script (i.e. flow direction 2,4,8 etc.), yet being totally misplaced and strangely warped. The lighter areas inbetween are filed with NoData. Map output "multiband_fd.png" shows the output from the script when using the multiband approach for flow directions - this is what it should look like. Another problem is how to get to my grid size of 5m - for i do not see how this is at all possible when using my supposed method.

script1.py - my glorious python script

out1.png - output of my script in ArcGIS

foo.nc - netCDF output file of my script

multiband_fd.png - output of my script, working fine

mbr.nc - netCDF file for the multiband output

Can anybody tell, maybe just from looking at the output "out1.png", whether this might be a problem with projections, assigning of dimensions (t, x, y) or something else?

I also put the script output into the tool "Panoply" which accesses and visualizes the netCDF data. I compared it to a preliminary output netCDF file that the programmer working on the model gave to me. I tried to reproduce his variable/dimension layout. "panoply_foo.png" shows my script output "panoply_merk1.png" and "panoply_merk2.png" show the output file given to me I realized that he is using netCDF3 instead of netCDF4 as I do.

panoply_foo.png - panoply output of my script

panoply_merk1.png - panoply output of the netCDF file given to me by the programmer

panoply_merk2.png - panoply output of the netCDF file given to me by the programmer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.