How to open Willmott and Matsuura climate data on ArcMap?

1172
17
Jump to solution
07-18-2017 08:57 PM
Highlighted
New Contributor II

  Hello folks,

I'm in a kind of bad situation and just couldn't find any solution. 

So I'm trying to open the "Willmott, C. J. and K. Matsuura (2001) Terrestrial Air Temperature and Precipitation: Monthly and Annual Time Series (1950 - 1999)" dataset, which is avaiable in this link, but I just can't do it on ArcMap. I downloaded 7Zip to un-tar the file (it comes as .tar), but even after uncompressed, each year comes in a .1996, .1997, .1998 and so on format. This leaves me no clue on how to open them in ArcMap.

I am annexing the files so you can take a look.

I really found no way to open them online, I even tried just renaming to a .nc file and them opening with one of the tools in ArcMap, but it just didn't work. Also, I tried to rename them to a .txt format and then just add them to ArcMap. Also didn't work - it was added as an attribute table but it crashed upon opening.

If anyone knows how I can import those files to ArcMap, I'd be very thankful.

Cheers,
Joao

1 Solution

Accepted Solutions
Highlighted
Esri Esteemed Contributor

The files are ASCII, but are not in ASCII raster. They are point and have a spacing of 0.5 degrees. You can use a little scripts to import the files and convert to raster afterwards:

Below the script that I used. In my case I copied the files in a folder called "C:\GeoNet\Willmot" and created an empty File Geodatabase in that same folder called "test.gdb". The script will not only create the points with the 12 months as attributes, but also calculate total, mean, min and max values.

def main():
    import arcpy
    import os

    sr = arcpy.SpatialReference(4326)
    arcpy.env.overwriteOutput = True

    # input
    years = [1996, 1997, 2013, 2014]
    for year in years:
        print 'year:', year
        precip = r'C:\GeoNet\Willmot\precip.{0}'.format(year)

        # create output fc
        print "create output fc..."
        fc = r'C:\GeoNet\Willmot\test.gdb\precip_{0}'.format(year)
        ws, fc_name = os.path.split(fc)
        arcpy.CreateFeatureclass_management(ws, fc_name, 'POINT', None, None, None, sr)

        print "add fields..."
        months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August',
                  'September', 'October', 'November', 'December', 'Total', 'Mean', 'Min', 'Max']
        for month in months:
            arcpy.AddField_management(fc, month, 'DOUBLE')

        # insert cursor
        print "read file and add features..."
        cnt = 0
        flds = ['SHAPE@']
        flds.extend(months)
        with arcpy.da.InsertCursor(fc, flds) as curs:
            with open(precip, 'r') as f:
                for line in f.readlines():
                    try:
                        cnt += 1
                        if cnt % 1000 == 0:
                            print " - processing point", cnt
                        while '  ' in line:
                            line = line.replace('  ', ' ')
                        line = line.strip()
                        lst = line.split(' ')
                        lst = [float(a) for a in lst]
                        pntg = arcpy.PointGeometry(arcpy.Point(lst[0], lst[1]), sr)
                        lst_row = [pntg]
                        data = lst[2:]
                        stats = [sum(data), sum(data) / 12.0, min(data), max(data)]
                        lst_row.extend(data)
                        lst_row.extend(stats)
                        row = tuple(lst_row)
                        curs.insertRow(row)

                    except Exception as e:
                        print line
                        print cnt, e

##    arcpy.PointToRaster_conversion(in_features="precip_2013",
##                value_field="January",
##                out_rasterdataset="C:/GeoNet/Willmot/test.gdb/precip_2013_jan",
##                cell_assignment="MEAN", priority_field="NONE", cellsize="0,5")

if __name__ == '__main__':
    main()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

To convert the point featureclasses to rasters you can use to tool Point to Raster—Conversion toolbox | ArcGIS Desktop. You will have to set the extent to X from -180 to 180 and Y from -90 to 90 and the pixel size to 0.5 (degrees). I have attached the FGDB to the thread, which contains some raster that you can use as snap raster and to set the size of the pixels and the extent.

Below a screenshot of the resulting raster and points on top of them: 

View solution in original post

17 Replies
Highlighted
MVP Esteemed Contributor

From the description and so I don't have to load a bunch of zips... is there any chance they are ASCII data files?

If so try Ascii to Raster to do the conversion

Reply
0 Kudos
Highlighted
New Contributor II

Hey Dan, how you doing? Thanks again for reaching out to help me. 

On the official website it says they're ASCII data files. When I use Ascii to Raster, however, ArcMap reports an instant failure. Even when I renamed the .1996 (and so on) files to .txt they still report the same error. Curiously though, the .txt file can be open in Notepad and displays values that seem proper for ArcGIS. I'm not home right now so I can't report the exact error reported by Ascii to Raster; I'll post an update with the error term later.

Thanks again and have a good day

Reply
0 Kudos
Highlighted
MVP Esteemed Contributor

