Convert shapefile to raster

Sep 12, 2014 at 6:52 AM
Hi, I'm new to DotSpatial.

I'm now trying to convert a shape file with lots of polygons to raster with special extent and cell size. and the raster pixel value should be calculated from the field of the polgons.

I'll be very grateful if you could answer with your code.
Sep 17, 2014 at 11:55 AM
You should have a look at the DotSpatial.Analysis.VectorToRaster.ToRaster method and the parameters it requires and the work is done

Oscar
Marked as answer by frozencookie on 9/23/2014 at 12:02 AM
Sep 18, 2014 at 1:37 AM
Thanks ~!!! I'll try!
Sep 18, 2014 at 7:06 AM
Oscarafone77 wrote:
You should have a look at the DotSpatial.Analysis.VectorToRaster.ToRaster method and the parameters it requires and the work is done

Oscar
Hi

This method doesn't work correctly. I special a cell size of 5, but the raster I got has a cell height of 4.69xxxx, and a cell width of 4.99xxxxx.

Here is my code:
            string shpFile = @"E:\Project\TestData\UpCatchment.shp";
            double xMin = 82292.591;
            double yMin = 81549.981;
            double xMax = 82851.763;
            double yMax = 81599.988;
            double cellsize = 5;
            IFeatureSet fs = FeatureSet.Open(shpFile);
            fs.FillAttributes();
            fs.DataTable.Columns.Add("myField");
            for (int i = 0; i < fs.Features.Count; i++)
            {
                fs.Features[i].DataRow.BeginEdit();
                fs.Features[i].DataRow["myField"] = 1;
                fs.Features[i].DataRow.EndEdit();
            }
            Extent extent = new Extent(new double[] {xMin, yMin, xMax,yMax});
            IRaster raster = VectorToRaster.ToRaster(fs, cellsize, "myField", "output.bgd", "", new string[] { }, null);
Sep 18, 2014 at 8:47 AM
Edited Sep 18, 2014 at 8:48 AM
Hi,

try to pass also the extent.

