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.

  • 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 gis.stackexchange.com/questions/187877/… 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 3: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 3:02

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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