if it opens in a text file, copy the first couple of lines.  unfortunately the term 'ascii' file is bantered around and sometimes doesn't mean the same thing.  You can compare what you have with the ascii file description in the help topic to see if it begins with a 6 row header of information on cell size and location followed by the data, or whether it is a stream of coordinates.

Reply
0 Kudos
Highlighted
New Contributor II

Ok Dan I'm home right now so I can provide you with better answers.

Beginning with the text file, I saw the help file in ArcMap help and it doesn't seem like it begins with this header of information. Actually, the lines are coordinate points (I've pasted the first line on the bottom of this reply).

When I try to open this .txt with Ascii to Raster, I get errors:

010328: Syntax error at or near symbol SPACE
010267: Syntax error in parsing grid expression

When I try to open with Add Data, I get this:

Either way the "original" .1996, .1997 format and so on does not enable any add option. So I really think it must be a .txt file, just need to figure out how to get ArcMap to read it, right?


P.S.

First line: 
-179.750 71.250 5.7 7.9 9.6 7.2 14.0 15.6 26.0 21.5 12.9 13.9 20.3 5.0
-179.750 68.750 8.0 4.8 10.5 7.3 4.9 14.6 59.0 55.5 8.9 28.3 27.3 4.8
-179.750 68.250 12.2 8.7 16.9 7.5 8.9 20.7 60.4 54.9 9.9 28.1 32.4 10.8
-179.750 67.750 19.8 16.6 33.6 7.8 16.9 33.9 68.5 57.4 13.5 28.1 46.4 22.2
-179.750 67.250 29.6 27.2 62.1 7.1 28.2 52.4 79.4 58.7 21.4 30.5 73.5 38.0
-179.750 66.750 30.2 33.1 94.5 1.2 42.3 72.6 95.2 62.3 36.8 35.7 105.2 47.1
-179.750 66.250 25.1 33.5 108.3 0.0 53.1 86.2 114.9 71.7 49.3 39.9 116.6 48.0
-179.750 65.750 30.3 32.5 85.4 0.0 39.2 65.0 87.3 72.6 33.4 26.8 86.2 40.0
-179.750 65.250 43.0 32.2 67.3 3.0 29.0 44.6 62.4 76.0 19.8 20.7 65.5 38.2
-179.75

Reply
0 Kudos
Highlighted
MVP Esteemed Contributor

well those look like 

Longitude, Latitude, Jan, Feb, Mar.... Dec

You could put those into excel as a space delimited text file, then add the headers... then add it to arcmap as an XY layer

http://desktop.arcgis.com/en/arcmap/latest/map/working-with-layers/adding-x-y-coordinate-data-as-a-l...

create a feature class out if the event layer, and define its coordinate system as a GCS WGS 1984, since the data appear to be in decimal degrees.

To get that into a raster format... don't know offhand, but worse comes to worse, you could create a fishnet so that the polygons of the fishnet have those points as their centroid and do a spatial join of the points to the fishnet.... I am sure there is a better way, but it is late

Reply
0 Kudos
Highlighted
New Contributor II

Wow Dan that was probably the solution... I'm trying it out but could you aid me with 2 more things? First, I know obviously the number of cells and columns and the nodata_value; but I'm not pretty sure I understood the other 3 criteria: XLLCENTER/XLLCORNER, YLLCENTER/YLLCORNER and CELLSIZE.

I read that in the XLL/YLL cell I must input the southwesternmost cell in my file, which would be -179,250/-89,75. 
The cellsize area, according to the README file, is 0,5 degrees. So I inputed the following header:

NCOLS18
NROWS85800
XLLCORNER-179.75
YLLCORNER-89.75
CELLSIZE0.5
NODATA_VALUE-9999

But, when opening in ArcMap by adding XY, it won't let me choose the right X and Y columns (which would be columns 1 and 2, respectively). The option only starts at column 3. I tried to move all the cells two units to the right; but again it would only let me choose from column 5 onwards.

Any idea of why this is happening?

Second thing: I want to add some statistics (mean, max, min, and the coef. of variation). Should I add them in Excel already or can it be problematic? The final use of the data is: I have polygons representing each brazilian municipality and I want to translate the data from the raster to each municipality.

But anyway I think that once I figure that out it's gonna work. Thanks Dan, you've been wonderful.

Reply
0 Kudos
Highlighted
MVP Esteemed Contributor

You are mixing the ascii format with a tabular data format.  the header format that you have reproduced is for a different type of 'ascii' data.  If you get that original file into excel, then do as I suggested.  Get it mapped correctly first before you try to join the data into another format.  

I would also go back to the site where you got the data... surely there is some documentation somewhere as to the exact format and how to process into a useable form

Reply
0 Kudos
Highlighted
New Contributor II

Wait now I'm lost again.

I got the precip.1996 file, renamed it to .txt, them opened it with Excel. I choose a delimited-space option and then it opened as a regular spreadsheet. I added the header in the Excel cells and then saved the file again in a .txt format. I see now that it saved a a tab-delimited text data whereas I wish it was space-delimited.  I used this Esri ASCII raster format—Help | ArcGIS for Desktop  as the example for the ASCII header. 

