Writing a Raster to a Personal Geodatabase

[BobCave]

Hello,

I am writing code that will export elevation data to a Geodatabase raster.  When I test the program writing output to a File Geodatabase, the output raster is correct.  Using the same test case, and writing to a Personal Geodatabase, the output raster contains only 'no data' values.  I haven't been able to find anything in the documentation that there are differences that I need to be aware of when writing to File vs. Personal Geodatabases.  Does anyone know of anything? 

I am testing using ArcGIS 10.1, but haven't tried this on other versions.

Here is the code I am using:

HRESULT hr;
IRasterWorkspaceExPtr rasterWs;  // Already opened for file or personal geodatabase

IRasterStorageDef2Ptr rasterStorageDef(CLSID_RasterStorageDef);
hr = rasterStorageDef->put_Tiled(VARIANT_FALSE);
IPntPtr cellSize(CLSID_Pnt);
hr = cellSize->SetCoords(aXSpacing, aYSpacing); // spacing values are passed in
hr = rasterStorageDef->putref_CellSize(cellSize);

IRasterDefPtr rasterDef(CLSID_RasterDef);
hr = rasterDef->putref_SpatialReference(aSpatialRef); // spatial ref is passed in and defined correctly

IRasterDataset3Ptr rasterDs;
hr = rasterWs->CreateRasterDataset("temp", 1, PT_FLOAT, rasterStorageDef, ("sde.DEFAULT"), rasterDef, NULL, &rasterDs);

IRasterPtr raster;
hr = rasterDs->CreateFullRaster(&raster);

IRasterPropsPtr rasterProps(raster);
hr = rasterProps->put_Height(aPixHeight); // Height and width are passed in
hr = rasterProps->put_Width(aPixWidth);

CComVariant varNoData;
hr = rasterProps->get_NoDataValue(&varNoData);
float defaultNoData = -FLT_MAX;
CComVariant varDefaultNoData(defaultNoData);    
hr = rasterProps->put_NoDataValue(varDefaultNoData);                

esriGeometry::IEnvelopePtr envelope(esriGeometry::CLSID_Envelope);
double minX, maxX, minY, maxY; // Values passed in
hr = envelope->PutCoords(minX, minY, maxX, maxY);
hr = rasterProps->put_Extent(envelope);

IPntPtr blockSize(CLSID_Pnt);
hr = blockSize->SetCoords((double)aPixWidth, (double)aPixHeight);
IPixelBlockPtr pixelBlock;
hr = raster->CreatePixelBlock(blockSize, &pixelBlock);

CComVariant varPixelArray;
hr = pixelBlock->get_SafeArray(0, &varPixelArray);
COleSafeArray pixels(varPixelArray);
long index[2];

// Get the elevation values
for(index[1] = 0; index[1] < aPixHeight; index[1]++)
{
    // Write out each elevation in the row
    float theElev;

    for(index[0] = 0; index[0] < aPixWidth; index[0]++)
    {
        if ( elevation value at pixel location is valid )
        {
            pixels.PutElement(index, &theElev);
        }
        else
        {
            pixels.PutElement(index, &noData);
        }
    }
}

hr = pixelBlock->put_SafeArray(0, pixels.Detach());

IRasterEditPtr rasterEdit(raster);
IPntPtr tlc(CLSID_Pnt);
hr = tlc->SetCoords(0.0, 0.0);
hr = rasterEdit->Write(tlc, pixelBlock);
hr = rasterEdit->Refresh();

ISaveAsPtr saveAs(raster);
IDatasetPtr newDs;
IWorkspacePtr ws(rasterWs);
hr = saveAs->SaveAs("table", ws, CComBSTR("GDB"), &newDs);

Does anyone see anything that looks wrong with respect to how ArcObjects works with File vs. Personal Geodatabases?

Another question related to writing a raster to a geodatabase:  I am having trouble getting the elevation values stored in the geodatabase.  The solution I came up with is to create a temporary raster dataset, use it to create a raster and write the elevation values to the raster, then use ISaveAs to write the raster to the geodatabase. (This is illustrated in the snippet above.) When the process is done, the code deletes the temporary raster dataset. 

