Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a mapbasic script that divides an irregular polygon into smaller polygons of 1 hectare each.

Would anyone help me to convert the following into an arcpy python script.

'  This program requires an open MapInfo Table named woodlots with fields lot (integer) and area (dec or float)

Declare Sub Main
Declare Function Calc_height(stand1 as object) as float
Declare Sub CookieVarH(ByVal width1, height1 as float,  sliver1 As Object, cookiecutter As Object)
Declare Sub setWoodlot()
Declare Sub setTimberlot()
Declare sub checker
global gs_last_band, gs_project as string

Sub Main 
INCLUDE "C:\program files\MAPINFO\MAPBASIC\MAPBASIC.DEF"
INCLUDE "C:\program files\MAPINFO\MAPBASIC\MENU.DEF"
Include "icons.def" 
Set CoordSys Earth 
Projection 8,33,"m",141,0,0.9996, 500000, 10000000
Dim stand1, rectangle1, band1, band2, sliver1, cookie1, sliver2 As Object 
Dim woodlot1, leftover1, cookie2, right_end1, left_end1 As Object
Dim minX, maxX, minY, maxY, area1, width1, height1, widthend, widthstart As Float
Dim diff_y, row_id As Float
Dim band_no, woodlot_no, max_band, i, lots as Integer
Dim sliver_area As Float
Dim mnX, mnY, mxX, mxY As Float

stand1 = selection.obj


band_no = 1
woodlot_no = 1
gs_last_band = "false"

while CartesianArea(stand1, "hectare") > 1
    minX = ObjectGeography(stand1, OBJ_GEO_MINX)
    minY = ObjectGeography(stand1, OBJ_GEO_MINY)
    maxX = ObjectGeography(stand1, OBJ_GEO_MAXX)
    maxY = ObjectGeography(stand1, OBJ_GEO_MAXY)
    width1 = 1
    area1 = 0
    height1 = 0
    create Rect into variable rectangle1 (minX, minY) (maxX, maxY)
    height1 = calc_height(stand1)                                                           '******   call height
    create Rect into variable band1 (minX, (maxY - height1)) (maxX, maxY)
    sliver1 = overlap(stand1, band1)
    sliver_area = CartesianArea(sliver1, "hectare")
    lots = round((sliver_area/1),1)

'   note " height1 = " + str$(height1) + " sliver_area = " + str$(sliver_area) + " lots = " + str$(lots)


    for i = 1 to lots
        width1 = 10
        area1 = 0
        if i < lots then
            mnX = ObjectGeography(sliver1, OBJ_GEO_MINX)
            mnY = ObjectGeography(sliver1, OBJ_GEO_MINY)
            mxX = ObjectGeography(sliver1, OBJ_GEO_MAXX)
            mxY = ObjectGeography(sliver1, OBJ_GEO_MAXY)
            do 
            create Rect into variable cookie1 ((mnX) , (mnY)) ((mnX + width1), mxY)
                woodlot1 = overlap(sliver1, cookie1)
                area1 = CartesianArea(woodlot1, "hectare")
                width1 = width1 + 0.1
            loop while Area1 < 1                                                ' if area is > 1 put woodlot in layer and remove from sliver
            insert into Woodlots (obj,area)
                values (woodlot1,  area1)
            woodlot_no = woodlot_no + 1

        else                                      'last lot in sliver cut out all
            woodlot1 = sliver1
            area1 = CartesianArea(woodlot1, "hectare")
            if gs_last_band = "false" then
                insert into Woodlots (obj,area)
                    values (woodlot1, area1)
                woodlot_no = woodlot_no + 1
            else
                insert into Woodlots (obj,lot,area)
                values (woodlot1, 9999, area1)              
            call CookieVarH(width1,height1, sliver1, cookie1)
        end if
        end if
        if lots = 1 then
            call CookieVarH(500,height1, sliver1, cookie1)
        end if
        sliver1 = Erase( sliver1, cookie1 )
        sliver_area = CartesianArea(sliver1, "hectare")

    next

create Rect into variable band2 (minX, (maxY - height1)) (maxX, (maxY + 1))
leftover1 = erase( stand1, band2 )
stand1 = leftover1
wend

'       insert remeinder into table
if gs_last_band = "false"     then          '   only do if there is a remeinder
    insert into Woodlots (obj,lot, area)
        values (leftover1, 9999, CartesianArea(stand1, "hectare"))
    stand1 = selection.obj
end if
' clean objects
select * from woodlots
' Objects Snap From Selection Tolerance Vector 1  Units "m" 
commit table woodlots
End Sub    'Main

function Calc_height (stand1 as object) As Float
    Dim sliver_area, trunc, height, stand_area, stand_trunc As Float
    dim bandx, sliver1 as object
    dim lots as float
    Dim mnX, mxX, mnY, mxY As Float
    stand_area = CartesianArea(stand1, "hectare")
    stand_trunc = int(stand_area)
    height = 65
    mnX = ObjectGeography(stand1, OBJ_GEO_MINX)
    mnY = ObjectGeography(stand1, OBJ_GEO_MINY)
    mxX = ObjectGeography(stand1, OBJ_GEO_MAXX)
    mxY = ObjectGeography(stand1, OBJ_GEO_MAXY)
    create Rect into variable bandx (mnX, (mxY - height)) (mxX, mxY)
    sliver1 = overlap(stand1, bandx)
    sliver_area = CartesianArea(sliver1, "hectare")

    lots = Round((sliver_area/1), 1) 
    trunc = int((sliver_area/1))
