Changeset 24879

Timestamp:

05/20/20 18:12:57 (5 years ago)

Author:

Eric.Larour

Message:

CHG: implemented read_contours.m routine.

File:

• issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m (modified) (1 diff)

issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m

@@ -76 +76 @@
         end
         %}}}
+                function readcontours(self) % {{{
+                        %we are reading the contours from each shape file for each region, reproject
+                        %the latlong in these contours to the local laea projection, and rewrite the shapefile
+                        %we a different name extension.
+
+                        for i=1:self.nregions,
+                                
+                                disp(['reading shapefile for region: '  self.regions(i).name]);
+                                contours=shaperead([self.root '/' self.regions(i).name '.shp']);
+
+                                disp(['concatenating contours (X,Y)']);
+                                count=0;
+                                for j=1:length(contours),
+                                        count=count+length(contours(j).X);
+                                end
+                                xc=zeros(count,1); 
+                                yc=zeros(count,1);
+                                count=0;
+                                for j=1:length(contours),
+                                        nj=length(contours(j).X);
+                                        xc((count+j):(count+j+nj-1))=contours(j).X;
+                                        yc((count+j):(count+j+nj-1))=contours(j).Y;
+                                        count=count+nj;
+                                end
+
+                                disp(['projecting (X,Y)']);
+                                lat0= fix(mean(self.regions(i).CenLat)); long0=fix(mean(self.regions(i).CenLon));
+                                proj=laea(lat0,long0); self.regions(i).proj=proj;
+                                [xc,yc]=gdaltransform(xc,yc,'EPSG:4326',proj);
+
+                                disp(['plugging back into contours']);
+                                count=0;
+                                for j=1:length(contours),
+                                        nj=length(contours(j).X);
+                                        contours(j).X=xc((count+j):(count+j+nj-1));
+                                        contours(j).Y=yc((count+j):(count+j+nj-1));
+                                        count=count+nj;
+                                end
+                                disp(['saving new contours to disk']);
+                                shapewrite(contours,[self.root '/' self.regions(i).name '.laeaproj.shp']);
+                        end
+                end
+                %}}}
                 function disp(self) % {{{
                         disp(sprintf('   Glacier inventory:'));