1

I am trying to create a script using arcpy to automate a serie of task in ArcGIS.

What I am trying to do is to interpolate rainfall data from a region and then, using hydrological basinshapefiles, clip the big raster file into small raster files. After that, I would like to get all the means values from the small basin raster files and save them into a CSV file together with 2 strings, which identify the raster name and the hydrological basin.

The problem that I am facing is to save the numpy arrays that have

[raster_name, day_interpolated, MEAN_value]

. One example of the output would be the following:

espacializacao_2_1994_D1.tif hydrological_basin_name 111,81740045547

Is numpy array the best approach to write these outputs in a CSV file? Because I already saw some tutorials using the "CSV module" to write directly to CSV files. However, I already tried it without success.

I am saying this because this code is going to be repeated a lot of times, something like 19.000 times. Because I am using it to interpolate daily rainfall data from a 54 years time series. So, I guess that an array of this size is not a good approach to solve this problem.

My code is the following:

[![import arcpy
import numpy  as np

from arcpy import sa 
from arcpy.sa import *
from calendar import monthrange


# Set environment settings
arcpy.env.workspace = "C:\\Projetos\\ArcGIS\\Teste9"  

arcpy.env.overwriteOutput = True 


#get the map document
mxd = arcpy.mapping.MapDocument("CURRENT")


#get the data frame
df = arcpy.mapping.ListDataFrames(mxd,"la*")\[0\]





#Months and years to interpolate
mes = \["2"\] #Months
ano = \["1994"\] #Years

#Days to interpolate based in the month lenght
coluna_interpolada_28 = \["D1", "D2", "D3", "D4", "D5", "D6", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14", "D15", "D16", "D17", "D18", "D19", "D20", "D21", "D22", "D23", "D24", "D25", "D26", "D27", "D28"\] 
coluna_interpolada_29 = \["D1", "D2", "D3", "D4", "D5", "D6", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14", "D15", "D16", "D17", "D18", "D19", "D20", "D21", "D22", "D23", "D24", "D25", "D26", "D27", "D28", "D29"\]
coluna_interpolada_30 = \["D1", "D2", "D3", "D4", "D5", "D6", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14", "D15", "D16", "D17", "D18", "D19", "D20", "D21", "D22", "D23", "D24", "D25", "D26", "D27", "D28", "D29", "D30"\]
coluna_interpolada_31 = \["D1", "D2", "D3", "D4", "D5", "D6", "D7", "D8", "D9", "D10", "D11", "D12", "D13", "D14", "D15", "D16", "D17", "D18", "D19", "D20", "D21", "D22", "D23", "D24", "D25", "D26", "D27", "D28", "D29", "D30", "D31"\]



#Interpolation extent
arcpy.env.extent = arcpy.env.workspace + "\\Shapefile\\" + "PB.shp"



#Final list
lista_final = np.array(\[\], dtype = object) 


#Counter that is going to be used to reshape the arrays
contador = 0 


#Loop to go through the time series
#For loop with the years
for i_ano in ano:

    #For loop with the months
    for i_mes in mes:

        #Month Range
        quantidade_dias = monthrange(int(i_ano), int(i_mes))

        #If clauses to define which columns it is goin to interpolate
        if quantidade_dias == (1, 28):
            coluna_interpolada = coluna_interpolada_28
        elif quantidade_dias == (1, 29):
            coluna_interpolada = coluna_interpolada_29
        elif quantidade_dias == (1, 30):
            coluna_interpolada = coluna_interpolada_30
        else:
            coluna_interpolada = coluna_interpolada_31


        #For loop with the days
        for i_dia in coluna_interpolada:

            tabela = i_mes + "_" + i_ano   #Exemplo "2_1994"

            in_Table = arcpy.env.workspace + "\\Dados\\" + tabela + ".csv" #Exemplo "2_1994.csv"
            x_coords = "LONG"
            y_coords = "LAT"
            z_coords = "POINT_Z"
            out_Layer = "espacializacao" + "_" + tabela + "_" + i_dia #nome da camada "lyr" que vai ser criada   
                                                                                        #NOME DO ARQUIVO QUE VAI SER SALVO. Exemplo "espacializacao_2_1994_D2"
            out_Layer_shp = out_Layer + "_shp" #nome da camada "shp" que vai ser criada
            out_Layer_tif = out_Layer + "_tif"
            tamanho_celula = "0,10" #precisar por as aspas, apesar de ser um número
            potencia_idw = "2" #precisa por as aspas, apesar de ser um número
            raio_de_busca = RadiusVariable(12, 12) #Critério para fazer a interpolação, raio variando até 12 quilômetros até conseguir englobar 12 pontos




            # Set the spatial reference
            spRef = arcpy.SpatialReference("WGS 1984")

            #Create event layer
            arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef, "")


            #Exporting event layer as shapefile
            arcpy.FeatureToPoint_management(out_Layer, arcpy.env.workspace + "\\" + "Shapefile\\Exportados\\" + out_Layer_shp + ".shp","")





            #Layer that is going to be deleted "lyr"
            lyr = arcpy.mapping.ListLayers(mxd, "espacializacao",df)



            #Deleting the layer
            for df in arcpy.mapping.ListDataFrames(mxd):
                 for lyr in arcpy.mapping.ListLayers(mxd, "", df): #O parâmetro que não foi passado foi o WildCard, não precisa
                     if lyr.name == out_Layer:
                         arcpy.mapping.RemoveLayer(df, lyr) #Removendo a camada da paleta lateral





            #Some variables to define some parameters to the software
            camada_editando = out_Layer_shp
            coluna_criada = "Media"
            tipo_campo = "FLOAT"
            precisao_campo = ""
            precisao_decimais = ""
            comprimento_campo = 50 #Tamanho qualquer suposto
            arcpy.AddField_management(camada_editando, coluna_criada, "LONG", precisao_campo,
                                          precisao_decimais, comprimento_campo, "", "NULLABLE",
                                          "NON_REQUIRED", "")   
            bacias = \["Abiai", "Camaratuba", "Curimatau", "Gramame", "Guaju", "Jacu", "Mamanguape", "Miriri", "Paraiba", "Piranhas", "Trairi"\]            




            #Code to interpolate
            arcpy.gp.Idw_sa(out_Layer_shp, i_dia, arcpy.env.workspace + "\\Raster\\" + out_Layer_tif + ".tif", tamanho_celula, potencia_idw, "VARIABLE 12 12", "") 


            #Deleting shapefile
            arcpy.Delete_management(out_Layer_shp, "")  






            #For loop to clip the raster file using the shapefiles
            for mascara in bacias:
                importar_camada = arcpy.env.workspace + "\\Shapefile\\Bacias\\" + mascara + ".shp"  #Importing shapefile to clip
                arcpy.MakeFeatureLayer_management(importar_camada, mascara)  



                #Some variables defining some parameters to using in the clip function
                camada_para_recortar = out_Layer_tif + ".tif"
                camada_resultante = out_Layer + "_recortada"
                nome_do_raster = camada_resultante + "_" + mascara + ".tif" 

                #Function to clip the raster file
                arcpy.Clip_management(camada_para_recortar, "", arcpy.env.workspace  + "\\Raster\\Recortes\\" + camada_resultante + "_" + mascara + ".tif", mascara, "-3,402823e+038", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
                media = arcpy.GetRasterProperties_management (nome_do_raster, "MEAN", "")


                lista_strings = np.array(\[out_Layer, mascara\])
                lista_medias = np.array(\[media\])


                #Name of the file to save the means values
                arquivo_com_as_medias = "medias 01" 


                lista_numpy_temporaria = np.append(lista_strings, lista_medias)
                lista_final = np.concatenate((lista_final, lista_numpy_temporaria))


                #Deleting the raster clipped
                arcpy.Delete_management(nome_do_raster, "") 

                #Counter
                contador = contador + 1 

                print(lista_final)

            #Reshaping the file 
            lista_final = lista_final.reshape(contador,3)       

            print(lista_final)

            #Saving the arrays to a CSV file
            np.savetxt(arcpy.env.workspace  + "\\Dados\\Exportados\\" + arquivo_com_as_medias + ".csv", lista_final, fmt="%10s %10s %10s", delimiter=";", header = "")




            #Deleting the original raster
            arcpy.Delete_management(out_Layer_tif + ".tif", "")][1]][1]

I added a picture with the CSV file that I get using this code. As we can see in the picture, it only saves the data related to the last day interpolated. It is overwriting the data.

enter image description here

Saving the arrays to a CSV file

np.savetxt(arcpy.env.workspace + "\Dados\Exportados\" + arquivo_com_as_medias + ".csv", lista_final, fmt="%10s %10s %10s", delimiter=";", header = "")

Improve this question
7
  • Your code sample is a bit long and the stated picture of the current output doesn't appear to exist. What are you currently using to write to a text file, is it on the line np.savetxt? on a quick skim I can't see an open statement. Have you tried the format command? This post might help stackoverflow.com/questions/35352045/…
    – Michael Stimson
    Commented Dec 19, 2018 at 1:20
  • I will try to upload again the picture. I will highlight the lines that I'm using to save the data to the CSV file. And I don't know this command that you said, I will read the page you suggested.
    – Elias
    Commented Dec 19, 2018 at 1:25
  • In the loop for mascara in bacias: you set arquivo_com_as_medias = "medias 01" so your output will forever overwrite medias 01.csv, surely you don't want to set the output name to a constant in the loop. How do you want the output file to be named? medias 01.csv medias 02.csv etc..? You also state your delimiter as a semicolon which isn't CSV (comma separated values), the delimiter should be a comma.
    – Michael Stimson
    Commented Dec 19, 2018 at 1:30
  • What I was trying to do is to keep saving in the same file, but not overwriting. It was supposed to add a new row every time and save the data. Because the time series that I'm using is a 54 years record, and I am making daily interpolation. So, if I have one file to each day, it would be hard to handle it later.
    – Elias
    Commented Dec 19, 2018 at 1:33
  • 1
    There isn't an option I can find to append with np.saveastxt, have a read of stackoverflow.com/questions/30376581/… which discusses this very topic. Personally I'd take control of the file I/O and use the standard with open(your_file,'a') as OutCSV: then iterate the rows, writing each one with OutCSV.write('{},{},{}\n'.format()), I've found this to be the best way to ensure the output is exactly what you want though it borders on micro-managing.
    – Michael Stimson
    Commented Dec 19, 2018 at 1:44

1 Answer 1

Reset to default
0

The solution that I found was to use the CSV File Reading and Writing library.

I just added the library in the beginning "import csv" and then added the following lines along the code.

First, in the begging, I needed to set the file directory to save the data.

#Opening the file to save the data
write_path = main_directory + "\\Dados\\Exportados\\" + "file_name" + ".csv"
file = open(write_path, 'a')
writer = csv.writer(file)

The 'a' parameter allows you to append to the end of the file

Then, in the step that I wanted the data to be saved, I added the following line:

#Writing to the file
writer.writerow([mes_pelo_nome, ano_pelo_nome, dia_pelo_nome, mascara, media])

In the end, I added the following line to close the file:

file.close()
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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