0

I am new to GIS processing, and I am using Python 3.6 to do so.
My current task is to read in several files of geographic data (specifically tiles of elevation data from JAXA), combine them for my area of interest, and then analyze them using Python.

Unfortunately, the data for each tile from AW3D30 come as a folder of six files: three tif files and two txt files. I haven't found any documentation about what each of these files contains, so I don't know which one(s) I need for the elevation data (although it must be at least the DSM file). And I don't know how to properly read them in and process them together, connect the proper lon/lat data to the pixels, etc.

For example, I read the *_AVE_DSM.tif file into Python using gdal, found out there is one band, and the band has a min and max value of None. 

altitudeData1 = gdal.Open('N034E138/N034E138_AVE_DSM.tif')
print("[ RASTER BAND COUNT ]: ", altitudeData1.RasterCount)       #-> 1
print("[ MIN ] = ", altitudeData1.GetRasterBand(1).GetMinimum())  #-> None
print("[ MAX ] = ", altitudeData1.GetRasterBand(1).GetMaximum())  #-> None

So I'm clearly missing something important.

My final desired output is geoPandas dataframe with a geometry column capturing the polygon of each pixel (i.e., not the center point, but the square's geometry) and a column for elevation. I am open to various workflows and intermediary steps.

Aaron Bramson
  • 925
  • 1
  • 9
  • 17
  • Be aware that doing a raster to vector conversion is super inefficient and memory intense. So try to avoid geopandas. What material are you following? Do you want gdal or is rasterio an option? Do you want to debug your code? What are you missing from readily available documentation on the net and answers on this site. – bugmenot123 May 29 '19 at 11:33
  • 1
    See also https://gis.stackexchange.com/questions/187877/how-to-polygonize-raster-to-shapely-polygons for how to polygonize raster data, if that is what you need. – joris May 29 '19 at 11:52
  • I can't imagine how this is off-topic, although perhaps it is a duplicate because the question/answer that @joris pointed to contained the needed bits of code (after updating them to current versions). I updated my question to reflect this solution for future noobs to find. – Aaron Bramson May 30 '19 at 03:50
  • 1
    All coding questions which do not contain code are likely to placed on hold. The off-topicness is an artifact of the way closure topics are organized by StackExchange. Please do not place answers in questions. If this is a duplicate then letting the question shift to Closed status is probably the best course.. – Vince May 30 '19 at 13:27
  • The question is not actually a duplicate because mine is about a specific dataset and a different desired output, but the solution uses (the newer version of) the same library. I added some worthless (i.e., not helpful to understand/answer the question) code to my question to show "what I tried so far" and moved my tentative answer to be an answer. – Aaron Bramson May 31 '19 at 03:02

2 Answers2

0

Based on @joris' comment, I updated the solution there to current versions:

thisFile = "yourFileRoot"
with rasterio.Env():
with rasterio.open(thisFile+'_AVE_DSM.tif') as src:
    with rasterio.open(thisFile+'_AVE_MSK.tif') as msk:        
        image = src.read(1)  ## first and only band      
        mask = msk.read(1)   ## first and only band
        mask[mask >= 4] = 0  ## Code 0 is valid, set other valid sources to 0
        mask[mask >= 1] = 1  ## Convert invalid, landwater, and sea locations to 1  
        mask = 1 - mask      ## Invert binary because this package excludes 0s in mask
        results = ({'properties': {'altitude': v}, 'geometry': s} for i, (s, v) in 
             enumerate(shapes(image, mask=mask, transform=src.transform)))

thisGeoData  = gp.GeoDataFrame.from_features(list(results))

Although this doesn't use 4 of the 6 files provided by Jaxa, it seems they do not include information necessary to get the altitude data.

The mask used by Jaxa and the mask required by rasterio are quite different, so this solution also includes mask edits to include only valid land data (from any source).

The GeoDataFrame allows convenient Pandas-style processing in Python as well as many output format options in addition to the geographic processing such as intersection, dissolve, etc.

Aaron Bramson
  • 925
  • 1
  • 9
  • 17
0

Before starting such a project, start by finding the metadata. That will have all the information on the model.

https://www.eorc.jaxa.jp/ALOS/en/aw3d30/aw3d30v22_product_e_a.pdf

You can use QGIS to figure out what is going on with the files.

Processing a raster using a geopandas dataframe is a very bad idea. You can used zonal or point statistics with raster data much more efficiently than by trying to convert this to a vector format. Zones may be rasters or vector file. See

https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#calculate-zonal-statistics

For the min max you need to compute the data. See https://gdal.org/tutorials/raster_api_tut.html

band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))
min = band.GetMinimum()
max = band.GetMaximum()
if not min or not max:
    (min,max) = band.ComputeRasterMinMax(True)
print("Min={:.3f}, Max={:.3f}".format(min,max))

Sorry, can’t get the python formatting correctly!!!

Nicolas Cadieux
  • 188
  • 7