nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

calculate # of features with distance and time constraints

QiuhuaMa

Hi,

Earlier I asked how to find the nearest fire that burned 7 years ago and get helped.

Here is the detail the discussion:

https://community.esri.com/message/491873?et=watches.email.thread#491873

Another way to check the effect of fire on house is to count number and size of fire within a certain distance in the last 7 years.

For example, I want to find number of fires within 5km of each house in the past 7 years, then calculate their average size!

How to code this in Python?

My sample data is attached. SHAPE_AERA is the size of fire.

Thanks so much,

Chelsea

datatest_update (2).zip

Find more posts tagged with

Accepted answers

XanderBakker

I guess I should have read you question better... there is no distance tolerance for the fires. However this is relatively easy to do. See the code below. I just did a run, but less than 200 house have at most a single fire within a distance of 5km. I am running a test with 25km to check the results.

Code:

The search distance can be specified on line 61
Since your data uses foot as linear unit and your search distance is metric, line 62 contains the conversion from feet to meter
Line 102 checks to see if the distance is within the search distance
minimum distance and oid of nearest fire is provided anyway.
To make it faster you could make a selection (select by location) of the fires that are with the search distance and use that as input.

import arcpy
from datetime import datetime, date, timedelta
from dateutil.relativedelta import relativedelta

def main():
    # input fc's
    fc_house = r"D:\Xander\GeoNet\DistFireHouse\shp\housesaledata.shp"
    fc_fire = r"D:\Xander\GeoNet\DistFireHouse\shp\firehistory_NAD_fips3002_feet.shp"

    # input field
    fld_date_sale = "ma_date_dt"

    # output fields
    fld_dist = "Dist2Fire"
    fld_fire_oid = "FireOID"

    # additional fields
    fld_FireMinArea = "FireMinAr"
    fld_FireMaxArea = "FireMaxAr"
    fld_FireMeanArea = "FireMeanAr"
    fld_FireCount = "FireCount"

    # add fields
    addField(fc_house, fld_dist, "DOUBLE")
    addField(fc_house, fld_fire_oid, "LONG")

    # add additional fields
    addField(fc_house, fld_FireMinArea, "DOUBLE")
    addField(fc_house, fld_FireMaxArea, "DOUBLE")
    addField(fc_house, fld_FireMeanArea, "DOUBLE")
    addField(fc_house, fld_FireCount, "LONG")

    # loop through houses
    cnt = 0
    flds = ("SHAPE@", fld_date_sale, fld_dist, fld_fire_oid, fld_FireMinArea,
            fld_FireMaxArea, fld_FireMeanArea, fld_FireCount)
    with arcpy.da.UpdateCursor(fc_house, flds) as curs:
        for row in curs:
            cnt += 1
            if cnt % 25 == 0:
                print "Processing house: {0}".format(cnt)
            pnt = row[0]
            date_sale = row[1]

            results = getNearestFireAndAreaStats(fc_fire, pnt, date_sale)
            row[2] = results[0] # dist nearest
            row[3] = results[1] # oid nearest
            row[4] = results[2] # Min Area
            row[5] = results[3] # Max Area
            row[6] = results[4] # Mean Area
            row[7] = results[5] # Count Fires
            curs.updateRow(row)


def addField(fc, fldname, fldtype):
    if len(arcpy.ListFields(fc, wild_card=fldname)) == 0:
        arcpy.AddField_management(fc, fldname, fldtype)