'   note "sliver_area = " + str$(sliver_area) + "lots = " + str$(lots) + " trunc = " + str$(trunc)
    if (lots = trunc and trunc <> 0 ) then   
'       note "height getting smaller"
        mnX = ObjectGeography(stand1, OBJ_GEO_MINX)
        mnY = ObjectGeography(stand1, OBJ_GEO_MINY)
        mxX = ObjectGeography(stand1, OBJ_GEO_MAXX)
        mxY = ObjectGeography(stand1, OBJ_GEO_MAXY)
        while CartesianArea((sliver1), "hectare")/1 > trunc
            height = height -0.05
            create Rect into variable bandx (mnX, (mxY - height)) (mxX, mxY)
            sliver1 = overlap(stand1, bandx)
        wend
        gs_last_band = "false"  
    elseif (stand_area < 1) or (stand_area < (lots*1)) then   ' last band
    '   note " stand_area = " + str$(stand_area) + "  lots = " + str$(lots)
        height = 300
        gs_last_band = "true"

    else              
    '   note "height getting larger"
    '   if (lots > stand_area) and (lots >1) then 
    '       lots  = lots -1
    '   end if
        if lots = 0 then lots = 1 end if
        gs_last_band = "false"  
        mnX = ObjectGeography(stand1, OBJ_GEO_MINX)
        mnY = ObjectGeography(stand1, OBJ_GEO_MINY)
        mxX = ObjectGeography(stand1, OBJ_GEO_MAXX)
        mxY = ObjectGeography(stand1, OBJ_GEO_MAXY)
        while CartesianArea(sliver1, "hectare")/1 < lots
            height = height + 0.05
            create Rect into variable bandx (mnX, (mxY - height)) (mxX, mxY)
            sliver1 = overlap(stand1, bandx)
        wend
    end if  
   Calc_height = height
    height = 65
'   note " calc_height = " + str$(height)
End Function 


Sub CookieVarH (ByVal width1, height1 as float,  sliver1 As Object, cookiecutter As Object) 
    Dim mnX1, mnY1, mxX1, mxY1 As Float
    mnX1 = ObjectGeography(sliver1, OBJ_GEO_MINX)
    mnY1 = ObjectGeography(sliver1, OBJ_GEO_MINY)
    mxX1 = ObjectGeography(sliver1, OBJ_GEO_MAXX)
    mxY1 = ObjectGeography(sliver1, OBJ_GEO_MAXY)
    create Rect into variable cookiecutter ((mnX1) , (mnY1)) ((mnX1 + width1), mxY1)
End Sub 
share|improve this question
    
Could you please post your MapBasic code so we could help you implement it in Arcpy? –  Fezter Sep 13 '12 at 5:41
    
@Fezter I have edited the question to include the code. Thanks a lot :-) –  Thiru Sep 13 '12 at 6:08

1 Answer 1

I'm going to attempt to answer this in more general way as I believe it's more useful long term to have a bit of a tutorial rather than a straight code dump (also I tend to think you should be paying someone for a straight up code translation, even if it is only a couple of hundred lines of code). If you don't understand Python I'd strongly recommend reading the Python Tutorial as a first step (I'm going to avoid things like loops below as they're describe much better in the tutorial).

  • Python is a dynamically-typed language. This means that you don't have to declare all your variables and functions up front. Functions and Subroutines in Python are the same thing (all functions return null unless declared otherwise).
  • In Python you can bring in other libraries with the import command, e.g., math, json, or in our case arcpy. In general you do this at the start of your script.
"""Split up some polygons"""
#ignore all the declarations

import arcpy

def main():
    #ignore all the declarations her also
    ...
  • In arcpy Spatial References are created from the Spatial Refernce class. See this GIS Stackexchange question + answer (/shameless plug) for a way of creating one from an EPSG code.
  • In ArcGIS feature selections are made on layers. This means that to access a selection through a layer we need to create an arcpy layer object:
lyr = "some_layer_name"
arcpy.MakeFeatureLayer_management("some dataset", lyr, "SQL where clause") #See http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000006p000000
  • To access a feature geometry in ArcGIS we have to use Cursors.
  • In arcpy there is no rectangle object, just a Polygon object. However, each polygon has a hullRectangle property, so we can loop through selected features and grab rectangles directly from the data.
shapeName = arcpy.Describe(lyr).shapeFieldName
for row in arcpy.SearchCursor(lyr, "", spatial_ref):
    rect = row.getValue(shapeName).hullRectangle.split() #hullRectangle returns a string (who knows why)
    point = arcpy.Point()
    array = arcpy.Array()

    for i in range(0, len(rect), 2):
        point.X = coordPair[i]
        point.Y = coordPair[i + 1]
        array.add(point)

    poly = arcpy.Polygon(array)
    print poly.area
    #do something with the polygon
  • By default arcpy geoprocessing tools create feature classes or layers as output. To return a result to a geometry object set the output to an empty geometry object.
g = arcpy.Geometry()
arcpy.Intersect_analysis([poly1, poly2], g)

From this point hopefully you will have enough information to write the Python script as a replacement to the MapBasic function. Please comment if there's not enough information here!

share|improve this answer
    
thank you very much @om_henners for your pointers. i wanted to find out if such a function can be written with python & arcpy. i am thinking of paying someone to write a python script; my requirements are bit different. –  Thiru Sep 18 '12 at 3:42

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.