CN117710606A - Urban area Gaussian projection remote sensing image acquisition method based on sky map - Google Patents
Urban area Gaussian projection remote sensing image acquisition method based on sky map Download PDFInfo
- Publication number
- CN117710606A CN117710606A CN202311514766.7A CN202311514766A CN117710606A CN 117710606 A CN117710606 A CN 117710606A CN 202311514766 A CN202311514766 A CN 202311514766A CN 117710606 A CN117710606 A CN 117710606A
- Authority
- CN
- China
- Prior art keywords
- image
- coordinate
- coordinates
- tile
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 49
- 239000011159 matrix material Substances 0.000 claims abstract description 41
- 238000004364 calculation method Methods 0.000 claims description 59
- 238000013507 mapping Methods 0.000 claims description 7
- 238000005516 engineering process Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 abstract description 5
- 238000012545 processing Methods 0.000 abstract description 2
- 238000013461 design Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 230000003442 weekly effect Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A30/00—Adapting or protecting infrastructure or their operation
- Y02A30/60—Planning or developing urban green infrastructure
Landscapes
- Image Processing (AREA)
Abstract
The invention relates to the technical field of remote sensing image processing, in particular to a method for acquiring a Gaussian projection remote sensing image of an urban area based on a sky map, which comprises the following steps: according to the central meridian of the city, the range from four coordinates to four coordinates and the zoom level of the sky map, calculating the Gaussian projection coordinates, the image resolution and the image size of the pixel origin of the city image, wherein the left and the downward directions of the image are respectively consistent with the east direction and the south direction of the Gaussian projection, and determining the pixel coordinate system of the city image matched with the Gaussian projection; dividing the urban image, numbering each divided square from left to right and from top to bottom; calculating and acquiring a single square image matrix according to the square numbers; and determining the grid range according to the requirement, and calculating all grid images to be solved to obtain the urban area image. The invention provides a method for establishing an urban image pixel coordinate system and dividing grids; the method for calculating the grid square Gaussian projection image by adopting an image remapping algorithm; the application of the sky map data in the urban planning and construction field is more convenient.
Description
Technical Field
The invention relates to the technical field of remote sensing image processing, in particular to a method for acquiring a Gaussian projection remote sensing image of an urban area based on a sky map.
Background
The national geographic information public service platform 'day map' (hereinafter referred to as 'day map') is a networked geographic information sharing and service portal, integrates geographic information departments of mapping of all levels of countries, provinces and cities (counties), and geographic information public service resources of related government departments, enterprises and public institutions, social groups and public, and provides authoritative, standard and unified online geographic information comprehensive services for various users. The sky map is widely applied in various fields of society, the data has very high situation and accuracy, taking the sky map Shandong as an example, the method realizes 15.8 ten thousand square kilometers full coverage of the land area with the resolution of better than 0.2 meter in the aspect of remote sensing images, the month coverage and weekly release of the 2 meter resolution images, and the image data supports the application of 3000 business systems in more than 50 fields in the whole province.
With the continuous deep application of the space map platform, the online and professional application mode of the space map platform also limits the application of the data in non-internet environment and non-GIS professional fields. Taking urban planning and construction fields as an example, the following problems mainly exist in the application of the sky map service: the working environment of the user is always off-line, and professional GIS software is not always available; equidistant cylindrical projection is adopted for the map tiles of the sky map, and an urban Gaussian projection coordinate system is generally adopted for urban reconnaissance design application, so that the map tiles and the urban Gaussian projection coordinate system cannot be directly matched; if the large scale remote sensing image of the whole urban area is used as an image, the large scale remote sensing image is not easy to use for a common computer. These factors limit the application of the space map remote sensing data.
Disclosure of Invention
The invention aims to overcome the defects in the prior art, and provides a method for acquiring the Gaussian projection remote sensing image of the urban area based on the sky map, so that the application of the sky map data in the urban planning and construction fields is more convenient.
The technical scheme adopted for solving the technical problems is as follows:
a method for acquiring Gaussian projection remote sensing images of urban areas based on a sky map comprises the following steps:
step 1, calculating Gaussian projection coordinates, image resolution and image size of an origin of pixels of a city image according to a central meridian of the city, a range from four coordinates to four coordinates and a zoom level of a sky map, wherein the left direction and the downward direction of the image are respectively consistent with the east direction and the south direction of Gaussian projection, so that a city image pixel coordinate system matched with the Gaussian projection is determined;
step 2, dividing the city image by taking 16384 multiplied by 16384 pixels as a unit to form a grid, numbering each divided grid from left to right and from top to bottom, wherein the format is as follows: column number_row number;
step 3, calculating and obtaining a single square image matrix according to the square numbers;
and 4, the whole urban image consists of all grids of the urban image grid, the grid range is determined according to the requirement, and all grid images to be solved are calculated to obtain the required urban area image.
Further, in step 1, the method for calculating the gaussian projection coordinates, the image resolution and the image size of the origin of the pixels of the urban image comprises the following steps:
step 11, setting the boundary redundancy value of the tile range to be 0 according to the four-to-range and the space map zoom level N of the city, and determining the calculation method of the space map tile range through two-point Gaussian projection coordinates to obtain the upper left corner tile number, the lower right corner tile number, the upper left corner vertex longitude and latitude coordinates and the lower right corner vertex longitude and latitude coordinates of the tile range;
step 12, taking the longitude and latitude of the central meridian and the latitude of the right lower corner vertex as coordinates to be combined with the longitude and latitude coordinates of the four vertices of the left upper corner vertex, the right lower corner vertex and the left lower corner vertex, and solving five-point Gaussian projection coordinates through a Gaussian projection forward calculation formula;
step 13, slave stepExtracting the minimum value X of the X coordinate from the calculation result in step 12 min Maximum X max And minimum Y of Y coordinate min Maximum value Y max Gaussian projection coordinates O (X 0 ,Y 0 ) Taking X 0 ＝X min ,Y 0 ＝Y max ；
Step 14, calculating the plane distance of the vertex line segments in the two meridian directions, wherein the image resolution dr is the quotient of the distance of the long side and the pixel number of the space map tile corresponding to the side, and the calculation formula is as follows:
dr＝md max /[(rb_row-lt_row+1)×256]；
wherein md max For a larger value of the plane distance of the vertex line segment in the meridian direction, rb_row is the number of the tile row in the lower right corner, lt_row is the number of the tile row in the upper left corner, and the size of the tile of the map is 256 multiplied by 256 pixels;
step 15, calculating the pixel size of the urban image, which is represented by the width W and the height H, wherein the calculation formula is as follows:
H＝[(Y max -Y min )/dr]+16384-[(Y max -Y min )/dr]％16384；
W＝[(X max -X min )/dr]+16384-[(X max -X min )/dr]％16384；
in [ of ]]For rounding,% is remainder, X min 、X max Respectively minimum and maximum of X coordinate, Y min 、Y max The minimum value and the maximum value of the Y coordinate are respectively, and the single grid size of the city image is 16384 multiplied by 16384.
Further, the step 11 includes:
step 111, taking the smaller value of the x-coordinate, the x-coordinate and the larger value of the y-coordinate, the x-coordinate and the y-coordinate as the upper left corner coordinates, taking the larger value of the x-coordinate, the x-coordinate and the y-coordinate as the lower right corner coordinates, and taking the smaller value of the x-coordinate, the y-coordinate and the y-coordinate as the lower right corner coordinates according to the known two-point gaussian projection coordinates (the x1, the y 1), (the x2, the y 2), and the x1, the y2, and the marking as the lower right corner coordinates: upper left corner coordinates (coordinate_x1, coordinate_y1) and lower right corner coordinates (coordinate_x2, coordinate_y2);
step 112, recording the x-coordinate and y-coordinate of four vertexes clockwise from the upper left corner point,
wherein, the x coordinate value is recorded as follows:
[coordinate_x1,coordinate_x2,coordinate_x2,coordinate_x1]；
the y coordinate values are noted:
[coordinate_y1,coordinate_y1,coordinate_y2,coordinate_y2]；
when the rectangle passes through the central meridian, taking a point (50000, chord_y1) as a fifth vertex, and calculating longitude and latitude coordinates of each vertex through a Gaussian projection back calculation formula;
step 113, taking out the minimum longitude value minl, the maximum longitude value maxl, the minimum latitude value minum and the maximum latitude value maxb from the calculation result of step 112 to obtain an upper left corner (minl, maxb) and a lower right corner (maxl, minub);
step 114, setting the zoom level to 18 or 20, and calculating the tile number (rb_col, rb_row) of the upper left corner (minl, maxb) and the lower right corner (maxl, minib) according to the tile number calculation formula;
step 115, setting a minimum redundancy value envelope of a tile range boundary, wherein a default value is 0, and finally determining a tile range of the sky map as follows: the longitude and latitude of the lower right corner tile of the upper left corner tile and the lower right corner tile are calculated through a tile longitude and latitude calculation formula, namely the longitude and latitude range of the tile range.
Further, in step 114, the tile number calculation formula is:
tileX＝[(λ+180)/360×2 N ]；
wherein [ the]Lambda is the rounding operationIn the form of a longitude and,for latitude, N is the zoom level, tileX is the tile column number, tileY is the tile row number.
Further, in step 115, the tile longitude and latitude calculation formula is:
λ＝(360×tileX)/2 N －180；
wherein lambda is the longitude and,for latitude, N is the zoom level, tileX is the tile column number, tileY is the tile row number.
Further, the step 3 includes:
step 31, calculating the city coordinate range of the square according to the square number col_row, wherein the city coordinate range is represented by the upper left corner vertex coordinates (lt_x, lt_y) and the lower right corner vertex coordinates (rb_x, rb_y), and the calculation formula is as follows:
lt_x＝X 0 +dr·col×16384；
lt_y＝Y 0 -dr·row×16384；
rb_x＝lt_x+dr·×16384；
rb_y＝lt_y-dr·×16384；
wherein (X) 0 ,Y 0 ) The pixel origin Gaussian projection coordinates of the urban image are represented by dr, and the resolution of the image is represented by dr;
step 32, a tile coordinate range is used for designating a zoom level as 18 or 20, a minimum redundancy value of a tile boundary as 1, a calculation method for determining a space map tile range by a two-point Gaussian projection coordinate is used for calculating a tile range of a minimum source image for mapping the square, and a calculation result comprises an upper left corner tile number (lt_col, lt_row) and a lower right corner tile number (rb_col, rb_row);
step 33, downloading the map tiles by adopting a crawler technology according to the tile range, splicing the tiles row by row into rows by adopting a numpy hstack, vstack method, and then splicing the tiles row by row into a row to acquire a source image matrix;
step 34, calculating a remapping matrix mapped from the source image matrix to the square image matrix;
and 35, setting an interpolation mode as bilinear interpolation according to the obtained remapping matrix, and carrying out remapping calculation on the source image through an image remapping algorithm to obtain a square image matrix.
Further, in step 35, the remapping calculation formula is:
dst_image＝remap(src_image,Matrix_X,Matrix_Y,Inter_Linear)；
wherein dst_image represents a target image Matrix, src_image represents a source image Matrix, matrix_x represents a remapping Matrix of each pixel of the target image in the direction of the source image u, matrix_y represents a remapping Matrix of each pixel of the target image in the direction of the source image v, inter_linear represents a bilinear interpolation mode, and remap represents a remapping calculation method.
Further, in step 34, the method for calculating the remapping matrix includes the following steps:
step 341, calculating the city coordinates (x, y) of any pixel (du, dv) of the target image, where the calculation formula is:
x＝x 0 +du·dr；
y＝y 0 -dv·dr；
in (x) 0 ,y 0 ) The top left corner vertex of the square image is matched with the coordinate of the urban Gaussian projection coordinate system;
step 342, calculating the earth coordinates (B, L) of the point by a gaussian projection back calculation formula;
step 343, calculating the pixel coordinates (su, sv) of the geodetic coordinates (B, L) in the source image:
su＝(L-L 0 )/sr；
sv＝(B 0 -B)/sr；
wherein B is 0 、L 0 The latitude and longitude coordinates of the top left vertex of the source image, sr is the image resolution of the source image, and the zoom level of the sky map is setSource image resolution sr= 5.3644180297851516 ×10 at 18 levels -6 The image resolution sr of the source image is 1/4 of the 18-level value when the zoom level is set to 20 levels;
in step 344, the mapping positions of all pixels of the target pixel Matrix in the source image are calculated according to the pixel coordinates by adopting the methods of steps 341 to 343, the corresponding longitudinal direction is marked as matrix_x (X, Y), and the latitudinal direction is marked as matrix_y (X, Y), so as to obtain the remapped matrices matrix_x and matrix_y.
The invention has the technical effects that:
compared with the prior art, the urban area Gaussian projection remote sensing image acquisition method based on the space map provides a method for establishing an urban image pixel coordinate system and dividing grids; the method for calculating the grid Gaussian projection image by adopting an image remapping algorithm is provided on the basis of a sky map (equidistant cylindrical projection) map service; the method has the advantages that the Gaussian projection image production in the urban area is independent of professional GIS software, the basic environment requirement is low, the operation speed is high, and the data can be conveniently obtained and updated in batches by building a grid; the application of the sky map data in the fields of city planning, construction and the like is more convenient, and the application scene of the sky map is enriched.
Drawings
FIG. 1 is a flow chart of a method for acquiring Gaussian projection remote sensing images of urban areas based on a sky map;
FIG. 2 is a schematic diagram of a coordinate representation of a determined rectangular vertex in accordance with the present invention;
FIG. 3 is a grid image of an urban image of Jinan city of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention are clearly and completely described below with reference to the accompanying drawings of the specification.
Example 1:
as shown in fig. 1, the method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to the embodiment includes the following steps: determining a pixel coordinate system of the urban image matched with Gaussian projection; dividing the city image by taking 16384 multiplied by 16384 pixels as a unit to form a grid and numbering; calculating and acquiring a single square image matrix according to the square numbers; and determining the grid range according to the requirement, and calculating all grid images to be solved to obtain the required urban area image.
Taking the eastern city of atanan as an example, the central meridian of the city is known to be the east meridian 117 DEG, the four-to-range is [429469.9436,4156123.4832,587153.4291,3984928.6951], the day map zoom level N is 18, and the tile boundary redundancy is 0. In the work of urban investigation design and the like, the adopted plane coordinate system is usually marked as x direction in the east-west direction and y direction in the north-south direction, so the embodiment adopts the metering method in the calculation process.
Step 1, calculating Gaussian projection coordinates, image resolution and image size of an origin of pixels of the urban image according to a central meridian of the urban, a range from four coordinates to four coordinates and a zoom level of a sky map, wherein the left direction and the downward direction of the image are respectively consistent with the east direction and the south direction of Gaussian projection, so that an urban image pixel coordinate system matched with the Gaussian projection is determined.
Specifically, the method for calculating the Gaussian projection coordinates, the image resolution and the image size of the pixel origin of the urban image comprises the following steps:
step 11, setting the boundary redundancy value of the tile range to be 0 according to the scale level N of the map of the city from four to four, and determining the calculation method of the map of the sky by two-point Gaussian projection coordinates to obtain the upper left corner tile number, the lower right corner tile number, the upper left corner vertex longitude and latitude coordinate and the lower right corner vertex longitude and latitude coordinate of the tile range (the Gaussian projection coordinate system which is suitable for translation of northern hemisphere and east-west coordinate for 500 km and is not added with a sign);
calculated by using an Jinan city example, the upper left corner tile number (215687,38202), the lower right corner tile number (216986,39328), the upper left corner vertex longitude and latitude coordinate (116.20101928710938,37.53753662109375) and the lower right corner vertex longitude and latitude coordinate (117.98629760742188,35.989837646484375);
step 12, taking the longitude and latitude of the central meridian and the latitude of the right lower corner vertex as coordinates to be combined with the longitude and latitude coordinates of the four vertices of the left upper corner vertex, the right lower corner vertex and the left lower corner vertex, and solving five-point Gaussian projection coordinates through a Gaussian projection forward calculation formula;
calculated by using Jinan city example, the X coordinate values of the Gaussian projection coordinates of the five points are as follows:
[429384.5615034835,587171.2477544363,588941.0089387904,427950.9820019439,500000.0], Y-coordinate values are:
[4156468.231220618,4156625.396751042,3984864.944948826,3984710.2832620037,3984415.0666628554]；
step 13, extracting the minimum value X of the X coordinate from the calculation result of step 12 min Maximum X max And minimum Y of Y coordinate min Maximum value Y max Gaussian projection coordinates O (X 0 ,Y 0 ) Taking X 0 ＝X min ,Y 0 ＝Y max ；
Calculated with the Jinan city example, the Gaussian projection coordinate of the pixel origin, O (427950.9820019439,4156625.396751042);
step 14, calculating the plane distance of the vertex line segments in the two meridian directions, wherein the image resolution dr is the quotient of the distance of the long side and the pixel number of the space map tile corresponding to the side, and the calculation formula is as follows:
dr＝md max /[(rb_row-lt_row+1)×256]；
wherein md max For a larger value of the plane distance of the vertex line segment in the meridian direction, rb_row is the number of the tile row in the lower right corner, lt_row is the number of the tile row in the upper left corner, and the size of the tile of the map is 256 multiplied by 256 pixels;
calculated with the atanan city example, the image resolution dr= 0.5953636904822707 meters/pixel;
step 15, calculating the pixel size of the urban image, which is represented by the width W and the height H, wherein the calculation formula is as follows:
H＝[(Y max -Y min )/dr]+16384-[(Y max -Y min )/dr]％16384；
W＝[(X max -X min )/dr]+16384-[(X max -X min )/dr]％16384；
in [ of ]]For rounding,% is remainder, X min 、X max Respectively minimum and maximum of X coordinate, Y min 、Y max The minimum value and the maximum value of the Y coordinate are respectively, and the single grid size of the city image is 16384 multiplied by 16384.
Calculated by using the Jinan city example, the city image pixel size width w= 278528 and the height= 294912;
the step 11 includes:
step 111, taking the smaller value of the x-coordinate, the x-coordinate and the larger value of the y-coordinate, the x-coordinate and the y-coordinate as the upper left corner coordinates, taking the larger value of the x-coordinate, the x-coordinate and the y-coordinate as the lower right corner coordinates, and taking the smaller value of the x-coordinate, the y-coordinate and the y-coordinate as the lower right corner coordinates according to the known two-point gaussian projection coordinates (the x1, the y 1), (the x2, the y 2), and the x1, the y2, and the marking as the lower right corner coordinates: upper left corner coordinates (coordinate_x1, coordinate_y1) and lower right corner coordinates (coordinate_x2, coordinate_y2);
step 112, the rectangle defined by the two points of the upper left corner and the lower right corner, records the x-coordinates and y-coordinates of the four vertices clockwise from the upper left corner point, as shown in figure 2,
wherein, the x coordinate value is recorded as follows:
[coordinate_x1,coordinate_x2,coordinate_x2,coordinate_x1]；
the y coordinate values are noted:
[coordinate_y1,coordinate_y1,coordinate_y2,coordinate_y2]；
when the rectangle passes through the central meridian, taking a point (50000, chord_y1) as a fifth vertex, and calculating longitude and latitude coordinates of each vertex through a Gaussian projection back calculation formula;
step 113, taking out the minimum longitude value minl, the maximum longitude value maxl, the minimum latitude value minum and the maximum latitude value maxb from the calculation result of step 112 to obtain an upper left corner (minl, maxb) and a lower right corner (maxl, minub);
step 114, setting the zoom level to 18 or 20, and calculating the tile number (rb_col, rb_row) of the upper left corner (minl, maxb) and the lower right corner (maxl, minib) according to the tile number calculation formula;
step 115, setting a minimum redundancy value envelope of a tile range boundary, wherein a default value is 0, and finally determining a tile range of the sky map as follows: the longitude and latitude of the lower right corner tile of the upper left corner tile and the lower right corner tile are calculated through a tile longitude and latitude calculation formula, namely the longitude and latitude range of the tile range.
In step 114, the tile number calculation formula is:
tileX＝[(λ+180)/360×2 N ]；
in step 115, the tile longitude and latitude calculation formula is:
λ＝(360×tileX)/2 N －180；
wherein [ the]For a rounding operation, λ is the longitude,for latitude, N is the zoom level, tileX is the tile column number, tileY is the tile row number.
Step 2, dividing the city image by taking 16384 multiplied by 16384 pixels as a unit to form a grid, numbering each divided grid from left to right and from top to bottom, wherein the format is as follows: column number_row number;
calculated by the Jinan city example, the grid number of the top left vertex is "0_0", the grid size is 278528 × 294912, and the grid numbers of the top right, bottom right and bottom left vertices are "16_0", "16_17" and "17_0" in sequence, as shown in fig. 3.
Step 3, calculating and obtaining a single square image matrix according to the square numbers; specifically, the step 3 includes:
step 31, calculating the city coordinate range of the square according to the square number col_row, wherein the city coordinate range is represented by the upper left corner vertex coordinates (lt_x, lt_y) and the lower right corner vertex coordinates (rb_x, rb_y), and the calculation formula is as follows:
lt_x＝X 0 +dr·col×16384；
lt_y＝Y 0 -dr·row×16384；
rb_x＝lt_x+dr·×16384；
rb_y＝lt_y-dr·×16384；
wherein (X) 0 ,Y 0 ) The pixel origin Gaussian projection coordinates of the urban image are represented by dr, and the resolution of the image is represented by dr;
step 32, a tile coordinate range is used for designating a zoom level as 18 or 20, a minimum redundancy value of a tile boundary as 1, a calculation method for determining a space map tile range by a two-point Gaussian projection coordinate is used for calculating a tile range of a minimum source image for mapping the square, and a calculation result comprises an upper left corner tile number (lt_col, lt_row) and a lower right corner tile number (rb_col, rb_row);
step 33, downloading the map tiles by adopting a crawler technology according to the tile range, splicing the tiles row by row into rows by adopting a numpy hstack, vstack method, and then splicing the tiles row by row into a row to acquire a source image matrix;
step 34, calculating a remapping matrix mapped from the source image matrix to the square image matrix;
and 35, setting an interpolation mode as bilinear interpolation according to the obtained remapping matrix, and carrying out remapping calculation on the source image through an image remapping algorithm to obtain a square image matrix.
In step 35, the remapping calculation formula is:
dst_image＝remap(src_image,Matrix_X,Matrix_Y,Inter_Linear)；
wherein dst_image represents a target image Matrix, src_image represents a source image Matrix, matrix_x represents a remapping Matrix of each pixel of the target image in the direction of the source image u, matrix_y represents a remapping Matrix of each pixel of the target image in the direction of the source image v, inter_linear represents a bilinear interpolation mode, and remap represents a remapping calculation method. The invention adopts a remap method in a python language openCV library during specific programming calculation.
In step 34, the method for calculating the remapping matrix includes the following steps:
step 341, calculating the city coordinates (x, y) of any pixel (du, dv) of the target image, where the calculation formula is:
x＝x 0 +du·dr；
y＝y 0 -dv·dr；
in (x) 0 ,y 0 ) The top left corner vertex of the square image is matched with the coordinate of the urban Gaussian projection coordinate system;
step 342, calculating the earth coordinates (B, L) of the point by a gaussian projection back calculation formula;
step 343, calculating the pixel coordinates (su, sv) of the geodetic coordinates (B, L) in the source image:
su＝(L-L 0 )/sr；
sv＝(B 0 -B)/sr；
wherein B is 0 、L 0 The latitude and longitude coordinates of the top left vertex of the source image (the value given by step 115), sr is the image resolution of the source image, and the source image resolution sr= 5.3644180297851516 ×10 when the scale level of the sky map is set to 18 levels -6 The image resolution sr of the source image is 1/4 of the 18-level value when the zoom level is set to 20 levels;
in step 344, the mapping positions of all pixels of the target pixel Matrix in the source image are calculated according to the pixel coordinates by adopting the methods of steps 341 to 343, the corresponding longitudinal direction is marked as matrix_x (X, Y), and the latitudinal direction is marked as matrix_y (X, Y), so as to obtain the remapped matrices matrix_x and matrix_y.
And 4, the whole urban image consists of all grids of the urban image grid, the grid range is determined according to the requirement, and all grid images to be solved are calculated to obtain the required urban area image.
The invention provides a method for establishing an urban image pixel coordinate system and dividing grids; the method for calculating the grid Gaussian projection image based on the equidistant cylindrical projection map service of the sky map by adopting an image remapping algorithm is provided; the method has the advantages that the production of the Gaussian projection image data in the urban area is independent of professional GIS software, the basic environment requirement is low, the operation speed is high, and the data can be conveniently obtained and updated in batches by building a grid; the application of the sky map data in the fields of city planning, construction and the like is more convenient, and the application scene of the sky map is enriched.
The above embodiments are merely examples of the present invention, and the scope of the present invention is not limited to the above embodiments, and any suitable changes or modifications made by those skilled in the art, which are consistent with the claims of the present invention, shall fall within the scope of the present invention.
Claims (8)
1. A method for acquiring a Gaussian projection remote sensing image of an urban area based on a sky map is characterized by comprising the following steps of: the method comprises the following steps:
step 1, calculating Gaussian projection coordinates, image resolution and image size of an origin of pixels of a city image according to a central meridian of the city, a range from four coordinates to four coordinates and a zoom level of a sky map, wherein the left direction and the downward direction of the image are respectively consistent with the east direction and the south direction of Gaussian projection, so that a city image pixel coordinate system matched with the Gaussian projection is determined;
step 2, dividing the city image by taking 16384 multiplied by 16384 pixels as a unit to form a grid, numbering each divided grid from left to right and from top to bottom, wherein the format is as follows: column number_row number;
step 3, calculating and obtaining a single square image matrix according to the square numbers;
and 4, the whole urban image consists of all grids of the urban image grid, the grid range is determined according to the requirement, and all grid images to be solved are calculated to obtain the required urban area image.
2. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 1, wherein the method comprises the following steps of: in step 1, the method for calculating the gaussian projection coordinates, the image resolution and the image size of the origin of the pixels of the urban image comprises the following steps:
step 11, setting the boundary redundancy value of the tile range to be 0 according to the four-to-range and the space map zoom level N of the city, and determining the calculation method of the space map tile range through two-point Gaussian projection coordinates to obtain the upper left corner tile number, the lower right corner tile number, the upper left corner vertex longitude and latitude coordinates and the lower right corner vertex longitude and latitude coordinates of the tile range;
step 12, taking the longitude and latitude of the central meridian and the latitude of the right lower corner vertex as coordinates to be combined with the longitude and latitude coordinates of the four vertices of the left upper corner vertex, the right lower corner vertex and the left lower corner vertex, and solving five-point Gaussian projection coordinates through a Gaussian projection forward calculation formula;
step 13, extracting the minimum value X of the X coordinate from the calculation result of step 12 min Maximum X max And minimum Y of Y coordinate min Maximum value Y max Gaussian projection coordinates O (X 0 ,Y 0 ) Taking X 0 ＝X min ,Y 0 ＝Y max ；
Step 14, calculating the plane distance of the vertex line segments in the two meridian directions, wherein the image resolution dr is the quotient of the distance of the long side and the pixel number of the space map tile corresponding to the side, and the calculation formula is as follows:
dr＝md max /[(rb_row-lt_row+1)×256]；
wherein md max For a larger value of the plane distance of the vertex line segment in the meridian direction, rb_row is the number of the tile row in the lower right corner, lt_row is the number of the tile row in the upper left corner, and the size of the tile of the map is 256 multiplied by 256 pixels;
step 15, calculating the pixel size of the urban image, which is represented by the width W and the height H, wherein the calculation formula is as follows:
H＝[(Y max -Y min )/dr]+16384-[(Y max -Y min )/dr]％16384；
W＝[(X max -X min )/dr]+16384-[(X max -X min )/dr]％16384；
in [ of ]]For rounding,% is remainder, X min 、X max Respectively minimum and maximum of X coordinate, Y min 、Y max The minimum value and the maximum value of the Y coordinate are respectively, and the single grid size of the city image is 16384 multiplied by 16384.
3. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 2, wherein the method comprises the following steps of: the step 11 includes:
step 111, taking the smaller value of the x-coordinate, the x-coordinate and the larger value of the y-coordinate, the x-coordinate and the y-coordinate as the upper left corner coordinates, taking the larger value of the x-coordinate, the x-coordinate and the y-coordinate as the lower right corner coordinates, and taking the smaller value of the x-coordinate, the y-coordinate and the y-coordinate as the lower right corner coordinates according to the known two-point gaussian projection coordinates (the x1, the y 1), (the x2, the y 2), and the x1, the y2, and the marking as the lower right corner coordinates: upper left corner coordinates (coordinate_x1, coordinate_y1) and lower right corner coordinates (coordinate_x2, coordinate_y2);
step 112, recording the x-coordinate and y-coordinate of four vertexes clockwise from the upper left corner point,
wherein, the x coordinate value is recorded as follows:
[coordinate_x1,coordinate_x2,coordinate_x2,coordinate_x1]；
the y coordinate values are noted:
[coordinate_y1,coordinate_y1,coordinate_y2,coordinate_y2]；
when the rectangle passes through the central meridian, taking a point (50000, chord_y1) as a fifth vertex, and calculating longitude and latitude coordinates of each vertex through a Gaussian projection back calculation formula;
step 113, taking out the minimum longitude value minl, the maximum longitude value maxl, the minimum latitude value minum and the maximum latitude value maxb from the calculation result of step 112 to obtain an upper left corner (minl, maxb) and a lower right corner (maxl, minub);
step 114, setting the zoom level to 18 or 20, and calculating the tile number (rb_col, rb_row) of the upper left corner (minl, maxb) and the lower right corner (maxl, minib) according to the tile number calculation formula;
step 115, setting a minimum redundancy value envelope of a tile range boundary, wherein a default value is 0, and finally determining a tile range of the sky map as follows: the longitude and latitude of the lower right corner tile of the upper left corner tile and the lower right corner tile are calculated through a tile longitude and latitude calculation formula, namely the longitude and latitude range of the tile range.
4. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 3, wherein the method comprises the following steps of: in step 114, the tile number calculation formula is:
tileX＝[(λ+180)/360×2 N ]；
wherein [ the]For a rounding operation, λ is the longitude,for latitude, N is the zoom level, tileX is the tile column number, tileY is the tile row number.
5. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 3, wherein the method comprises the following steps of: in step 115, the tile longitude and latitude calculation formula is:
λ＝(360×tileX)/2 N －180；
wherein lambda is the longitude and,for latitude, N is the zoom level, tileX is the tile column number, tileY is the tile row number.
6. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 3, wherein the method comprises the following steps of: the step 3 comprises the following steps:
step 31, calculating the city coordinate range of the square according to the square number col_row, wherein the city coordinate range is represented by the upper left corner vertex coordinates (lt_x, lt_y) and the lower right corner vertex coordinates (rb_x, rb_y), and the calculation formula is as follows:
lt_x＝X 0 +dr·col×16384；
lt_y＝Y 0 -dr·row×16384；
rb_x＝lt_x+dr·×16384；
rb_y＝lt_y-dr·×16384；
wherein (X) 0 ,Y 0 ) The pixel origin Gaussian projection coordinates of the urban image are represented by dr, and the resolution of the image is represented by dr;
step 32, a tile coordinate range is used for designating a zoom level as 18 or 20, a minimum redundancy value of a tile boundary as 1, a calculation method for determining a space map tile range by a two-point Gaussian projection coordinate is used for calculating a tile range of a minimum source image for mapping the square, and a calculation result comprises an upper left corner tile number (lt_col, lt_row) and a lower right corner tile number (rb_col, rb_row);
step 33, downloading the map tiles by adopting a crawler technology according to the tile range, splicing the tiles row by row into rows by adopting a numpy hstack, vstack method, and then splicing the tiles row by row into a row to acquire a source image matrix;
step 34, calculating a remapping matrix mapped from the source image matrix to the square image matrix;
and 35, setting an interpolation mode as bilinear interpolation according to the obtained remapping matrix, and carrying out remapping calculation on the source image through an image remapping algorithm to obtain a square image matrix.
7. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 6, wherein the method comprises the following steps of: in step 35, the remapping calculation formula is:
dst_image＝remap(src_image,Matrix_X,Matrix_Y,Inter_Linear)；
wherein dst_image represents a target image Matrix, src_image represents a source image Matrix, matrix_x represents a remapping Matrix of each pixel of the target image in the direction of the source image u, matrix_y represents a remapping Matrix of each pixel of the target image in the direction of the source image v, inter_linear represents a bilinear interpolation mode, and remap represents a remapping calculation method.
8. The method for acquiring the urban area gaussian projection remote sensing image based on the sky map according to claim 6, wherein the method comprises the following steps of: in step 34, the method for calculating the remapping matrix includes the following steps:
step 341, calculating the city coordinates (x, y) of any pixel (du, dv) of the target image, where the calculation formula is:
x＝x 0 +du·dr；
y＝y 0 -dv·dr；
in (x) 0 ,y 0 ) The top left corner vertex of the square image is matched with the coordinate of the urban Gaussian projection coordinate system;
step 342, calculating the earth coordinates (B, L) of the point by a gaussian projection back calculation formula;
step 343, calculating the pixel coordinates (su, sv) of the geodetic coordinates (B, L) in the source image:
su＝(L-L 0 )/sr；
sv＝(B 0 -B)/sr；
wherein B is 0 、L 0 The latitude and longitude coordinates of the top left vertex of the source image, sr is the image resolution of the source image, and when the scale level of the sky map is set to 18 levels, sr= 5.3644180297851516 ×10 is the image resolution of the source image -6 A degree, when the zoom level is set to 20, the source image resolution sr is 1/4 of the 18-level value;
in step 344, the mapping positions of all pixels of the target pixel Matrix in the source image are calculated according to the pixel coordinates by adopting the methods of steps 341 to 343, the corresponding longitudinal direction is marked as matrix_x (X, Y), and the latitudinal direction is marked as matrix_y (X, Y), so as to obtain the remapped matrices matrix_x and matrix_y.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311514766.7A CN117710606A (en) | 2023-11-15 | 2023-11-15 | Urban area Gaussian projection remote sensing image acquisition method based on sky map |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311514766.7A CN117710606A (en) | 2023-11-15 | 2023-11-15 | Urban area Gaussian projection remote sensing image acquisition method based on sky map |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117710606A true CN117710606A (en) | 2024-03-15 |
Family
ID=90150500
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311514766.7A Pending CN117710606A (en) | 2023-11-15 | 2023-11-15 | Urban area Gaussian projection remote sensing image acquisition method based on sky map |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117710606A (en) |
-
2023
- 2023-11-15 CN CN202311514766.7A patent/CN117710606A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109493407A (en) | Realize the method, apparatus and computer equipment of laser point cloud denseization | |
| CN108921943B (en) | Road three-dimensional model modeling method based on lane-level high-precision map | |
| CN109801219A (en) | GIS data bearing calibration and device towards Online Map superposition | |
| EP0439543A1 (en) | A method of storage and retrieval of digital map data based upon a tessellated geoid system. | |
| CN107967702A (en) | The thermodynamic chart method for visualizing of reverse color applying drawing | |
| Chen et al. | Geographical data acquisition | |
| Khakimova et al. | DEVELOPMENT OF CADASTRAL MAPS AND PLANS IN THE GEOINFORMATION SYSTEM | |
| CN109801371B (en) | Network three-dimensional electronic map construction method based on Cesium | |
| CN110516015B (en) | Method for manufacturing geographical PDF map based on map graphic data and DLG | |
| CN103150328A (en) | Integration method for automatically drawing maps by utilizing database | |
| CN105334525A (en) | Geographic information display method based on augmented reality technology | |
| CN103700063B (en) | Topography integration quick mapping method based on high definition satellite image | |
| US20100066740A1 (en) | Unified spectral and Geospatial Information Model and the Method and System Generating It | |
| CN115471619A (en) | City three-dimensional model construction method based on stereo imaging high-resolution satellite image | |
| CN107451958A (en) | A kind of angular transition algorithm towards ArcGIS polar regions vector field visualization | |
| CN117710606A (en) | Urban area Gaussian projection remote sensing image acquisition method based on sky map | |
| CN115186347B (en) | Building CityGML modeling method combining house type plan view and inclination model | |
| CN110634156A (en) | Online registration method of indoor plane graph | |
| CN115859414A (en) | Cross-coordinate system use method for global scale geographic information base map | |
| CN112525208B (en) | Method, device and equipment for quickly constructing urban road map | |
| CN115409962A (en) | Method for constructing coordinate system in illusion engine, electronic equipment and storage medium | |
| Petrie | Digital mapping technology: procedures and applications | |
| CN111339232B (en) | Method for applying weather radar standard format base data in geographic information system | |
| CN113066181A (en) | Terrain simulation method based on satellite images and digital elevation data | |
| CN105279341B (en) | A method of interception Google Earth image directly applies to engineering design |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |