POST
|
Hello, I have a series of aerial photos of a flood. They are not orthorectified and have basically no metadata even though I know roughly the location. There are simple pictures takes from an airplane. I would like to superimpose my flood shape file to them. Any idea about how it can be done? Thanks Fabio
... View more
07-24-2015
02:38 AM
|
0
|
1
|
3004
|
POST
|
Hello, thanks a lot for replying. Attached the file (named ras1.tif instead of line1.tif). This is the code: def GetPixelValueForPercentile(dctper, percentile):
"""will return last pixel value
where percentile LE searched percentile."""
pix_val = sorted(dctper.keys())[0] # initially assign lowest pixel value
for k in sorted(dctper.keys()):
perc = dctper
if perc <= percentile:
pix_val = k
else:
break
return pix_val
def FieldExist(featureclass, fieldname):
"""Check if field exists"""
import arcpy
fieldList = arcpy.ListFields(featureclass, fieldname)
return len(fieldList) == 1
import arcpy, sys
from arcpy import env
from arcpy.sa import *
env.workspace = r"C:\Users\Flood Depth"
f = open('log_file.txt', 'w')
# Enable Spatial analyst
arcpy.CheckOutExtension("Spatial")
# settings
ras = r"C:\Users\Dataset\Flood_Depth\LiD_MA_Sal_cm.tif"
fc = r"C:\Users\Dataset\Flood_Depth\C_Saletto.shp"
# ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster"
# fc = r"C:\Path\To\Folder\Or\Geodatabase\With\PolygonFeatureclass"
lst_perc = [2, 5, 10, 90, 95, 98] # list of percentiles to be calculated
fld_prefix = "Perc_"
f.write("DEM File input: " + ras + "\n")
f.write("Shapefile input: " + fc + "\n" + "\n")
# arcpy.env.workspace = "IN_MEMORY"
arcpy.env.overwriteOutput = True
# add fields if they do not exist in fc
# be aware, field that exist will be overwritten!
print "filling field list and adding fields"
flds = []
for perc in lst_perc:
fld_perc = "{0}{1}".format(fld_prefix, perc)
flds.append(fld_perc)
if not FieldExist(fc, fld_perc):
arcpy.AddField_management(fc, fld_perc, "LONG")
print "flds={0}".format(flds)
# loop through polygons
print "loop through polygons"
flds.append("SHAPE@")
i = 0
with arcpy.da.UpdateCursor(fc, flds) as curs:
for row in curs:
i += 1
print "Processing polygon: {0}".format(i)
f.write("Processing polygon: {0}".format(i) + "\n")
polygon = row[flds.index("SHAPE@")]
# Execute ExtractByMask (you can't send the polygon object)
try:
print " - ExtractByMask..."
ras_pol = ExtractByMask(ras, polygon)
outname = "ras{0}".format(i) + ".tif"
ras_pol.save(outname)
print " - saved raster as {0}".format(outname)
except:
print arcpy.GetMessages(2) #"Unexpected error:", sys.exc_info()[0]
# create dictionary with value vs count
print " - fill dict with Value x Count"
flds_ras = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(outname, flds_ras)}
# calculate number of pixels in raster
print " - determine sum"
cnt_sum = sum(dct.values())
print " - sum={0}".format(cnt_sum)
# loop through dictionary and create new dictionary with val vs percentile
print " - create percentile dict"
dct_per = {}
cnt_i = 0
for val in sorted(dct.keys()):
cnt_i += dct[val]
dct_per[val] = float(cnt_i) / float(cnt_sum)
sort_dct_per = sorted(dct_per)
#print sort_dct_per
print dct_per
# loop through list of percentiles
print " - iterate perceniles"
f.write(" - iterate perceniles" + "\n")
for perc in lst_perc:
# use dct_per to determine percentiles
perc_dec = float(perc) / 100
print " - Perc_dec for is {0}".format(perc_dec)
f.write(" - Perc_dec for is {0}".format(perc_dec) + "\n")
pixval = GetPixelValueForPercentile(dct_per, perc_dec)
print " - Perc for {0}% is {1}".format(perc, pixval)
f.write(" - Perc for {0}% is {1}".format(perc, pixval)+ "\n")
# write pixel value to percentile field
print " - Store value"
fld_perc = "{0}{1}".format(fld_prefix, perc)
row[flds.index(fld_perc)] = pixval
# update row
print " - update row"
curs.updateRow(row)
del(outname)
# return SA license
arcpy.CheckInExtension("Spatial")
f.close()
... View more
07-06-2015
09:11 AM
|
0
|
3
|
800
|
POST
|
Hi everyone, I'm running a script that it's basically zonal statistic over polylines (I modified one script that I found in this forum Zonal statistics - percentiles ). The code takes a shape file containing the lines along which I want to analyse the elevation. For every line, it creates a raster (let's say it converts the line into a raster) where in every pixel it writes the corresponding elevation value. It analyses the values of the elevation along the line and returns the max, min and percentiles values. It repeats the same thing for all the lines in the shape file. The code is running without problems with a certain digital elevation model. Recently I received a new LiDAR file and I tried to run the code but I got an error: RuntimeError: cannot open "line1.tif" line1.tif is the line converted to raster. I checked the file and it is exactly the line in tif format with the right values of elevation. The folder in which it is stored, is the same as using the previous DEM. The DEM and LiDAR have different resolution and therefore line1.tif in case of DEM in something like 2 KB and in case of LIDAR around 60 KB. This is the only difference I can notice. Anyone has a clue of what is happening? Thanks a lot Fabio
... View more
07-06-2015
07:39 AM
|
0
|
6
|
5681
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|