geodesymiami / insarmaps

3 stars 0 forks source link

intermediate pixel sizes for display #55

Closed falkamelung closed 1 year ago

falkamelung commented 2 years ago

We have the "actual pixel size" option, but it is not quite right. When I select this option there are white space in between although there should not.

It would be good to have some more options between the actual pixel size and the default size. E.g. a button to half or double the pixel size and/or to give the radius in meters, as well as smaller than actual pixel size.

It also will be good to display the pixel dimension (in x-direction) somewhere.

It would be good to also have the ability to switch from circle to rectangle so that the whole area is covered when the actual pixel size is selected.

We never talked about map projectjections. Are we using Equirectangular? If so a pixel should show up as a circle and a square pixel (we have lat-on squares) should show up as a square. MapBox has options for that (see links)

https://www.mapbox.com/blog/adaptive-projections https://en.wikipedia.org/wiki/Equirectangular_projection#/media/File:Plate_Carrée_with_Tissot's_Indicatrices_of_Distortion.svg

https://insarmaps.miami.edu/start/25.8786/-80.1209/16.8865?flyToDatasetCenter=false&startDataset=S1_IW3_048_0081_0082_20160927_20210703_PS040

image image

XXXXXXXXXXXXXXXXXXXX

image image

https://insarmaps.miami.edu/start/-7.5369/110.4420/15.6673?flyToDatasetCenter=false&startDataset=S1_IW2_076_0618_0621_20180110_XXXXXXXX_S07624_S07465_E110371_E110513&pointLat=-7.54070&pointLon=110.44117&minScale=-10&maxScale=10&startDate=20180110&endDate=20220113

stackTom commented 2 years ago

We have the "actual pixel size" option, but it is not quite right. When I select this option there are white space in between although there should not.

It would be good to have some more options between the actual pixel size and the default size. E.g. a button to half or double the pixel size and/or to give the radius in meters, as well as smaller than actual pixel size.

It also will be good to display the pixel dimension (in x-direction) somewhere.

It would be good to also have the ability to switch from circle to rectangle so that the whole area is covered when the actual pixel size is selected.

We never talked about map projectjections. Are we using Equirectangular? If so a pixel should show up as a circle and a square pixel (we have lat-on squares) should show up as a square. MapBox has options for that (see links)

https://www.mapbox.com/blog/adaptive-projections

https://en.wikipedia.org/wiki/Equirectangular_projection#/media/File:Plate_Carrée_with_Tissot's_Indicatrices_of_Distortion.svg

https://insarmaps.miami.edu/start/25.8786/-80.1209/16.8865?flyToDatasetCenter=false&startDataset=S1_IW3_048_0081_0082_20160927_20210703_PS040

image image

XXXXXXXXXXXXXXXXXXXX

image image

https://insarmaps.miami.edu/start/-7.5369/110.4420/15.6673?flyToDatasetCenter=false&startDataset=S1_IW2_076_0618_0621_20180110_XXXXXXXX_S07624_S07465_E110371_E110513&pointLat=-7.54070&pointLon=110.44117&minScale=-10&maxScale=10&startDate=20180110&endDate=20220113

I believe mapbox uses Mercator projection, not that one you mention. In regards to the squares, I haven't looked at mapbox in a long time, but last checked (over a year ago) you couldn't do squares with mbtiles. If it's not supported, we could add another layer on top of squares. The problem is that adding custom layers is slow and memory intensive. Would only be feasible when user is zoomed in a lot and not many points are displayed.

Also, what is the size of each pixel in the h5 file. Is this constant among files, or is there a variable in the file that tells us the size? (I vaguely remember x step and y step. Is it that?)

falkamelung commented 2 years ago

I checked Google maps. It uses transverse Mercator (UTM) which we probably use as well. I don't know much about map projections but I think Pixel with identical Lat/Lon spacing should show up in Equirectangular as a square and in UTM as a rectangle. But if mapbox can't do squares or rectangles that is not relevant.

In UTM X and Y are equidistant so "actual size" can't work. We would need "Actual size-long" and "Actual size-lat". In this example we have a slight overlap in Long direction but too wide spacing in Lat direction. A complication is that we can have different spacing in lat and long direction, but if you have this information that should not matter.

image