def getNearestFireAndAreaStats(fc, pnt, date_sale):
    # maximum distance for fire
    search_dist = 5000.0 # meters
    ft_to_m = 0.30480061 # conversion feet to meter

    # input fields
    fld_year = "FIRE_YEAR"
    fld_month = "FIRE_MONTH"
    fld_day = "FIRE_DAY"
    fld_area = "SHAPE_AREA"

    # defaults
    min_dist = 999999
    oid_fire = -1
    min_area = 9999999999
    max_area = 0
    sum_area = 0
    mean_area = 0
    cnt_fires = 0

    # results
    res = []
    flds = ("SHAPE@", "OID@", fld_year, fld_month, fld_day, fld_area)
    with arcpy.da.SearchCursor(fc, flds) as curs:
        for row in curs:
            polygon = row[0]
            oid = row[1]
            year = int(row[2])
            try:
                month = int(row[3])
            except:
                month = 6
            try:
                day = int(row[4])
            except:
                day = 15

            if chkDate(date_sale, year, month, day):
                # these are the fires we are looking for
                dist = pnt.distanceTo(polygon)
                if dist < min_dist:
                    min_dist = dist
                    oid_fire = oid
                if (dist * ft_to_m) < search_dist:
                    area = row[5]
                    if cnt_fires == 0:
                        min_area = area
                        max_area = area
                        sum_area = area
                    else:
                        if area < min_area: min_area = area
                        if area > max_area: max_area = area
                        sum_area += area
                    cnt_fires += 1

    if cnt_fires > 0:
        mean_area = sum_area / float(cnt_fires)
    else:
        min_area = 0
    return [min_dist * ft_to_m, oid_fire, min_area, max_area, mean_area, cnt_fires]


def chkDate(date_sale, year_fire, month_fire, day_fire):
    date_fire = datetime(year=year_fire, month=month_fire, day=day_fire)
    date_2monthsbefore_sale = date_sale + relativedelta(months=-2)
    year_sale = date_sale.year
    test = (year_fire >= year_sale - 7) and (date_fire < date_2monthsbefore_sale)
    return test

if __name__ == '__main__':
    main()

All comments

XanderBakker

Is it an empty geometry? In that case I suppose it can't calculate distances and errors will occur. I you exclude it from the dataset it should work. If not post it and I will have a look.

Kind regards, Xander

QiuhuaMa

I think I figured out the problem. Input Field Shape_area is 0 in my new

data.

thanks,

qiuhua

XanderBakker

Are those results based on the dataset you posted earlier or are you using new set of data? If so please post a sample of this data to be able to reproduce the problem.

QiuhuaMa

Xander,

is there a problem with the code?

I want to calculate number of fires and average size of fires within 5km of a houses in the past 7 years. I saw that number of fire is 1 but max, min and average size of fire is 0.

if there is only one fire, these sizes should be the size of that fire, right?

QiuhuaMa

Yes, you are right! Most houses in my data are located in the city and far

away from the fires. With 25km, 70% of the houses should have more than one

fire!

Qiuhua

XanderBakker

Code:

The search distance can be specified on line 61
Since your data uses foot as linear unit and your search distance is metric, line 62 contains the conversion from feet to meter
Line 102 checks to see if the distance is within the search distance
minimum distance and oid of nearest fire is provided anyway.
To make it faster you could make a selection (select by location) of the fires that are with the search distance and use that as input.

import arcpy
from datetime import datetime, date, timedelta
from dateutil.relativedelta import relativedelta

def main():
    # input fc's
    fc_house = r"D:\Xander\GeoNet\DistFireHouse\shp\housesaledata.shp"
    fc_fire = r"D:\Xander\GeoNet\DistFireHouse\shp\firehistory_NAD_fips3002_feet.shp"

    # input field
    fld_date_sale = "ma_date_dt"

    # output fields
    fld_dist = "Dist2Fire"
    fld_fire_oid = "FireOID"

    # additional fields
    fld_FireMinArea = "FireMinAr"
    fld_FireMaxArea = "FireMaxAr"
    fld_FireMeanArea = "FireMeanAr"
    fld_FireCount = "FireCount"

    # add fields
    addField(fc_house, fld_dist, "DOUBLE")
    addField(fc_house, fld_fire_oid, "LONG")

    # add additional fields
    addField(fc_house, fld_FireMinArea, "DOUBLE")
    addField(fc_house, fld_FireMaxArea, "DOUBLE")
    addField(fc_house, fld_FireMeanArea, "DOUBLE")
    addField(fc_house, fld_FireCount, "LONG")

    # loop through houses
    cnt = 0
    flds = ("SHAPE@", fld_date_sale, fld_dist, fld_fire_oid, fld_FireMinArea,
            fld_FireMaxArea, fld_FireMeanArea, fld_FireCount)
    with arcpy.da.UpdateCursor(fc_house, flds) as curs:
        for row in curs:
            cnt += 1
            if cnt % 25 == 0:
                print "Processing house: {0}".format(cnt)
            pnt = row[0]
            date_sale = row[1]

            results = getNearestFireAndAreaStats(fc_fire, pnt, date_sale)
            row[2] = results[0] # dist nearest
            row[3] = results[1] # oid nearest
            row[4] = results[2] # Min Area
            row[5] = results[3] # Max Area
            row[6] = results[4] # Mean Area
            row[7] = results[5] # Count Fires
            curs.updateRow(row)