So what should I exactly do and how do I get an example of the ASCII header? Sorry for taking your time...

Edit: I just tried two things:

a) opening the file in excel, no header, just saved it as a .csv file. Them I could add as a X, Y layer but all the data were points. I exported the data inside ArcMap but as they were points, I think it is kind of useless, right?
b) opening the file in excel, adding the header, saved it as a .csv file. Opened word, replaced all ";" by spaces " ", and saved it as a .txt. So now it was a space-delimited text file with a header according to ArcGis specs. I tried using ASCII to Raster with this space-delimited text file but got the error 010328 (syntax error at or near symbol SPACE) and 010267 (syntax error in parsing grid expression).


And only now I realized the fishnet thing you were talking about. Got to take a look into that. So if I get points like I did in a), I can make them cells, that's it right? I just wonder why inputting the header like I did in b) didn't work...

*As for the website, all I get in the README files is a technical description of the data and this is in the header of the downloads page: "Important note: all files available here for downloading have been "compressed" under Unix in order to save space. A single compressed file has ".Z" extension. You may need to respecify the ".Z" extension when you save a compressed file. Some of our archives also have multiple files (such as time series) that are "tar"ed under Unix and have a ".tar" extension. As a result, files must be "untared" and/or "uncompressed" to return them to their original ASCII format, before they may be used. 

How to cite our data: we recommend that users cite pertinent information contained in the README files associated with the data sets. A citation should contain the creators' (our) names and the creation and/or most-recent-modification date, as well as the title and URL associated with the data set. An example of how to cite one of our data sets is:

Willmott, C. J. and K. Matsuura (2001) Terrestrial Air Temperature and Precipitation: Monthly and Annual Time Series (1950 - 1999), http://climate.geog.udel.edu/~climate/html_pages/README.ghcn_ts2.html. 

Please change the content accordingly."

Reply
0 Kudos
Highlighted
Esri Esteemed Contributor

The files are ASCII, but are not in ASCII raster. They are point and have a spacing of 0.5 degrees. You can use a little scripts to import the files and convert to raster afterwards:

Below the script that I used. In my case I copied the files in a folder called "C:\GeoNet\Willmot" and created an empty File Geodatabase in that same folder called "test.gdb". The script will not only create the points with the 12 months as attributes, but also calculate total, mean, min and max values.

def main():
    import arcpy
    import os

    sr = arcpy.SpatialReference(4326)
    arcpy.env.overwriteOutput = True

    # input
    years = [1996, 1997, 2013, 2014]
    for year in years:
        print 'year:', year
        precip = r'C:\GeoNet\Willmot\precip.{0}'.format(year)

        # create output fc
        print "create output fc..."
        fc = r'C:\GeoNet\Willmot\test.gdb\precip_{0}'.format(year)
        ws, fc_name = os.path.split(fc)
        arcpy.CreateFeatureclass_management(ws, fc_name, 'POINT', None, None, None, sr)

        print "add fields..."
        months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August',
                  'September', 'October', 'November', 'December', 'Total', 'Mean', 'Min', 'Max']
        for month in months:
            arcpy.AddField_management(fc, month, 'DOUBLE')

        # insert cursor
        print "read file and add features..."
        cnt = 0
        flds = ['SHAPE@']
        flds.extend(months)
        with arcpy.da.InsertCursor(fc, flds) as curs:
            with open(precip, 'r') as f:
                for line in f.readlines():
                    try:
                        cnt += 1
                        if cnt % 1000 == 0:
                            print " - processing point", cnt
                        while '  ' in line:
                            line = line.replace('  ', ' ')
                        line = line.strip()
                        lst = line.split(' ')
                        lst = [float(a) for a in lst]
                        pntg = arcpy.PointGeometry(arcpy.Point(lst[0], lst[1]), sr)
                        lst_row = [pntg]
                        data = lst[2:]
                        stats = [sum(data), sum(data) / 12.0, min(data), max(data)]
                        lst_row.extend(data)
                        lst_row.extend(stats)
                        row = tuple(lst_row)
                        curs.insertRow(row)

                    except Exception as e:
                        print line
                        print cnt, e

##    arcpy.PointToRaster_conversion(in_features="precip_2013",
##                value_field="January",
##                out_rasterdataset="C:/GeoNet/Willmot/test.gdb/precip_2013_jan",
##                cell_assignment="MEAN", priority_field="NONE", cellsize="0,5")

if __name__ == '__main__':
    main()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

To convert the point featureclasses to rasters you can use to tool Point to Raster—Conversion toolbox | ArcGIS Desktop. You will have to set the extent to X from -180 to 180 and Y from -90 to 90 and the pixel size to 0.5 (degrees). I have attached the FGDB to the thread, which contains some raster that you can use as snap raster and to set the size of the pixels and the extent.

Below a screenshot of the resulting raster and points on top of them: 

View solution in original post