https://insarmaps.miami.edu/start/25.7814/-80.1900/19.5558?flyToDatasetCenter=false&startDataset=S1_IW3_048_0080_0083_20160927_20211112&minScale=1.5&maxScale=-1.5

So, whenever you have time, it would be goof to have these two options but at the same time also have the option to change the pixel size, either to give it in meter or to double/half it with a button. And display the pixel size in meters somewhere.

stackTom commented 2 years ago

Question, I am looking at the code. I see we have x_step and y_step. Are the units of these in meters?

stackTom commented 2 years ago

Also, I can do actual lat and actual lon. However, I don't know if it will affect the results much? I also don't know much about map projections, but I imaging the fundamental problem is how to display a spherical earth onto a 2d screen. At the high zoom levels this feature is bound to be used in, is there really a big difference among the different projections and between lat and long?

falkamelung commented 2 years ago

What is the value of those attributes? Is it of the order of what s shown below? That is degree. It should be easy to roughly convert to meters (depends on latitude).

//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/mintpy[1007] info.py S1_IW2_076_0618_0621_20180110_XXXXXXXX_S08141_S06941_E109662_E110704.he5 | grep STEP
  X_STEP                                      0.0008411197140170134
  Y_STEP                                      -0.0006718467894873353
//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/mintpy[1008] cd ../minopy/single_reference_network/

//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/minopy/single_reference_network[1010] info.py S1_IW2_076_0618_0621_20180110_XXXXXXXX_S07624_S07465_E110371_E110513.he5 | grep STEP
  X_STEP                                         4.1558472324803656e-05
  Y_STEP                                         -0.00015023347896117587

Do you mean you can display a marker given delta_lat and delta_long? I came to the conclusion that the map projection does not matter.

stackTom commented 2 years ago

What is the value of those attributes? Is it of the order of what s shown below? That is degree. It should be easy to roughly convert to meters (depends on latitude).


//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/mintpy[1007] info.py S1_IW2_076_0618_0621_20180110_XXXXXXXX_S08141_S06941_E109662_E110704.he5 | grep STEP

  X_STEP                                      0.0008411197140170134

  Y_STEP                                      -0.0006718467894873353

//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/mintpy[1008] cd ../minopy/single_reference_network/

//tg851601-api-u-1/data/HDF5EOS/MerapiSenDT76/minopy/single_reference_network[1010] info.py S1_IW2_076_0618_0621_20180110_XXXXXXXX_S07624_S07465_E110371_E110513.he5 | grep STEP

  X_STEP                                         4.1558472324803656e-05

  Y_STEP                                         -0.00015023347896117587

Do you mean you can display a marker given delta_lat and delta_long? I came to the conclusion that the map projection does not matter.

Yes. The code assumes it's degrees. I was confused but it makes sense since you say the units are degrees.

No, the map shows a circle given a radius. What I meant is if you zoom in enough to a sphere, it's basically a planar surface.

falkamelung commented 2 years ago

So, you can do only circles and not rectangles, can you?

The third answer in this thread gives the formula on how to convert lat/long to meters. This is good enough for now.

https://stackoverflow.com/questions/639695/how-to-convert-latitude-or-longitude-to-meters/39540339#39540339

stackTom commented 2 years ago

Mapbox requires us to give the circle size in pixels. So, my approach is to get the amount of degrees per pixel at the current zoom level. Then, we can get the amount of pixels by dividing DEGREES / degrees per pixel.