def addField(fc, fldname, fldtype):
    if len(arcpy.ListFields(fc, wild_card=fldname)) == 0:
        arcpy.AddField_management(fc, fldname, fldtype)

def getNearestFireAndAreaStats(fc, pnt, date_sale):
    # maximum distance for fire
    search_dist = 5000.0 # meters
    ft_to_m = 0.30480061 # conversion feet to meter

    # input fields
    fld_year = "FIRE_YEAR"
    fld_month = "FIRE_MONTH"
    fld_day = "FIRE_DAY"
    fld_area = "SHAPE_AREA"

    # defaults
    min_dist = 999999
    oid_fire = -1
    min_area = 9999999999
    max_area = 0
    sum_area = 0
    mean_area = 0
    cnt_fires = 0

    # results
    res = []
    flds = ("SHAPE@", "OID@", fld_year, fld_month, fld_day, fld_area)
    with arcpy.da.SearchCursor(fc, flds) as curs:
        for row in curs:
            polygon = row[0]
            oid = row[1]
            year = int(row[2])
            try:
                month = int(row[3])
            except:
                month = 6
            try:
                day = int(row[4])
            except:
                day = 15

            if chkDate(date_sale, year, month, day):
                # these are the fires we are looking for
                dist = pnt.distanceTo(polygon)
                if dist < min_dist:
                    min_dist = dist
                    oid_fire = oid
                if (dist * ft_to_m) < search_dist:
                    area = row[5]
                    if cnt_fires == 0:
                        min_area = area
                        max_area = area
                        sum_area = area
                    else:
                        if area < min_area: min_area = area
                        if area > max_area: max_area = area
                        sum_area += area
                    cnt_fires += 1

    if cnt_fires > 0:
        mean_area = sum_area / float(cnt_fires)
    else:
        min_area = 0
    return [min_dist * ft_to_m, oid_fire, min_area, max_area, mean_area, cnt_fires]


def chkDate(date_sale, year_fire, month_fire, day_fire):
    date_fire = datetime(year=year_fire, month=month_fire, day=day_fire)
    date_2monthsbefore_sale = date_sale + relativedelta(months=-2)
    year_sale = date_sale.year
    test = (year_fire >= year_sale - 7) and (date_fire < date_2monthsbefore_sale)
    return test

if __name__ == '__main__':
    main()

QiuhuaMa

Xander,

Thanks so much! however I cannot see which line restrict distance, # of

fires and average size of fires within 5kms of a house in the past 7 years!

Qiuhua

XanderBakker

Hi Qiuhua Ma ,

Have a look at the changed code below. To run it, make sure you reference your shapefiles on lines 07 and 08.

Changes are:

additional output field names (lines 17 to 21)
adding the additional fields (lines 27 to 31)
on line 35 and 36 add the fields to the list of fields used in the cursor
line 44 calls a function getNearestFireAndAreaStats (defined on line 58) which replaces the function getNearestFire. This new function returns a list with 6 results.
Lines 45-50 reads the results and sets them to the current row
On line 58 the function getNearestFireAndAreaStats starts:
- additional variables on lines 67 - 71
- lines 92 - 101 determines the min, max, sum and count of the relevant fires
- line 108 - 109 calculates the mean area (if there are fires found)
- line 110 returns the list of results for the current house

import arcpy
from datetime import datetime, date, timedelta
from dateutil.relativedelta import relativedelta

