brycefrank / pyfor

Tools for analyzing aerial point clouds of forest data.
MIT License
93 stars 19 forks source link

Flipping the canopy height model? #26

Closed bw4sz closed 5 years ago

bw4sz commented 5 years ago

Hey Bryce,

A tiny question today, i'll really try not to lean too much on ya.

I'm working on deep learning on joint RGB and lidar imagery, so I need to crop and overlay both data at the same window. I'm all good to go, but it looks like the canopy height model is mirrored over the y axis when plotting.

When overlaying the RGB (image) and LIDAR (CHM), its easy to see.

        fig, ax = pyplot.subplots()
        ax.imshow(image[:,:,::-1])
        ax.matshow(CHM.array,alpha=0.4)
        pyplot.show()
screen shot 2018-11-20 at 9 32 34 am

Then flip the CHM over the y axis

        fig, ax = pyplot.subplots()
        ax.imshow(np.flip(image[:,:,::-1],0))
        ax.matshow(CHM.array,alpha=0.4)
        pyplot.show()
screen shot 2018-11-20 at 9 28 06 am

I think you already noticed this in pyfor:

https://github.com/brycefrank/pyfor/blob/eb158be09b197a8a1f7c1d9f3e00dced79936a25/pyfor/rasterizer.py#L170

Do you know what's happening here? Is it

1) That matplotlib "shows" the y axis inverted. So you are fixing it visually (this might be dangerous?). 2) The matrix itself is flipped, because of the way pyfor calculates the CHM, so you invert it for visualization 3) Neither 1 or 2, meaning my image is flipped (seems less likely).

It matters because the next step is to concatenate the CHM to the 3 channel rgb to create a 4 channel input to the deep net. I don't want to stick it on upside down.

brycefrank commented 5 years ago

Ben,

You have found a bit of a sore spot for me. I have been looking to get this consistent myself. I am wrapping up a small assignment right now, but I can hop on this in a few minutes.

bw4sz commented 5 years ago

Sure, i'm trying to think about how to make a reproducible example. its within the deep learning data generator that would take alot for you to bury into. Trying now. My labmate thinks its option 2. The rgb image plot shows the UTM N increasing as you go up, which i guess is arbitrary, but we should be consistent.

brycefrank commented 5 years ago

To be honest there is a deep nest of flips and unflips in pyfor (one of those things, you know). The least I can tell you is that line 170 only flips for visualization. It may be worthwhile to matshow the image and CHM separately in two different plot for debugging. You could probably take advantage of the top left corner of the image you have to make sure they are correct. You can access the chm array like so:

pc = Cloud('my_las.las')
my_chm = pc.chm(1)
chm.array

My other though is that when I generate the bins for the points to convert to a raster is that, I think, they start from the bottom left, rather than the top left. Hence all the flipping/flopping. I will come up with a more formal solution soon. Keep me posted on here (or gitter if you want more instant feedback)!

bw4sz commented 5 years ago

Okay this should be reproducible.

Data: Archive.zip

import pyfor
from PIL import Image
from matplotlib import pyplot
import numpy as np

pc=pyfor.cloud.Cloud("/Users/ben/Documents/DeepLidar/data/SJER/SJER_002.laz")
im=Image.open("/Users/ben/Documents/DeepLidar/data/SJER/SJER_002.tif")
#Maybe numpy inverts the image from PIL?
image = np.array(im)    

pc.normalize(2)
chm = pc.chm(cell_size = 0.1 , interp_method = "nearest", pit_filter = "median", kernel_size = 9)

fig, ax = pyplot.subplots()
ax.imshow(image)
ax.matshow(chm.array,alpha=0.4)
fig.title="Original"

#versus
fig, ax = pyplot.subplots()
ax.imshow(image)
ax.matshow(np.flip(chm.array,0),alpha=0.4)
fig.title="CHM Flipped"
pyplot.show()

figure_1 figure_2