However, how do we calculate DEGREES? The code currently gets the larger of x_step and y_step (let's call it largerStep). It then does: DEGREES = SquareRoot(largerStep^2 + largerStep^2). I don't understand why we do this. Do you know? I looked at the git commit comments for this, and it says "calculate per dr falk's formula.", so maybe it's something you suggested that I don't understand.

I believe DEGREES should just be x_step or y_step, no? Is the formula above a way to account for the fact that we use circles instead of rectangles?

stackTom commented 2 years ago

I've analyzed everything more. There seems to be a new way to draw rectangles. However, it would require one of different modifications: 1) modify the ingest scripts. right now, these ingest the data as individual points. Would need to modify it so that the data is ingested as polygons (starting coordinates with x_step width and y_step height). Not sure the repercussions of this in terms of the performance of the website. Might slow it down, might not. 2) I can add a new custom set of polygons to the site. It would work exactly like on the fly json works, where I take the currently rendered points, and I add a customized layer (not dependent on mbtiles).

I think option 2 is the better solution in terms of simplicity. Although it might be slower than option 1. What do you think?

falkamelung commented 2 years ago

To answer the first question, this formula does not make much sense to me. The X_STEP and Y_STEP are the DEGREES?

stackTom commented 2 years ago

To answer the first question, this formula does not make much sense to me. The X_STEP and Y_STEP are the DEGREES?

I completely agree. I have no clue where I got this formula from. Maybe we decided on it at some point to account for the fact that we are dealing with circles and not triangles.

falkamelung commented 2 years ago

I just send you an email. We better talk as I am not sure I understand.

If we can have a polygon for each point and then fill it, that would be fantastic.

falkamelung commented 2 years ago

So there are three issues with the high-res Miami data

  1. The dot/rectangle size
  2. Sara: The problem is that during the conversion of data to insarmaps format, it transfers them to a regular grid with a specific spacing which is greater than the original file, so you loose the resolution.
  3. occasional ingest error. I will be happy to cut a data set into chunks, but we have not figured out whether the file size, number of pixels, or number of dates is the problem.

Miami data on jetstream (PS file has less dots)

ssh -YC -o ServerAliveInterval=60  centos@129.114.104.223

/data/HDF5EOS/Miami2SenAT48/minopy_SM
drwxrwxr-x. 2 centos centos      4096 Apr 20 02:32 JSON
-rw-rw-r--. 1 centos centos 364802200 Apr 19 18:20 S1_IW3_048_0081_0082_20150921_20211112.he5
-rw-rw-r--. 1 centos centos 364742393 Apr 19 18:27 S1_IW3_048_0081_0082_20150921_20211112_PS.he5

https://insarmaps.miami.edu/start/25.8774/-80.1215/17.1694?flyToDatasetCenter=false&startDataset=S1_IW3_048_0081_0082_20150921_20211112_PS&minScale=0.8&maxScale=-0.8

Google Earth Insarmaps
image image
stackTom commented 2 years ago

1) Im currently workin on the dot/rectangle size. 2) I'm not sure I understand. I thought the he5 file was a regular grid with a specific spacing between each point. Also, what does Sara mean? 3) Will need to investigate

stackTom commented 2 years ago

Regarding 2). I just saw your email - yes, the original ingest script takes the X_FIRST and Y_FIRST points and then just iterates with X_STEP and Y_STEP. If a point has nan values, it is not displayed. What is the format of this irregular data h5 file? We can do it, I would just need to know the format.

falkamelung commented 2 years ago

Sorry. I am online now. Sara is my student who tried to ingest.

The high-resolution data don't have equal spacing, in particular not in x-direction. Does insarmaps need equal spacing? If we could support unequal spacing that would be perfect. The actual lat long of each pixel is somewhere in the S1 file.

Assuming that this is not possible, the way around would be to manually specify the spacing to be used so that (almost) all points are captured. That could be lots of NaNs, though. As files could get large we have to figure out the reasons for the occasional ingest errors. If a file gets to big it should not allow you to create it and/or ingest.

Other people use Google Earth or Arcview. this works well for irregular speccing but we don't have the recoloring option which is the big strength of insarmaps.

stackTom commented 2 years ago

Yes, irregular spacing shouldnt be an issue. Its just a matter of me learning the format of the irregular h5 file and modifying the ingest script appropriately. Shouldnt be difficult. I am currently not home for the rest of the day (hopefully the server is back when I am home). Do you have a document for the format for irregular spacing, so I can start thinking about how to do it?

falkamelung commented 2 years ago

That would be fantastic. The S1 file has fields/arrays with the lat and long information., I think. If not I have to look. There is a google earth function in MintPy that does this. Save_timeseries_kmz or similar. Let me know if it is not obvious and I will look for the function.

falkamelung commented 2 years ago

I checked it. It is right in the S1* file:


HDF5 dataset "/HDFEOS/GRIDS/timeseries/geometry/latitude              ": shape (565, 1364)         , dtype <float32>
  MissingValue    0.0
  Title           latitude
  Units           degrees
  _FillValue      0.0

HDF5 dataset "/HDFEOS/GRIDS/timeseries/geometry/longitude             ": shape (565, 1364)         , dtype <float32>
  MissingValue    0.0
  Title           longitude
  Units           degrees
  _FillValue      0.0