Is there a better way to do this?

Thanks,

Bob

[ErikaDade]

I am following the code in this thread, but I am unable to write rasters to File GDB (have written GRID rasters successfully for years).

I cannot get the nodata value to stick in the output raster (blank in properties) and the data values are reported as nodata when interrogating with the Identify tool in ArcMap.  With GRID rasters, it was necessary to compute statistics on the raster before writing.  Something is clearly wrong here because the code fails when computing statistics. If I skip that step then I get a garbage raster.

Any ideas would be much appreciated!

(All result code testing removed for brevity)

HRESULT WriteGDBRasterTest(CString rasterParentDir,
     CString rasterName, 
     CString& sDatatype,
     int nCols, int nRows, double Xmin, double Ymin, 
     float pixWidth, float pixHeight,
     float noDataValue,
     float* rData)
{
 AFX_MANAGE_STATE(AfxGetStaticModuleState());
 HRESULT hr = S_OK;

 CString sRaster = ExtractFileName(rasterName);
 CString sWorkspaceDir = rasterParentDir;
 long Width = nCols;
 long Height = nRows;
 double cellWidth = (double)pixWidth;
 double cellHeight = (double)pixHeight;
 double dXmin = (double)Xmin;
 double dYmin = (double)Ymin;
 double dXmax = dXmin + cellWidth*Width;
 double dYmax = dYmin + cellHeight*Height;

 IPointPtr ipPointOrigin(CLSID_Point);
 ipPointOrigin->put_X(dXmin);
 ipPointOrigin->put_Y(dYmin);

 IPntPtr ipPntCellsize(CLSID_Pnt);
 ipPntCellsize->put_X(cellWidth);
 ipPntCellsize->put_Y(cellHeight);

 IWorkspaceFactoryPtr ipWorkspaceFactory;
 IWorkspacePtr ipWorkspace;
 IRasterWorkspaceExPtr ipRasterWorkspaceEx;

 hr = ipWorkspaceFactory.CreateInstance(CLSID_FileGDBWorkspaceFactory);

 CComBSTR bWorkspaceDir = sWorkspaceDir.AllocSysString();
 hr = ipWorkspaceFactory->OpenFromFile(bWorkspaceDir, 0, &ipWorkspace);

 hr = ipWorkspace->QueryInterface(IID_IRasterWorkspaceEx, (void**)&ipRasterWorkspaceEx);

 VARIANT_BOOL bPermanent = TRUE;
 CComBSTR bsRaster = sRaster.AllocSysString();
 IRasterDatasetPtr ipRasterDataset;
 IRasterStorageDef3Ptr ipRasterStorageDef(CLSID_RasterStorageDef);
 hr = ipRasterStorageDef->putref_CellSize(ipPntCellsize);
 hr = ipRasterStorageDef->putref_Origin(ipPointOrigin);

 CString ConfigKeyword = "";
 CComBSTR bsConfigKeyword(ConfigKeyword);
 IRasterDefPtr ipRasterDef(CLSID_RasterDef);

 IGeometryDefPtr ipGeomDef = 0;
 rstPixelType rasterPixelType =  PT_FLOAT;

 // Start with temp raster because need to use SaveAs to write raster correctly to FileGDB?
 CComBSTR bsTempRaster = "temp";
 hr = ipRasterWorkspaceEx->CreateRasterDataset(bsTempRaster, 1,  rasterPixelType, ipRasterStorageDef, 
      bsConfigKeyword, ipRasterDef, ipGeomDef, &ipRasterDataset);
 bsConfigKeyword.Empty();

 ipRasterDef = 0;
 ipRasterStorageDef = 0;
 ipPntCellsize = 0;
 ipPointOrigin = 0;

 IRasterDataset3Ptr ipRasterDataset3(ipRasterDataset);
 IRasterPtr ipRaster;
 hr = ipRasterDataset3->CreateFullRaster(&ipRaster);

 // If set props on raster then raster cursor works properly (but cannot set nodata)
 IRasterPropsPtr ipRasterProps_fromRaster;
 hr = ipRaster->QueryInterface(IID_IRasterProps, (void **) &ipRasterProps_fromRaster);
 //// If set on dataset then raster cursor will NOT work  (but can set nodata)
 //IRasterPropsPtr ipRasterProps_fromDataset;
 //IRasterBandPtr ipRasterBand;
 //IRasterBandCollectionPtr ipRasterBandCollection;
 //hr = ipRasterDataset->QueryInterface(IID_IRasterBandCollection, (void **)&ipRasterBandCollection);
 //hr = ipRasterBandCollection->Item(0, &ipRasterBand);
 //hr = ipRasterBand->QueryInterface(IID_IRasterProps, (void **) &ipRasterProps_fromDataset);

 ipRasterProps_fromRaster->put_PixelType(rasterPixelType);
 //ipRasterProps_fromDataset->put_PixelType(rasterPixelType);

 IEnvelopePtr ipEnv(CLSID_Envelope);
 ipEnv->PutCoords(dXmin, dYmin, dXmax, dYmax);
 ipRasterProps_fromRaster->put_Extent(ipEnv);
 ipRasterProps_fromRaster->put_Height((long)Height);
 ipRasterProps_fromRaster->put_Width((long)Width);
 CComVariant nodata(noDataValue);

 //hr = ipRasterProps_fromDataset->put_NoDataValue(nodata);
 hr = ipRasterProps_fromRaster->put_NoDataValue(nodata);
 ::VariantClear(&nodata);

 // Why isn't nodata set in output layer????  

 //CComVariant nodata2;
 //hr = ipRasterProps_fromDataset->get_NoDataValue(&nodata2);
 //float test2 = nodata2.fltVal;
 //VARTYPE vType = nodata2.vt;     // This returns correctly {-9999., R4}

 CComVariant nodata3((float)0);  // Initializes vtype R4
 hr = ipRasterProps_fromRaster->get_NoDataValue(&nodata3);  
 float test3 = nodata3.fltVal;  // Comes back as safearray of R4, (with correct value in element 0)
 VARTYPE vType3 = nodata3.vt;   //  8196 = 0x2000 + 0x0004 (SafeArray + R4)

 IRaster2Ptr ipRaster2(ipRaster);
 IRasterCursorPtr ipRasterCursor;
 hr = ipRaster2->CreateCursorEx(0, &ipRasterCursor);  // Use default of full width x 128 rows

 // Set up raster for writing
 IRasterEditPtr ipRasterEdit(ipRaster2);

 IPixelBlockPtr ipPixelBlock;
 IPixelBlock3Ptr ipPixelBlock3;
 long pixelBlockWidth = 0;
 long pixelBlockHeight = 0;
 IPntPtr ipTLC;
 
 VARIANT_BOOL vNotFinished;

 // Rotate through blocks using Raster Cursor
 do {
  hr = ipRasterCursor->get_PixelBlock(&ipPixelBlock);

  hr = ipPixelBlock->QueryInterface(IID_IPixelBlock3, (void**)&ipPixelBlock3);

  hr = ipPixelBlock3->get_Width(&pixelBlockWidth);  // Should be same as ncols

  hr = ipPixelBlock3->get_Height(&pixelBlockHeight);  // Should be lesser of nrows and 128

  hr = ipRasterCursor->get_TopLeft(&ipTLC);

  double dTLRow;
  hr = ipTLC->get_Y(&dTLRow);
  int iTLRow = (int) dTLRow;

  CComVariant varArray; 
  ::VariantInit(&varArray);

  hr = ipPixelBlock3->get_PixelDataByRef(0, &varArray);  // Single band

  SAFEARRAY *saArray = *(varArray.pparray); 
  hr = SafeArrayLock(saArray);

  long long ij = 0;
  long long lWidth = static_cast<long long>(Width);
  
  for(int i = 0; i < pixelBlockHeight; i++) {
   for(int j = 0; j < Width; j++)  {
    long ind[2];
    ind[0]=j; 
    ind[1]=i; 
    ij = (i+iTLRow)*lWidth + j;
    float val = rData[ij];
    hr = ::SafeArrayPutElement(saArray, ind, (void*)&val);
   }
  }

  hr = SafeArrayUnlock(saArray);

  // Write out the new raster block
  hr = ipRasterEdit->Write(ipTLC, ipPixelBlock);
  hr = ipRasterEdit->Refresh();

  // Advance the Raster Cursor
  hr = ipRasterCursor->Next(&vNotFinished);

  ipPixelBlock3 = 0;
  ipPixelBlock = 0;

  ::VariantClear(&varArray);

 } while (vNotFinished);
 ipRasterCursor = 0;
 ipRasterEdit = 0;

 // New method to save to File GDB
 ISaveAs2Ptr ipSaveAs(ipRaster2);
 IDatasetPtr ipDS;
 IWorkspacePtr ipWS(ipRasterWorkspaceEx);
 hr = ipSaveAs->SaveAs(bsRaster, ipWS, CComBSTR("GDB"), &ipDS);
 bsRaster.Empty();
 ipDS = 0;
 ipSaveAs = 0;
 ipWS = 0;

 //// Calculate RasterBand stats. MUST do to complete the band.(Was true for GRID anyway)
 //hr = ipRasterBand->ComputeStatsAndHist(); // E_FAIL HERE  (Also, NoData not getting set correctly)

 //IRasterStatisticsPtr ipStats;
 //hr = ipRasterBand->get_Statistics(&ipStats);

 //double rMin, rMax;
 //ipStats->get_Minimum(&rMin);
 //ipStats->get_Maximum(&rMax);
 //ipStats = 0;

 // MUST release these pointers BEFORE creating and adding layer to map
 ipRasterProps_fromRaster = 0;
 //ipRasterProps_fromDataset = 0;
 //ipRasterBand = 0;
 //ipRasterBandCollection = 0;
 ipRaster = 0;
 ipRaster2 = 0;
 ipRasterDataset = 0;
 ipRasterDataset3 = 0;
 ipRasterWorkspaceEx = 0;
 ipWorkspace = 0;

 // Must create new workspace and rasterdataset connections to get layer
 hr = ipWorkspaceFactory->OpenFromFile(bWorkspaceDir, 0, &ipWorkspace);
 bWorkspaceDir.Empty();
 hr = ipWorkspace->QueryInterface(IID_IRasterWorkspaceEx, (void**)&ipRasterWorkspaceEx);
 bWorkspaceDir.Empty();

 // Delete temp raster
 hr = ipRasterWorkspaceEx->DeleteRasterDataset(bsTempRaster);
 bsTempRaster.Empty();

 // Open saved raster
 IRasterDatasetPtr ipRDS;
 CComBSTR bRName = sRaster.AllocSysString();
 hr = ipRasterWorkspaceEx->OpenRasterDataset(bRName, &ipRDS);
 bRName.Empty();

 IRasterLayerPtr ipRasterLayer(CLSID_RasterLayer);
 hr = ipRasterLayer->CreateFromDataset(ipRDS);

 //  Add Raster Layer to current map
 IMapPtr ipMap; 
 hr = m_ipMxDoc->get_FocusMap(&ipMap);

 hr = ipMap->AddLayer(ipRasterLayer);

 // Refresh TOC
 IActiveViewPtr ipActiveView;
 hr = m_ipMxDoc->get_ActiveView(&ipActiveView);
 hr = ipActiveView->Refresh();
 hr = ipActiveView->ContentsChanged();

 // MUST Release smart pointers
 ipRasterLayer = 0;
 ipActiveView = 0;
 ipMap = 0;
 ipRasterWorkspaceEx = 0;
 ipWorkspace = 0;
 ipWorkspaceFactory = 0;
 ipRDS = 0;

 return hr;
}