If you download dotspatial sources and open the vectortoraster.cs file you find the following code in the toraster method:
        public static IRaster ToRaster(IFeatureSet fs, Extent extent, double cellSize, string fieldName,
                                       string outputFileName,
                                       string driverCode, string[] options, IProgressHandler progressHandler)
        {
            Extent env = extent;
            if (cellSize == 0)
            {
                if (env.Width < env.Height)
                {
                    cellSize = env.Width / 256;
                }
                else
                {
                    cellSize = env.Height / 256;
                }
            }
            int w = (int)Math.Ceiling(env.Width / cellSize);
            if (w > 8000)
            {
                w = 8000;
                cellSize = env.Width / 8000;
            }
            int h = (int)Math.Ceiling(env.Height / cellSize);
            if (h > 8000)
            {
                h = 8000;
            }
           ......
you should test the following instructions to understand if the Toraster method will return to you the results you are looking for:
int w = (int)Math.Ceiling(env.Width / cellSize);
int h = (int)Math.Ceiling(env.Height / cellSize);
I didn't test it but I guess that
since your envelope height and width might be not a multiple of the cellsize you specified, it might turn out that the newly created raster will not have the discretization you want.

Have a look into the ToRaster method and see what it does to understand it. The cellsize is never used after the evaluation of w and h, hence I guess that the output raster cellsize is re-evaluated from w,h and the envelope.

A simple test would be to pass an envelope which has a width and height exactly multiple of the cellsize you want. Then the results should be what you are looking for.

Hope I've been clear

Oscar
Sep 19, 2014 at 2:20 AM
hi,

this is weird. I tried to test the instructions. Here is my code:
            Extent extent = new Extent(new double[] {xMin, yMin, xMax,yMax});

            int w = (int)Math.Ceiling(extent.Width / cellsize);
            int h = (int)Math.Ceiling(extent.Height / cellsize);

            IRaster raster = VectorToRaster.ToRaster(fs, cellsize, "unUsed", "output.bgd", "", new string[] { }, null);

            int rw = raster.NumColumns;
            int rh = raster.NumRows;
            int crw = (int)Math.Ceiling(extent.Width / raster.CellWidth);
            int crh = (int)Math.Ceiling(extent.Height / raster.CellHeight);
And I get :

w = 112, h =11, rw = 135, rh = 8, crw = 112, crh = 11

the output.bdg has the sanme size of my input polygon, not the extent I specified. :(
Sep 19, 2014 at 2:52 PM
check if the Extent you are building is the same as the extent of fs

to create the raster the method uses fs.extent

Oscar
Sep 20, 2014 at 1:56 AM
Hi,

the extent is not the same as fs. I want to convert the feature to raster with special extent. In fact, I have several shape files. I want to convert them to rasters with a same extent.
Sep 22, 2014 at 11:51 AM
Maybe you can do a work around by coupling a raster conversion with a clip raster method

First you convert your shapes to raster, then clip those rasters with a feature of desired extent.

To clip a raster with a polygon I found this old sub I made a long ago (maybe there is one already implemented in dotspatial)

I don't think it is very optimized and also I don't remember what I was doing exactly, I guess that if you are clipping a raster with a larger polygon it adds nodata value on the outer border...but check it

if you want to have a look, if I remember it was pretty slow

    Public Shared Function ClipGrid2Featureset(inputRast As IRaster, inputFS As Envelope, cellWidth As Double, cellHeight As Double, Filename As String, cancelProgressHandler As ICancelProgressHandler) As IRaster
        'Validates the input and output data
        If inputRast Is Nothing OrElse inputFS Is Nothing Then
            Return Nothing
        End If

        'Calculate the Intersect Envelope
        Dim envelope1 As IEnvelope = inputRast.Bounds.Extent.ToEnvelope
        Dim envelope2 As IEnvelope = inputFS
        envelope1 = envelope1.Intersection(envelope2)
        Dim envelopr2use As IEnvelope = envelope1

        Dim noOfRow As Integer = Convert.ToInt32(envelope2.Height / cellHeight)
        Dim noOfCol As Integer = Convert.ToInt32(envelope2.Width / cellWidth)

        'create output raster
        Dim output As IRaster
        output = Raster.Create(Filename, inputRast.DriverCode, noOfCol, noOfRow, 1, inputRast.DataType, New String() {""})
        output.NoDataValue = inputRast.NoDataValue
        Dim extOLD As Extent = envelope2.ToExtent

        Dim extNEW As Extent
        Dim dx As Double = ((extOLD.Width) / 2) / noOfCol
        Dim dy As Double = ((extOLD.Height) / 2) / noOfRow
        extNEW = New Extent(extOLD.MinX + dx, extOLD.MinY - dy, extOLD.MaxX + dx, extOLD.MaxY - dy)
        Dim bound As New RasterBounds(noOfRow, noOfCol, extNEW)

        output.Bounds = bound

        output.NoDataValue = inputRast.NoDataValue
        Dim v1 As RcIndex

        For i As Integer = 0 To output.Bounds.NumRows - 1
            For j As Integer = 0 To output.Bounds.NumColumns - 1
                output.Value(i, j) = output.NoDataValue
            Next
        Next

        Dim imin, imax, jmin, jmax As Integer
        Dim xmin, xmax, ymin, ymax As Double
        Dim x1, y1 As Double

        Dim coord2 As New Coordinate
        Dim coord3 As New Coordinate
        Dim i1, j1, i2, j2 As Integer
        Dim rc1i, rc2i As RcIndex
        Dim sizemax, sizeX1, sizeY1, sizemin As Integer
        Dim startX, startY As Integer
        If envelopr2use.IsNull = False Then
            coord2 = envelopr2use.TopLeft
            coord3 = envelopr2use.BottomRight
            rc1i = inputRast.ProjToCell(coord2)
            If rc1i.Column < 0 Then rc1i.Column = 1
            If rc1i.Row < 0 Then rc1i.Row = 1

            rc2i = inputRast.ProjToCell(coord3)
            sizeX1 = rc2i.Column - rc1i.Column + 1
            sizeY1 = rc2i.Row - rc1i.Row + 1

            startX = rc1i.Column
            startY = rc1i.Row
        Else
            startX = 1
            startY = 1
            sizeX1 = 1
            sizeY1 = 1

        End If
        Dim sizemin2 As Integer
        Dim s1_x, s2_y
        sizemin2 = Math.Min(sizeX1, sizeY1)
        sizemin = Math.Min(sizemin2, 1000)
        sizemax = Math.Max(sizeX1, sizeY1)

        For s2 As Integer = startY To startY + sizeY1 - 1 Step sizemin
            For s1 As Integer = startX To startX + sizeX1 - 1 Step sizemin

                s1_x = s1
                s2_y = s2

                Dim inputrast2 As IRaster
                Dim inputrast3 As IRaster

 
                If inputRast.IsInRam = False Then
                    s1_x = Math.Min(s1, inputRast.NumColumns - 1000 - 6)
                    s2_y = Math.Min(s2, inputRast.NumRows - 1000 - 6)
                    inputrast2 = inputRast.ReadBlock(s1_x, s2_y, sizemin + 5, sizemin + 5)
                Else
                    inputrast2 = inputRast
                End If

                xmin = inputrast2.Bounds.Extent.MinX 'inputrast2.Bounds.Extent.MinX
                ymin = inputrast2.Bounds.Extent.MinY
                xmax = inputrast2.Bounds.Extent.MaxX
                ymax = inputrast2.Bounds.Extent.MaxY



                Dim coord1 As New Coordinate

                imin = 0
                For i As Integer = 0 To output.Bounds.NumRows - 1
                    coord1 = output.CellToProj(i, 1)
                    If ymax > coord1.Y Then
                        imin = i
                        Exit For
                    End If
                Next

                imax = -1
                For i As Integer = output.Bounds.NumRows - 1 To 0 Step -1
                    coord1 = output.CellToProj(i, 1)
                    If ymin < coord1.Y Then
                        imax = i
                        Exit For
                    End If
                Next

                jmin = 0
                For j As Integer = 0 To output.Bounds.NumColumns - 1
                    coord1 = output.CellToProj(1, j)
                    If xmin < coord1.X Then
                        jmin = j
                        Exit For
                    End If
                Next

                jmax = -1
                For j As Integer = output.Bounds.NumColumns - 1 To 0 Step -1
                    coord1 = output.CellToProj(1, j)
                    If xmax > coord1.X Then
                        jmax = j
                        Exit For
                    End If
                Next

                Dim previous As Integer = 0
                Dim cellcenter As Coordinate
                Dim max As Integer = (output.Bounds.NumRows + 1)
                For i As Integer = imin To imax
                    For j As Integer = jmin To jmax
                        cellcenter = output.CellToProj(i, j)

                        If output.Value(i, j) = output.NoDataValue Then
                            output.Value(i, j) = inputrast2.GetNearestValue(cellcenter)

                        End If

                    Next

                Next
            Next
        Next

        output.GetStatistics()
        output.Save()
        Return output
    End Function
Oscar
Sep 23, 2014 at 2:50 AM
Hi Oscar,

the clip raster method is a bit slow for what I need. I'm still working on this.

Anyway, Thank you !!!
Sep 23, 2014 at 7:02 AM
Edited Sep 23, 2014 at 7:09 AM
Hi Oscar,

the DotSpatial.Analysis.VectorToRaster.ToRaster method works. I down load the source code and rebuild the DotSpatial.Analysis.dll, then it works correctly. But when I use DotSpatial_Full.1.7, I didn't the correct result.

Thank you, again!!