Extracting Pixel Data

StephenFricke · ‎08-07-2012

Hey Everyone,

I am working on a script right now that extracts climate data for one particular pixel over a period of time. I have written a script which works, but I feel like there is a much faster and efficient way to go about it. Right now I have a script which downloads a NetCDF file, creates a numpy 2d array, converts this array to a raster, then I create a point which is fed x and y coordinates, and the raster is then clipped by this point so that I am left with the pixel in the raster that this point falls within. I feel like an easier way to do this is if I could specify the xy coordinate of the pixel of interest, and then simply read this pixel value from the numpy array. This way I wouldn't be creating a numpy array for entire climate dataset (the entire United States) and then converting this huge dataset to a raster. I would be very grateful if anyone knows how to do this. Thanks!

Here is a bit of my code where I create the point which is used to clip the climate raster to my pixel of interest:

def ClipPix(output,pixelX,pixelY):
    arcpy.CreateFolder_management(output,'/clipped climate rasters')
    pointFC= arcpy.CreateFeatureclass_management(output, "Point.shp", "POINT")
    inPoint = arcpy.CreateObject("Point")
    inPoint.X = float(pixelX)
    inPoint.Y = float(pixelY)
    cursor= arcpy.InsertCursor(pointFC)
    feature= cursor.newRow()
    feature.shape= inPoint
    cursor.insertRow(feature)
    del cursor
    arcpy.env.workspace= output+'/climate rasters'
    ClimRasters= arcpy.ListRasters()
    for ClimRaster in ClimRasters:
        arcpy.Clip_management(output+'/climate rasters/'+ClimRaster, "", output+"/clipped climate rasters/"+ClimRaster, pointFC, "", "ClippingGeometry")

ChristopherThompson · ‎08-08-2012

I haven't worked with netCDF data but it looks like you could avoid a lot of what you are doing by using the Make NetCDF Table View (arcpy.MakeNetCDFTableView_md) - the help/tutorial for this show how you can create a table view for a single location, which is what you want. I guess i'm curious how you are obtaining the coordinates of your pixel in the first place - do you have multiple locations you are workinig through to build your dataset? Ultimately, i'm curious what the output is you need to work with - is it something you can show in a map? the tabular data? Lots of questions I know but those help understand the context of what you're trying to accomplish at the end of the day.

StephenFricke · ‎08-08-2012

I haven't worked with netCDF data but it looks like you could avoid a lot of what you are doing by using the Make NetCDF Table View (arcpy.MakeNetCDFTableView_md) - the help/tutorial for this show how you can create a table view for a single location, which is what you want. I guess i'm curious how you are obtaining the coordinates of your pixel in the first place - do you have multiple locations you are workinig through to build your dataset? Ultimately, i'm curious what the output is you need to work with - is it something you can show in a map? the tabular data? Lots of questions I know but those help understand the context of what you're trying to accomplish at the end of the day.

Hey, so my ultimate goal is to create a tool that can access climate data at the smallest resolution possible. I have climate datasets for the entire United States. I would like a tool where you could enter in the x and y coordinates for any location in the United States and then be able to extract climate information about the particular pixel within the climate dataset in which that x and y coordinate lies.

ChristopherThompson · ‎08-08-2012

i think you could easily wrap that tool into a python function that includes the netCDF file and the x,y coordinate pair as arguments - from there you could probably create a feature layer based on the x,y data or a table. so that might look something like:

import arcpy
climvar = 'xxxx' #variable to contain the climate variable you want to use
def get_dat(netcdf,pixelX,pixelY):
    tbl = arcpy.MakeNetCDF(netcdf,climvar,'',[[lon,{pixelY}],[lat,{pixelX}]])
    arcpy.CopyRows_management(tbl,'path to saved table.. ')

again, this is highly conceptual but seems like it would work given the information in the tutorial for working with netcdf data that you can see here:
http://resources.arcgis.com/en/help/main/10.1/index.html#/Exercise_2_Creating_a_temperature_table_at...

The help is for 10.1 and shows the entries you'd make in the GUI for the tool but those elements are all available for scripting as well.

StephenFricke · ‎08-09-2012

Hey thanks again for the reply. Ideally, I don't want to have to download a netcdf file because it takes a while to download. I would simply like to be able to read a netcdf file and create a numpy array at a specific pixel. Not sure how to do this though. Here is my script that creates a numpy array and then a raster for the entire USA. I would really like to know a way to specify the array to simply be the one pixel at a specific x and y coordinate.

