You should be able to use the method you were describing above, but create your vector layer using the 'vector grid' tool under the vector>research tools menu.
This has a checkbox to align the extents and resolution of the vector grid to a raster layer, so should give you a perfect match with your original raster, and also remember to check the option to output as polygons (it's set to lines by default)!
You should then be able to calculate the area of each cell and convert back to a raster as you described above.
Since you put GRASS into the tags, here's a solution based on GRASS:
First, you need to know exactly what coordinate system the original data is in (as always with GRASS). I see that the *.prj file contains "TRANSVERSE MERCATOR" but it's not one of the standard UTM zones. Since you have not mentioned the EPSG code of this data, here's the proj4 string for creating a new matching GRASS Location:
+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs
Now start GRASS, create a LOCATION (based on the above) and a MAPSET.
And now import the shapefile into GRASS by running:
v.in.ogr input="sp_data.shp" output=sp_data
Add a column for the area of each feature, and calculate its area:
v.db.addcolumn map=sp_data columns="area_sqm INT"
v.to.db sp_data option=area column=area_sqm unit=meters
Now convert to raster, using the area_sqm column as the raster value.
But first set the GRASS computational region:
g.region -p vector=sp_data res=100
Now create a vector grid of polygons at the raster resolution:
v.mkgrid map=grid
Intersect with the forest polygons, and calculate areas of these intersect polygons:
Intersect with the forest polygons, dissolve the polygons inside the same grid cell, calculate areas of these intersect polygons, and move the area value to the grid cells:
v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect
v.dissolve --overwrite input=forest_grid_intersect@jn column=a_cat output=forest_grid_intersect_dissolve
v.db.addtable map=forest_grid_intersect_dissolve columns="area_sqm INT"
v.to.db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters
v.db.join map=grid@jn column=cat other_table=forest_grid_intersect_dissolve other_column=cat
v.db.addcolumn map=forest_grid_intersect columns="area_sqm INT"
v.to.db forest_grid_intersect option=area column=area_sqm unit=meters
and convert to raster, using the area attribute as the raster value:
v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm
And you're done. Hope I have it right this time :-) .
Best Answer
As pointed out in this geonet discussion (https://geonet.esri.com/thread/43270), which I think the OP was referring to, you could project your data to an cylindrical Equal Area projection. Your raster would then have all cells with same size regardless of the proximity to the equator. You could look at the cell size in the raster layer properties (the unit would be the one of your chosen projection), convert it to hectares, calculate the area of a cell (cellsize^2), and eventually perform zonal statistics ad you need. Of course, there's surely be some distortion, you have to choose carefully the projection which fits your needs the most, or subdivide your area in many parts and project each of them in a different projection (but this would be very tricky).