def main():
    # input fc's
    fc_house = r"D:\Xander\GeoNet\DistFireHouse\shp\housesaledata.shp"
    fc_fire = r"D:\Xander\GeoNet\DistFireHouse\shp\firehistory_NAD_fips3002_feet.shp"

    # input field
    fld_date_sale = "ma_date_dt"

    # output fields
    fld_dist = "Dist2Fire"
    fld_fire_oid = "FireOID"

    # additional fields
    fld_FireMinArea = "FireMinAr"
    fld_FireMaxArea = "FireMaxAr"
    fld_FireMeanArea = "FireMeanAr"
    fld_FireCount = "FireCount"

    # add fields
    addField(fc_house, fld_dist, "DOUBLE")
    addField(fc_house, fld_fire_oid, "LONG")

    # add additional fields
    addField(fc_house, fld_FireMinArea, "DOUBLE")
    addField(fc_house, fld_FireMaxArea, "DOUBLE")
    addField(fc_house, fld_FireMeanArea, "DOUBLE")
    addField(fc_house, fld_FireCount, "LONG")

    # loop through houses
    cnt = 0
    flds = ("SHAPE@", fld_date_sale, fld_dist, fld_fire_oid, fld_FireMinArea,
            fld_FireMaxArea, fld_FireMeanArea, fld_FireCount)
    with arcpy.da.UpdateCursor(fc_house, flds) as curs:
        for row in curs:
            cnt += 1
            if cnt % 25 == 0:
                print "Processing house: {0}".format(cnt)
            pnt = row[0]
            date_sale = row[1]
            results = getNearestFireAndAreaStats(fc_fire, pnt, date_sale)
            row[2] = results[0] # dist nearest
            row[3] = results[1] # oid nearest
            row[4] = results[2] # Min Area
            row[5] = results[3] # Max Area
            row[6] = results[4] # Mean Area
            row[7] = results[5] # Count Fires
            curs.updateRow(row)


def addField(fc, fldname, fldtype):
    if len(arcpy.ListFields(fc, wild_card=fldname)) == 0:
        arcpy.AddField_management(fc, fldname, fldtype)

def getNearestFireAndAreaStats(fc, pnt, date_sale):
    # input fields
    fld_year = "FIRE_YEAR"
    fld_month = "FIRE_MONTH"
    fld_day = "FIRE_DAY"

    # defaults
    min_dist = 999999
    oid_fire = -1
    min_area = 9999999999
    max_area = 0
    sum_area = 0
    mean_area = 0
    cnt_fires = 0

    # results
    res = []

    flds = ("SHAPE@", "OID@", fld_year, fld_month, fld_day)
    with arcpy.da.SearchCursor(fc, flds) as curs:
        for row in curs:
            polygon = row[0]
            oid = row[1]
            year = int(row[2])
            try:
                month = int(row[3])
            except:
                month = 6
            try:
                day = int(row[4])
            except:
                day = 15

            if chkDate(date_sale, year, month, day):
                # these are the fires we are looking for
                if cnt_fires == 0:
                    min_area = polygon.area
                    max_area = polygon.area
                    sum_area = polygon.area
                else:
                    if polygon.area < min_area: min_area = polygon.area
                    if polygon.area > max_area: max_area = polygon.area
                    sum_area += polygon.area
                cnt_fires += 1

                dist = pnt.distanceTo(polygon)
                if dist < min_dist:
                    min_dist = dist
                    oid_fire = oid

    if cnt_fires > 0:
        mean_area = sum_area / float(cnt_fires)
    return [min_dist, oid_fire, min_area, max_area, mean_area, cnt_fires]


def chkDate(date_sale, year_fire, month_fire, day_fire):
    date_fire = datetime(year=year_fire, month=month_fire, day=day_fire)
    date_2monthsbefore_sale = date_sale + relativedelta(months=-2)
    year_sale = date_sale.year
    test = (year_fire >= year_sale - 7) and (date_fire < date_2monthsbefore_sale)
    return test

if __name__ == '__main__':
    main()

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.