The MintPy command does not use the S1* file but I guess that was just not implemented:

save_kmz_timeseries.py --vel geo/geo_velocity.h5 --tcoh geo/geo_temporalCoherence.h5 --mask geo/geo_maskTempCoh.h5 geo/geo_timeseries_demErr.h5
stackTom commented 2 years ago

Sounds good, I will modify the ingest script. Ill make a new issue.

Regarding the actual size button - do you think it would be worth it to ingest all data as polygons instead of rectangles like I mentioned previously?

1) I don't know how performance will be. 2) It would require us to re ingest all data sets.

I could ingest a couple of test data set and we can try them on my development server. What do you think?

Or should we just keep the ingest the same but have a separate "actual size" mode? Ive already begun work on this. I thought this would be easier but its actually just as much if not more work than if we pursued the polygon method...

stackTom commented 2 years ago

I think it will be good to start from scratch and re ingest everything anyways (even if we don't pursue the polygon method). I am now realizing there is a "get_lat_long" function in mintpy. I think I should use this to get lat and long of points for ingest. It seems to handle both grid mode (currently supported by insar ingest scripts) and lat/long mode (the mode I need to add to support the surfside building data). Grid mode seems similar to the code I use in the insar ingest scripts, but I see one minor difference.

Will make sense to re-ingest everything, and use this function so that everything aligns perfectly on the site and mintpy... Let me know what you think.

falkamelung commented 2 years ago

So we can’t do both ? Keeping the old way and add an option for using exact lat/long? If not and/or if there are performance issues we may want to do two servers. Keep this old one and have another one for high resolution.

There is actually a new jetstream cloud which is much better then the old one. The old one will be discontinued and of this month so I have to switch anyway. Just have to figure out ow to do it:

https://jetstream-cloud.org

stackTom commented 2 years ago

We can do exact/lat long right now. No need to change anything.

We'd need to change if we want to try uploading datasets as polygons (squares) instead of circles.

stackTom commented 2 years ago

Don't think we can ingest as polygons. Please see here: http://129.213.120.104/start/19.4777/-155.6024/17.5354?flyToDatasetCenter=false&startDataset=S1_IW12_087_0527_0531_20180908_XXXXXXXX_N19390_N19549_W155679_W155507&minScale=3&maxScale=-3

When we ingest as points (as we currently do), mapbox is able to only render some points, and it then interpolates/resizes the points so the image looks good from a zoomed out view. If I render polygons, it seems to be unable to do this. So, although the pixels are correct when you are zoomed in, there appear empty spaces between them when you are zoomed out (as not all of them are rendered).

So, we need to ingest as points, and then I need to add custom code for the actual pixel size button as I was previously doing.

falkamelung commented 2 years ago

At high resolution is looks good to me? (but the timeseries display does not work. I only see problems (white areas) when you zoom out but for the Miami data people would always look at high resolution???

image image

I would be OK to add a suffix to the filename for those I need in high resolution (Miami only at the moment), and ingest both, the regular way and with rectangles.

stackTom commented 2 years ago

Yes, there is empty space when zoomed out. The timeseries doesn't work as that was just a quick hack to get rectangles working, but they will work.

There are many ways I can tackle/design this high resolution/actual size button problem.

Each has their own pros/cons.

I basically need to know - does every dataset need to have an actual pixel size button? If not, we can ingest two versions for each dataset - low resolution (circles) and actual size/high resolution (rectangles). And then I can switch between the two depending on zoom level/user clicking a button. Benefits is less code, possibly better performance. Drawback is more space on our server (two mbtiles per area).

If however, every dataset needs high res/actual size functionality, it's impractical to upload two mbtiles per area... I will need to manually create "on the fly" rectangles in the browser. Benefits is no need for two mbtiles. Drawback is more code and complexity, and possibly slightly slower performance (might need to reload when panning for example).

Not sure which is best at this stage. Which do you think?

falkamelung commented 2 years ago

No we don't need it for all datasets. But I am a bit hesitant to eliminate the current "actual pixel size" functionality before we have verified that the Miami data work well. Can't we add a new option to display the rectangles if available?

stackTom commented 1 year ago

Please try this again. I have added support for it.

stackTom commented 1 year ago

Completed.