def AddBcsd(output,GcmModel,GcmScenario,GcmMonth,GcmVar,GcmYear):
    arcpy.CreateFolder_management (output, "climate rasters")
    dataset =('http://climatedata_NetCDF/')
    fn = ('dcsbcsd_'+str(GcmModel)+'_'+str(GcmScenario)+'_'+str(GcmMonth)+'_'+str(GcmVar)+'_run1_conus_epscor.nc')
    bcsd= open_url(dataset + fn)
    arcpy.AddMessage("\nDownloading Climate Data for "+str(GcmMonth)+"/"+str(GcmYear)+"...")
    lt = bcsd['lat']
    ln = bcsd['lon']
    lat,lon = np.meshgrid(ln.data,lt.data)
    vars = bcsd.keys()
    years = range(1900,2100)
    data = bcsd[vars[1]][years.index(int(GcmYear)),:,:]
    indata = np.float32(data)
    array= np.array(indata).reshape(621,1405, order='F').copy()
    ycellsize = .041667938#float(lat[0][0] - lat[1][0])
    xcellsize = .041671753#float(abs(lon[0][0] - lon[0][1]))
    llc = arcpy.Point(float(lat.min().tolist())-(ycellsize/2),float(lon.min().tolist())-(xcellsize/2))
    myRaster = arcpy.NumPyArrayToRaster(array,llc,xcellsize,ycellsize,-9999)
    myRaster.save(output+'/climate rasters/'+str(GcmMonth)+'_'+str(GcmYear))
    coordinateSystem = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    arcpy.DefineProjection_management(myRaster, coordinateSystem)

ChristopherThompson · ‎08-09-2012

so you are reading the url directly in these lines:

dataset =('http://climatedata_NetCDF/')
fn = ('dcsbcsd_'+str(GcmModel)+'_'+str(GcmScenario)+'_'+str(GcmMonth)+'_'+str(GcmVar)+'_run1_conus_epscor.nc')
bcsd= open_url(dataset + fn)

Could you not supply the variable bcsd as the netCDF argument in this line:

tbl = arcpy.MakeNetCDF(netcdf,climvar,'',[[lon,{pixelY}],[lat,{pixelX}]])

instead that might look like:

tbl = arcpy.MakeNetCDF(bcsd,climvar,'',[[lon,{pixelY}],[lat,{pixelX}]])

I feel like you're using the numpy array to create a dataset you may not need to create in this instance as there are tools already for creating vector and raster data from a netCDF file, or a table which is the suggestion above. Did you take a look at these tools to see if they will work for you?

What I'm not sure of is if you can reference a file through the url_open method you are using.. if this works currently then it may well work for what you need. Then the issue is how to supply the pixelX and pixelY values - where do those come from or how are they derived? If you truly have a single or even a finite known set of locations you need this data for then that's easy. Discovering them through other means.. well, thats a different set of issues and tangential to this thread.

StephenFricke · ‎08-09-2012

so you are reading the url directly in these lines:
dataset =('http://climatedata_NetCDF/')
fn = ('dcsbcsd_'+str(GcmModel)+'_'+str(GcmScenario)+'_'+str(GcmMonth)+'_'+str(GcmVar)+'_run1_conus_epscor.nc')
bcsd= open_url(dataset + fn)
Could you not supply the variable bcsd as the netCDF argument in this line:
tbl = arcpy.MakeNetCDF(netcdf,climvar,'',[[lon,{pixelY}],[lat,{pixelX}]])
instead that might look like:
tbl = arcpy.MakeNetCDF(bcsd,climvar,'',[[lon,{pixelY}],[lat,{pixelX}]])
I feel like you're using the numpy array to create a dataset you may not need to create in this instance as there are tools already for creating vector and raster data from a netCDF file, or a table which is the suggestion above. Did you take a look at these tools to see if they will work for you?

What I'm not sure of is if you can reference a file through the url_open method you are using.. if this works currently then it may well work for what you need. Then the issue is how to supply the pixelX and pixelY values - where do those come from or how are they derived? If you truly have a single or even a finite known set of locations you need this data for then that's easy. Discovering them through other means.. well, thats a different set of issues and tangential to this thread.

Hey thanks for the suggestion. The problem is, I am having an error pop up when I use the bcsd as the input netcdf. I don't know the specifics of the error, all it says is "RuntimeError: Object: Error in executing tool"... The pixelY and pixelX are simply GetParameterAsText values which the user specifies on the scripts user interface.

ChristopherThompson · ‎08-13-2012

Hey thanks for the suggestion. The problem is, I am having an error pop up when I use the bcsd as the input netcdf. I don't know the specifics of the error, all it says is "RuntimeError: Object: Error in executing tool"... The pixelY and pixelX are simply GetParameterAsText values which the user specifies on the scripts user interface.

Yeah, I gave that a whirl and can't seem to make it work either.. you may have to download the file, even to a scratch location that you can clear out later. Otherwise, have you looked at the NetCDF module that is part of the numpy package? you might be able to get what you are looking for using those tools.

StephenFricke · ‎08-13-2012

Yeah, I gave that a whirl and can't seem to make it work either.. you may have to download the file, even to a scratch location that you can clear out later. Otherwise, have you looked at the NetCDF module that is part of the numpy package? you might be able to get what you are looking for using those tools.

Hey thanks for trying it out. I am going to research more into Python's NetCDF4 module.