I'm going to assume that the CHM is flipped for some reason (UTM N facing down, not up) and just move forward.

brycefrank commented 5 years ago

Ben,

I just arrived to the same conclusion. The index [0,0] of chm.array refers to the bottom left of the CHM (in real space). When we display this in matplotlib using plt.matshow(chm.array) it is displayed upside-down for this reason.

I will need to think about the implications down the line for flipping it up in rasterizer.Raster. It seems to be an easy fix, but I need to write some tests for this sort of thing anyway. I

I am committing the rest of the day (barring any other unforeseen work) to the testing suite and other fixes you have presented. I will let you know when I push to laz_fix, and I will include a UserWarning in this push when the array is constructed such that [0,0] refers to the top left when rasterizer.Raster is initialized.

bw4sz commented 5 years ago

Great. I think that makes sense, numpy arrays use the top left origin, and it probably is the most pythonic thing to do. I will say that the reverse is true in R, requiring some impromptu math on occasion. I'll make a quick fix and pull when ready.

On Tue, Nov 20, 2018 at 10:25 AM Bryce Frank notifications@github.com wrote:

Ben,

I just arrived to the same conclusion. The index [0,0] of chm.array refers to the bottom left of the CHM (in real space). When we display this in matplotlib using plt.matshow(chm.array) it is displayed upside-down for this reason.

I will need to think about the implications down the line for flipping it up in rasterizer.Raster. It seems to be an easy fix, but I need to write some tests for this sort of thing anyway. I

I am committing the rest of the day (barring any other unforeseen work) to the testing suite and other fixes you have presented. I will let you know when I push to laz_fix, and I will include a UserWarning in this push when the array is constructed such that [0,0] refers to the top left when rasterizer.Raster is initialized.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brycefrank/pyfor/issues/26#issuecomment-440380719, or mute the thread https://github.com/notifications/unsubscribe-auth/ABJwrPrUOd4kyEQUldpNY5UAHOF6jB2Gks5uxEkPgaJpZM4Yrm5J .

-- Ben Weinstein, Ph.D. Postdoctoral Fellow University of Florida http://benweinstein.weebly.com/

brycefrank commented 5 years ago

Ben,

A quick update. I am debating fixing this way up at the top (when the bins themselves are assigned). This will likely take a day or two to solidify. At least we know the temporary solution works for now. There are a lot of loose threads as well that I think would benefit from fixing all the way back in Grid.__init__ specifically:

https://github.com/brycefrank/pyfor/blob/laz_fix/pyfor/rasterizer.py#L28-L30

Edit

After spending a bit of time with it it is much more maintainable to do the fix all the way up in Grid. All of the calls to np.flipud are no longer needed. I have to un-do a few more of these types of calls, but I am much more happy with how it is turning out. I will push to laz_fix in an hour or so.

brycefrank commented 5 years ago

This fix (along with the laz_fix branch) has been merged to master. If you update, make sure to remove any flipping you have done yourself. It should no longer be needed (a UserWarning is issued any time a Grid object is initialized, I'll leave this in until 0.3.2 is merged to master, which will be a while).

bw4sz commented 5 years ago

great. pulling now

On Wed, Nov 21, 2018 at 11:53 AM Bryce Frank notifications@github.com wrote:

This fix (along with the laz_fix branch) has been merged to master. If you update, make sure to remove any flipping you have done yourself. It should no longer be needed (a UserWarning is issued any time a Grid object is initialized, I'll leave this in until 0.3.2 is merged to master, which will be a while).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/brycefrank/pyfor/issues/26#issuecomment-440790279, or mute the thread https://github.com/notifications/unsubscribe-auth/ABJwrKQX7zwC3x85uoMrMpjY4vLeoNPJks5uxa9JgaJpZM4Yrm5J .

-- Ben Weinstein, Ph.D. Postdoctoral Fellow University of Florida http://benweinstein.weebly.com/

brycefrank commented 5 years ago

Looking forward to the next challenge ;)