texts: gdal cannot work when using openMP in c++ for reading a dataset?

texts:
void generate( GDALDataset  *ds)
{
    int tcount = (1 + abs(tmaxx - tminx))*(1 + abs(tmaxy - tminy));
    GDALDriver* jpegDriver = GetGDALDriverManager()->GetDriverByName("PNG");
    int tz = tmaxz;
    for (int ty = tmaxy; ty >= tminy; ty--)
    {
        #pragma omp parallel for schedule(dynamic, 1)
        for (int tx = tminx; tx <= tmaxx; tx++){

            // TODO tilefilename is exits
            int googleY = (pow(2, tz) - 1) - ty;
            string tileDir = "D:\\test2\\" + to_string(tz) + "\\" + to_string(tx);
            string tilefilename = tileDir + "\\" +to_string(googleY) + "." + tileext;
            create_dir(tileDir.c_str());
            Envelope b = mercator->TileBounds(tx, ty, tz);
            int * rwxy;
            rwxy = geo_query(b,querysize);
            int rx, ry, rxsize, rysize, wx, wy, wxsize, wysize;
            rx = *(rwxy+0);
            ry = *(rwxy + 1);
            rxsize = *(rwxy + 2);
            rysize = *(rwxy + 3);
            wx = *(rwxy + 4);
            wy = *(rwxy + 5);
            wxsize = *(rwxy + 6);
            wysize = *(rwxy + 7);

            GDALDataset *dstile;
            char **papszOptions = NULL;
            dstile = mem_drv->Create("temp", tilesize, tilesize, 3, GDT_Byte, papszOptions);

            int bandNum;
            bandNum = ds->GetRasterCount();    //波段数
            bandNum = 3;
            int depth = GDALGetDataTypeSize(ds->GetRasterBand(1)->GetRasterDataType());//8;    //图像深度
            size_t imgBufNum = (size_t)wxsize * wysize * bandNum*depth;
            GByte *data = new GByte[imgBufNum];

            ds->RasterIO(GF_Read, rx, ry, rxsize, rysize, data, wxsize, wysize,
                GDT_Byte, bandNum, nullptr, bandNum*depth, wxsize*bandNum*depth, depth);
            if (tilesize == querysize)
            {
                //写入
                dstile->RasterIO(GF_Write, wx, wy, wxsize, wysize, data, wxsize, wysize,
                GDT_Byte, bandNum, nullptr, bandNum*depth, wxsize*bandNum*depth, depth);
                dstile->GetRasterBand(4)->SetNoDataValue(0);
            }
            else
            {

                GDALDataset *dsquery;
                dsquery = mem_drv->Create("temp", querysize, querysize, 3, GDT_Byte, papszOptions);
                dsquery->RasterIO(GF_Write, wx, wy, wxsize, wysize, data, wxsize, wysize,
                    GDT_Byte, bandNum, nullptr, bandNum*depth, wxsize*bandNum*depth, depth);
                scale_query_to_tile(dsquery, dstile);
                GDALClose((GDALDatasetH)dsquery);

            }

            delete[] data;
            data = NULL;
            clock_t start, finish;
            start = clock();
            GDALDataset *poDstDS;
            poDstDS = jpegDriver->CreateCopy(tilefilename.c_str(), dstile, FALSE, NULL, NULL, NULL);
            //Once we're done, close properly the dataset */
            if (poDstDS != NULL)
                GDALClose((GDALDatasetH)poDstDS);
            GDALClose((GDALDatasetH)dstile);
        }
    }
    //GetLocalTime(&st);
    }

The above is the main code, and the error happen in"ds->RasterIO(GF_Read, ...)". and the error is " The error is “An access violation occurs when the location 0x000001D93A661000 is written”. "  How can I fixed it or any other better solution?


texts:
I am developing earthquake visualization web-app for arctic zone using Google Maps JavaScript API. During the development I found Google Maps doesn't display map above approximately N85 degree, if you build Google maps in your site (on their site arctic displays correctly in satellite mode). It's fine if you don't work with arctic.

In my case this is a problem, because I have earthquake data above N85 to display. You can see my current situation on a picture below.

Maybe someone already had same problem and knows certain how to fix it. I tried to set another projection for map (it also fix object distortion problem), but don't find high resolution tile for arctic. Hope there is another solutions.


texts: Calculating area of solar radiation using ArcGIS Desktop?

texts: How can I obtain area solar radiation using ArcGIS Desktop?


texts: Having different annotations for same feature in ArcGIS Desktop?

texts:
I've been given a dataset that has two different annotations set for a single feature which I can flip between. See below screenshot which allows me to set either a blue or red leader line for the same annotation.

I need to create my own data set from scratch that allows me to do this. I only seem to have the default annotation. How do I create an annotation other than the default.


texts:
I am looking for a way to convert H3 hex-ids (e.g., '87422c2a9ffffff') to Polygon geometry. In R, we would use:
h3_to_polygon()

But I am having a hard time finding the equivalent in Python.


texts: Using GDAL to strip nodata in border regions only?

texts:
I am trying to create a script that given a tiff file with islands of data (which can be disjoint) and black pixels of nodata, can generate polygons (or new tiffs) that only wrap around the islands of data. The issue I am facing is that my nodata value is automatically 0 (0 0 0 I think actually, or black), and some of my islands also contain the value 0 inside of them, what is the best way to accomplish this in gdal?

For example, if I run a command like:

gdalwarp -srcnodata 0 -dstnodata 0 orig-ignore-black.tif black-nodata.tif

my resulting tiff has pockets inside the data regions that have been stripped because they had a 0 value.


texts: georeferencing raster data with unknown CRS

texts: I’m new to QGIS and georeferencing raster data. I have a raster dataset with unknown CRS. The coordinates for goerefencing the raster data with reference points (x = 589 Y = 1483 / x = 500 Y = 1500) from an archaeological digging raster. The digging site is in Switzerland. I should georefernce that raster data, how do I have to set the proper CRS? I tried several time but failed and I have no clue how to proceed.


texts:
I'm trying to query my PostGIS database for ID of people who did X but did not do Y. I've attempted the query below, but it does not return any results. I'm very new to SQL so I may be overlooking something. Any thoughts? Thanks for your help.

SELECT *
FROM Table1 INNER JOIN Table2
ON Table1.ID = Table2.ID
WHERE Table1.Date = '2008-11-04'
AND NOT EXISTS (
    SELECT *
    FROM Table1 INNER JOIN Table2
    ON Table1.ID = Table2.ID
    WHERE Table1.Date = '2012-11-06')


texts: Creating custom SLDs for Heatmaps using GeoServer

texts:
im trying to figure out whether it is possible to create a custom SLD to style heatmaps provided by GeoServer depending on dynamic data input. Already checked out the Rest api and some related documentations, but i did not find a clear answer yet. I am thinking about stuff like a dynamic colour ramp depending on how many features there are to be displayed, also like changing the colours of the table by an input database.

Are there any ways this can be done?


texts: Adding or removing values in style tab (layer properties) of QGIS?

texts:
I'm new to QGIS and have recently been testing the software as my org plan to migrate from ArcGIS to QGIS in the near future.

I often use the 'Add values', 'Remove values'  and the group values feature in ArcGIS (Symobology¬ Categories¬ Unique Values) to display Symbology for selected values.

Is it possible to do this in QGIS?

With categorized selected in the style menu the only options I can see is to add or delete.


texts: Saving after using fet.setFeatureId in PyQGIS?

texts:
By using the following code, I want to change the id of the line when saving the line. But the setFeatureId command, gives a  false when doing the commitChanges().

Any idea how to change the id?

vl = QgsVectorLayer("MultiLineString", layerName, "memory")
pr = vl.dataProvider()
vl.startEditing()
fet = QgsFeature()
fet.setGeometry( QgsGeometry.fromMultiPolyline( [[ QgsPoint(x3,y3),QgsPoint(x4,y4) ]] ))
fet.setFeatureId(10001)
pr.addFeatures( [ fet ] )
vl.commitChanges()

I added an attribute with the id i wanted instead of changing the existing id.


texts: QGIS 3.10 atlas / composer, avoid missing text edge of map?

texts:
How to avoid missing texts near edge between map images?

I have vector topographic map which need to cut through map names.
When I create atlas project using Composer, those cutting text are missing from both maps. How to solve this problem?

I am using 3.10 version


texts:
I'm running QGIS 2.2 on Ubuntu 14.04. Trying to project data gives me a hard time right now.

The last time I did that was in 2.0 and I basically changed the projection by saving the data as a new file with a certain layer CRS and opening it in another QGIS project. That doesn't work right now, even when on-the-fly transformation is disabled in the project properties.

I've been trying project some vector data with WGS 84 geographical coordinates on WGS 84 UTM 32 N (Germany) but see no results from changes I make in the project or the layer properties.


texts: Route creation with specific length and elevation change

texts:
I'm currently working on a project for my computer science study. The goal is to make a service that allows the user to enter a distance and elevation change which will then result in a random route (either with the starting point being the endpoint at the same time or not).

Our "stakeholder" wants the routes to always be different so it doesn't get boring while he's jogging / biking haha.

My current progress:

I wrote a small library in golang to calculate the elevation at a certain gps coordinate (with SRTM data). -> It's not relevant right now, just FYI :)
I loaded the OSM file of switzerland into a local postgres DB, which is working. (I'm using the SRID 21781 for switzerland).
--> I was able to calculate a distance between two points but it's just a direct line without respecting streets etc.

My idea is to make some kind of a grid with a distance of 20 - 50 meters between all the points (looking like a chess board haha). Then I calculate the distance between those points by choosing different lines (I'm thinking of the table planet_osm_line). Like that I can slowly build up a route. Since the elevation change in those 20-50 meters isn't rapidly growing, I could take the average elevation between A and B to get an overall elevation difference.

I saw that there is a tool called pgRouting. But as far as I've seen, it calculates a route between two points with Dijkstras algorithm (which wouldn't help me get a route with length x, but the shortest one).

I'm not sure if this is a good approach and how I should actually realize it with the PostGIS tool.

Do you guys have some tips that can guide me in the right direction?


texts: Visualizing results of an aggregate function query

texts:
I want to visualize the results of an aggregate function in CartoDB. For example I want to understand and visualize the top 20 owners of land by acres with the following query.  As single owner often owns many parcels, so an aggregate function is required.  I do the following query to return a list of owners and the summed value of the acres.

SELECT owner, SUM(acres) as totalacres

FROM owner_types

WHERE county ='Monterey County' and detailed_ownership_type = 'company'
GROUP BY owner
ORDER BY totalacres DESC LIMIT 20

Now, I want to create a map that visualizes the location of all of the polygons associated with the owners generated by the above query.  I have tried a subquery like such:

       SELECT * FROM owner_types
       WHERE owner =
           (SELECT owner, SUM(acres) as totalacres

FROM owner_types

WHERE county ='Monterey County' and detailed_ownership_type = 'company'
GROUP BY owner
ORDER BY totalacres DESC LIMIT 20)

But I am violating some subquery rules including more than one column and aggregate functions within subqueries.  Can someone guide me to the best way to match the owner values of the first query with their original geometries to make a new map?


texts: Negative OSM ID values - building polygons

texts:
The use case is simple: extract all the buildings from OSM via osm2pgsql.

In the process, it seems other entities make it into the polygons list, such as "islands", only because they have the "building=NULL" tag available, such as:

I am reading this documentation: http://wiki.openstreetmap.org/wiki/Osm2pgsql/schema

This table contains all polygons (closed ways) which were imported. For polygons which are members in one or more relations, multiple rows will be created: one with the polygon's own ID and tags, and one more for each relation membership. These additional rows contain the tags of the relation and the negative ID of the relation in the osm_id column (since a single relation may contain multiple polygons, negative IDs are not necessarily unique).

I've created a style file which allows me to extract all the building polygons, but it seems some areas, which have the tag 'building=NULL' are not actually buildings, but describing an island in the water, or the actual country itself (which also comes with the 'building=NULL' tag).

Would it be safe to simply discard all the polygons which have a negative ID?

At the moment, I collected all the OSM IDs of these "islands" and other elements tagged wrong, and I'm deleting them via PGSQL:

delete from planet_osm_polygon where  osm_id in (-6040654, -6045552, -6045364, -6188599, -5763894 );

Is there a best practice in relation to this, should I go on the map and feed these fixes in, or is there a way to avoid this type of query where I want to look at building polygons without getting these "odd" entities?


texts: Why is a normalized parameter involving the Green band and Red band not used for Vegetation determination?

texts:
The green region reflects solar energy sufficiently well. So does the NIR region, much more than the green region, indeed. However, reflection in Green region is based on the chlorophyll concentration, to the best of my knowledge. Thus, it is a great indicator of leaf's health and can be used to locate vegetation.

Is it so that the green band is attenuated too much by the atmosphere or something along those lines, for which the green band is not chosen for calculating a normalized parameter? Basically, why is NIR preferred over green region?


texts:
In QGis by pressing shift or ctrl keys a user can select multiple layers in toc.
How can a python script iterate over selected layers?

I'd expect a iface.selectedLayers() but I can find only iface.activeLayer() which return only one layer.


texts:
A project I am working on contains 15 different layers with 5 containing over 50,000 polygons in each layer. Some of these layers must be "queryable" for data (ex: click to get more information of layer). An exported image of the layers, basemap, and client side drawings is a requirement. I am at a crossroads of selecting client side or server side rendering for this project. I need to know if my logic is correct on this? Any libraries that could be used to solve issues I have (Ex: see Click Query or Image Export)? This question focuses on the client side. See Leaflet Server Side rendering for server side perspective.

Client Side Rendering

Layers

To handle this many layers with some larger polygon files we will probably export the layers as topoJSON. Perhaps we could use vector tiles because certain layers contain labels that need to follow polylines (ex: street names) or points that needs to be labeled with the label matching a rotation angle from the JSON.

Chunked Client Side Rendering Layers

We could load the geoJSON that is visible in the browser as the map moves around and remove the geoJSON that is no longer visible using a code like this:

geojson7 = L.geoJSON(null, {
 table: 'bore_lines4326'
}).addTo(map);

geojson8 = L.geoJSON(null, {
 table: 'facility_lines4326'
}).addTo(map);

var overlays = {
    "Bore Lines": geojson7,
    "Facility Lines": geojson8
};

map.on('moveend', function() {
  for(i in overlays){
    if (map.hasLayer(overlays[i]) && overlays[i].options.table){
      findgJSON(overlays[i], printgJSON);
      };
  }
});

function findgJSON(selectedlayer, callback) {
  var dreak = map.getBounds();
  $.ajax({ ***call stuff here***
      },
      success: function (drawpoly){
            callback(selectedlayer, drawpoly);
        }
    });
}

function printgJSON(selectedlayer, drawpoly) {

catdd = [];
for(n in selectedlayer._layers){
   //pid is a primary key
   catdd.push(selectedlayer._layers[n].feature.properties.pid);
 }

 var deletes = catdd.filter(function(o1){
       // filter out (!) items in result2
       return !drawpoly.features.some(function(o2){
           return o1 === o2.properties.pid;          // assumes unique id
       });
  });

  var addpear = drawpoly.features.filter(function(o1){
    // filter out (!) items in result2
     return !catdd.some(function(o2){
     return o1.properties.pid === o2;          // assumes unique id
    });
   });

    selectedlayer.eachLayer(function(layer){
      for(m = 0; m < deletes.length; m++){
       if (layer.feature.properties.pid === deletes[m]){
         selectedlayer.removeLayer(layer);
       }
      }
    });

    selectedlayer.addData(addpear);
};

Click Query

The click function will have to be a point with a buffer around it to grab anything close by the clicked point. Have you ever tried clicking a polyline in Leaflet especially when it is thin ... :/. Also this will allow us to grab overlapping layers. ex: https://github.com/openplans/Leaflet.FeatureSelect or https://github.com/mapbox/leaflet-pip or maybe turfjs.org

Image Export

To export an image client side we would need to use a library like github.com/niklasvh/html2canvas or github.com/mapbox/leaflet-image

*Apparently there has been a bug introduced by the latest version of leaflet that causes the exported image to be misaligned if the map is moved from its original loaded state github.com/mapbox/leaflet-image/issues/66 ... if i remember correctly the same issue occurs with html2canvas also. Away around this issue would be reload the map in another window and fire off the javascript to create the image once the canvas or geoJSON has fully rendered. How to determine this or perhaps a solution to the problem? The problem might lie with tranform: tranlate3d(xx px, xx px, xx px) style of leaflet


texts: I'm wondering if there have been any studies which attempt to analyze or characterize GPS error distribution over different situations/conditions (e.g. surrounded by tall buildings/mountains, precipitation/clouds/fog, rural areas, etc.) to understand how these errors may change as compared to 'perfect' conditions.


texts: Why Autocad map consumes so slow WMS comparing with Arcgis or QGIS, and what can be done about that?

texts:
I created a WMS for aerial photo using Geoserver. Trying to use it with Autocad Map 3d 2013-2014-2015

Why this product consume so bad and slow this and every other WMS?

While the bing aerial used by the same products its very fast. Can be done something about this, or its just the product issue?


texts:
I am trying to add great circle routes between various regions in R. Here's the sample data:

structure(list(CommencingRegion = c("RedSea", "EastAfrica", "GulfofMexico",
"FarEast", "RedSea", "GulfofMexico"), LoadRegion = c("RedSea",
"RedSea", "GulfofMexico", "FarEast", "RedSea", "GulfofMexico"
), DischargeRegion = c("NorthWestAfrica", "WestMedditerranean",
"WestCoastLatinAmerica", "AustraliaNewZealand", "WestMedditerranean",
"WestCoastCentralAmerica"), Count = c(1L, 1L, 2L, 1L, 2L, 5L),
    AvgTCE = c(38879.53, 31783.55, 28520.79, 26068.8, 26054.28,
    25883.81), CLon = c(37.8274879335485, 47.0791103099334, -90.9633701509553,
    146.2727573458, 37.8274879335485, -90.9633701509553), CLat = c(21.4460561443517,
    -12.9570828789565, 25.2035802054683, 47.6530892619773, 21.4460561443517,
    25.2035802054683), LLon = c(37.8274879335485, 37.8274879335485,
    -90.9633701509553, 146.2727573458, 37.8274879335485, -90.9633701509553
    ), LLat = c(21.4460561443517, 21.4460561443517, 25.2035802054683,
    47.6530892619773, 21.4460561443517, 25.2035802054683), DLon = c(-17.1597117430475,
    7.03639948506481, -73.4238157230412, 151.051220297802, 7.03639948506481,
    255.83509305644), DLat = c(24.2308740597312, 38.8907379374372,
    -25.8934046406896, -25.1880219406131, 38.8907379374372, 21.8130318388702
    )), row.names = c(NA, -6L), class = c("tbl_df", "tbl", "data.frame"
))

I would like to plot great circle routes between CommencingRegion to LoadingRegion and then from LoadingRegion to DischargeRegion for every row. Additionally, every row needs to be in a different color and the thickness needs to be proportional to AvgTCE. The last 6 variables in the data above are the coordinates in Lat Long for the commencing, loading and discharging regions respectively. I am quite clueless on how to go about achieving this. This is what I have tried and failed:

leaflet() %>%
  addTiles() %>%
  for(i in 1:6){
    addPolylines(data=gcIntermediate(c(ByRoute$CLon[i],ByRoute$CLat[i]),c(ByRoute$LLon[i],ByRoute$CLat[i]),n=100,addStartEnd = T,sp=T))
  }


texts: qgis2web pop up image -- path definition

texts:
I have consulted the post on GIS stack exchange and read the docs, however I cannot find a solution that works.

The idea is to put a .png image of a site report within the maplayer for web output.  The points show up fine but the popup image does not.

The column that holds the file name for the image has been formatted in the following ways (plus many more varieties):

"filename.png" #<---------- images in the same folder as the project
"images/filename.png"#<-------- images in an images "folder" in the
    samefolder as the project
"../images/filename.png"

Mulitple variations on this have been tried by loading the different methods in the attribute table and having as pop up value(output from most recent try):

The field attributes are set to photo:

I have taken a look at the output generated by the script in the Data folder that holds the css, js and other files associated with the webpage export(in my case these files go to 'users>me>AppData>Local>Temp>qgis2web>folder'

The data is in a dictionary of dicts I put the first one here:

{ "type": "Feature", "properties": { "field_1": "Aare_Bern_CaveltiN", "Density": 1.61, "lon_x1": 7.4733190977753496, "lat_y1": 46.923578788505381, "City": "Muri bei Bern", "File": "Aare_Bern_CaveltiN.png", "FileName": "\"\\images\\Aare_Bern_CaveltiN.png\"", "Source": "", "File2": "\"Aare_Bern_CaveltiN.png\"", "FileNameThree": "\"images\\Aare_Bern_CaveltiN.png\"", "FileNameFour": "", "FileNameFive": "\"\"" }, "geometry": { "type": "Point", "coordinates": [ 7.47331909777535, 46.923578788505381 ] } }

I cannot figure out exactly what needs to be passed on to the plugin for this to work.
Nor are the examples I have found clearing this up.

Is it an absolute path or relative path?
Path to where? Although the project is one folder the output goes to temp.

Or maybe even more precise: When I call the column value from python (before exporting csv for the map layer) like this: print(df['filename'][0])

where df=layer data and filename = column that holds the files.

Should I see this: "filename.png" (in quotes) or this: filename.png (no quotes)?


texts:
I'm trying to create a new surface in Autocad Civil 2012. I have elevations for a landfill, and it seems like I'm following all the correct steps. But when I check the surface out under Surface Properties/Statistics, it appears the X,Y & Z coordinates have been scaled down by 3 decimal places. I'm using UTM coordinates to create the points.

Example:

Min elevation = 0.071m
Min X = 481.510m
etc.

Any help greatly appreciated!

-jason m.


texts: Remove Errant LAS Points Outside Tile Boundary

texts:
I have a series of LAS tiles for a project area. One of those tiles contains points outside of the tile area (out to 0,0 in the coordinate system). When we convert the las tile to a raster the result is insanely large, in the multi-GB range, when it should be only a few MB.
Is there a way to filter the errant LAS points to only those in a specified rectangular area, and create a new LAS tile as a result?


texts: InitSpatialMetaData errors with Spatialite

texts:
This Stackexchange may not be the best place for this question but I will start here.

I am trying to create a spatially-enabled, sqlite3 database.  The following code is a short example of what I am trying to do.

#include <sqlite3.h>
#include <iostream>

//extern InitSpatialMetaData; // Always is an incomplete type no matter which spatialite headers I include.

int main()
{
  sqlite3* db;
  char* errMsg;
  int rc {0};

  rc = sqlite3_open_v2(":memory:", &db, SQLITE_OPEN_READWRITE|SQLITE_OPEN_CREATE, nullptr);
  rc = sqlite3_exec(db, "SELECT InitSpatialMetaData();", nullptr, 0, &errMsg);
  if (rc != SQLITE_OK)
  {
    std::cout << "CreateSpatialMetaData error: " << errMsg << std::endl;
    sqlite3_free(errMsg);
  }
  sqlite3_close(db);
  return 0;
}

I have tried compiling with the following:

g++ example.cpp --std=c++11 -lspatialite -lsqlite3

g++ example.cpp --std=c++11 `pkg-config --libs spatialite` -lsplite3

g++ example.cpp --std=c++11 -Wl,--whole-archive -lspatialite -lsqlite3 -Wl,--no-whole-archive

Everything will compile and link.  But at runtime the following error appears:

CreateSpatialMetaData error: no such function: InitSpatialMetaData

How can I include the InitSpatialMetaData symbol?  Where is this even defined?

I tried going through the spatialite_gui source to see how this is implemented (which is not much different than mine).  For Alessandro's work everything is in three files so tracking dependences is a bit difficult for me as there so much stuff in-scope.

I am building on Debian stable with supported packages:

g++ 4.9.2
libspatialite-dev 4.1.1
libsqlite-dev 3.8.7.1

I am also unable to start the sqlite3 plugin at runtime, this results in an unauthorised error message.


texts: Georeferencing declassified CORONA images of coastlines?

texts:
I'm trying to georeference declassified CORONA images of coastlines in order to asses coastline changes and human impact.

In most case, georeferencing is easy, but since images of coastlines can only have GCP's inland of the object of interest (in rare occasions, islands an artificial features seaward), namely, shore line and strip of coastal sands, there will be inherent georeferencing errors.

How can I minimize errors?

I've thought of two options, but have yet to try them:

Georeferencing long swathes of coastline at a time (roughly 2Km long by 200m wide) with a linear/1st order equation (Corona images aren't orthorectified), so any localized errors will have minute effect on the outcome. possible problems: not enough GCP in 200 meters wide coastal area.
Short swathes of coastline (roughly 5oom * 1Km), thus probably increasing the potential GCP's inland and reducing local errors. possible problems: cumulative error if I try to use one image as GCP for another.

The reason I am inclined to georeference one complete CORONA strip is to save space, and ease workload for my computer.

I am using ArcMap 10.2, and have a recent orthorectified high resolution image of the terrain


texts: Identifying the nearest center of a buffer for an observation when buffers overlap

texts:
I am using ArcGIS 10.5. Desktop version and would like to do the following:

I have a 150 points with coordinates over a map of a city and would like to create 300m buffers around each of them. I also have coordinates for several thousands of observations spanned across the city. After creating the buffers I would like to compare observations inside the buffer zones to those outside of the zone but still within the same zip code area.

My problem is this:

Many of the buffers are overlapping so how do I attach each observation to the nearest center point of a buffer it belongs to?

If possible I would also like to count the amount of buffers each observations is in (ranging from zero to several).

The ideal form of my data would therefore indicate:

whether an observation is within a buffer zone?
which buffer center point is the nearest to each observation that is within a buffer?
how many buffers overlap each observation point?


texts: How to set the value of a pixel based on its current value?

texts:
I am attempting to set the value of a pixel based on its current value.  I am using the Raster Calculator in Quantum GIS dev 1.7 trunk.

For example, if the value is 0, I want to set it to 255.

Based on an example that I found, something like this should work, but it doesn't.

myRaster@1=0 255

The RasterCalc plugin appears to use a function syntax, but that doesn't appear to be valid either

eq( myRaster@1, 0, 255)

Any suggestions?  Thanks.


texts: Computation Time Out in Complex Multi-Band Normalization of Sentinel 2 Image Composite

texts:
I am trying to implement Pseudo-Invariant Feature Normalization for Sentinel 2 image in 2015 based on the median composite acquired in 2019 (script can be found here). I once ran the script without an issue, but this time, the script returned the “User memory limit” error. Usually exporting the result to Google Drive would solve the issue, but not this time. I would suspect that this chunk below (line 41-62 in the code) is the culprit:
var normalized_image = medImage_2015.bandNames().map(function(name){
  name = ee.String(name);
  var band = medImage_2015.select(name);
  var bandMask = band.updateMask(PIFmask);
  // var bandRef = medImage_2019.select(name.cat('_median')).updateMask(PIFmask);
  var bandRef = medImage_2019.select(name).updateMask(PIFmask);
  var bandReg = bandMask.addBands(bandRef);
  bandReg = bandReg.reduceRegion({
    reducer: ee.Reducer.linearFit(),
    geometry: studyBound_poly,
    scale: 10,
    bestEffort: true
  });
  return band.expression(
    '(x * slope) + intercept', {
      'x': band,
      'slope': ee.Number(bandReg.get('scale')),
      'intercept': ee.Number(bandReg.get('offset'))});
});
normalized_image = ee.ImageCollection(normalized_image).toBands().rename(medImage_2015.bandNames());

Can you please share some tips to overcome this issue?


texts: Processing toolbox: random points inside polygon Parameter meaning

texts:
I am using the random points inside polygon Tool from the QGIS-processing toolbox. Unfortunately I dont know want the following parameter Exponent means.
Also I would like to know a little bit about the optional usage of a transformation curve. Obviosly in "default usage" (unticked) there is a linear relationship. Any cases where the usage of applying a transformation curve is useful?
The documentation https://docs.qgis.org/3.10/en/docs/user_manual/processing_algs/qgis/vectorcreation.html#random-points-inside-polygons does not provide any hints to EXPONENT and the usage of the TRANSFORMATION CURVE.


texts: Extracting Soil moisture value using Google Earth Engine?

texts:
I'm writing a script to extract the soil moisture content of a specific polygon. The polygon covers multiple pixels, I need a weighted average of the values corresponding to those pixels. When I execute the script, the output is 'null'. I'm relatively new to Google earth engine, and haven't been able to figure out what I've done wrong.

// which dataset + where + when
var image = ee.ImageCollection('NASA_USDA/HSL/SMAP_soil_moisture')
                .filterBounds(Zutendaal)
                .filterDate('2018-03-05','2019-03-06');

// select bands -> ssm = surface soil moisture (mm)
var soilmoisture = image.select('ssm');

// Take average
var meanDict = image.reduce(ee.Reducer.mean());

// Get the mean from the dictionary and print it.
var mean = meanDict.get('ssm');
print('Gemiddelde bodemvochtgehalte (mm)', mean);

Could someone help me?

Here's the link: https://code.earthengine.google.com/c37820b890bec18142baa2b475dd7dd2


texts: QGIS OpenGeo Explorer plugin and https geoserver

texts:
Is there a way to connect to a secure geoserver with OpenGeo Explorer plugin in QGIS?  When I connect to geosever with a secure protocol I get the error Warning: Cannot connect using the provided certificate/key values  and the site does not load in the catalog list and when I try to connect via the ip or hostname on my vpn with x.x.x.x/geoserver or hostname/geoserver it connects and I see it in the catalog  but I get an error in QGIS.

Traceback (most recent call last):
  File "C:/Users/xxx/.qgis2/python/plugins\geoserverexplorer\gui\explorer.py", line 101, in run
    command(*params)
  File "C:/Users/xxx/.qgis2/python/plugins\geoserverexplorer\gui\gsexploreritems.py", line 345, in populate
    layerItem.populate()
  File "C:/Users/xxx/.qgis2/python/plugins\geoserverexplorer\gui\gsexploreritems.py", line 714, in populate
    if layer.default_style is not None:
  File "C:\Users\xxx\.qgis2\python\plugins\geoserverexplorer\ext-libs\gsconfig-1.0.6-py2.7.egg\geoserver\layer.py", line 127, in _get_default_style
    return self._resolve_style(element) if element is not None else None
  File "C:\Users\xxx\.qgis2\python\plugins\geoserverexplorer\ext-libs\gsconfig-1.0.6-py2.7.egg\geoserver\layer.py", line 135, in _resolve_style
    return self.catalog.get_style_by_url(style_workspace_url)
  File "C:\Users\xxx\.qgis2\python\plugins\geoserverexplorer\ext-libs\gsconfig-1.0.6-py2.7.egg\geoserver\catalog.py", line 797, in get_style_by_url
    dom = self.get_xml(style_workspace_url)
  File "C:\Users\xxx\.qgis2\python\plugins\geoserverexplorer\ext-libs\gsconfig-1.0.6-py2.7.egg\geoserver\catalog.py", line 211, in get_xml
    response, content = self.http.request(rest_url)
  File "C:/Users/xxx/.qgis2/python/plugins\geoserverexplorer\geoserver\retry.py", line 19, in decorator
    raise e
SSLError: [Errno 185090050] _ssl.c:340: error:0B084002:x509 certificate routines:X509_load_cert_crl_file:system lib

How do I get QGIS and the plugin to accept a cert?


texts: Create points symbology with more than one parameter

texts:
In QGIS i would like to create symbology for points layer where the symbology will take into account more than one parameter (currently two).
i.e.:
My points layer has street names and house numbers and i would like to have each street with a different symbol and each number with a different color. Each street_name+house_number is unique. I wouldn't mind if distant streets share symbols.

What will be the best way to achieve it in QGIS?
Consider that Rule-based symbology is manual and i have many streets, i don't see it as a good solution for this specific problem.


texts: WMS and WFS access in ArcGIS Pro localhost

texts:
I've set up a localhost GeoServer instance and connected without issue to QGIS.  I've tried to connect to the localhost instance in ArcGIS Pro and get errors (Resource not found 404 error).  Is there a trick to getting the software to recognize a GeoServer instance that is local?   Does ArcGIS Pro only allow connections to public facing sites?
I've been trying to connect using this URL:
http://localhost:8080/geoserver/<*workspace*>/wms?

Is there something wrong here? The workspace I input exists and I've logged in as GeoServer admin (full access to the workspace).  If I mess with the URL I can get it to connect but it doesn't show any of the existing default services.


texts:
We are working with GeoServer 2.16 version, and we have installed the AWS S3 Bucket plugin, to load GeoTIFF from our public S3 bucket.
We are working with GeoServer REST API to load a GeoTIFF from S3 on a specific workspace using POST request, this is the XML parsed:

<coverageStore>
     <name>test_dem</name>
     <workspace>Sentinel_1</workspace>
     <enabled>true</enabled>
     <type>S3GeoTiff</type>
     <url>s3://xxxxxxxxx/public/test/geoserver/xxxxxxxxx_cog.tif?useAnon=true&awsRegion=EU_CENTRAL_1</url>
</coverageStore>

and this is the endpoint of REST Request:

http://xxxxxxxxxxx:8080/geoserver/rest/workspaces/Sentinel_1/coveragestores?configure=all

This request gives us this error:

---- Debugging information ----
cause-exception     : com.thoughtworks.xstream.io.StreamException
cause-message       :
class               : java.lang.String
required-type       : java.lang.String
converter-type      : com.thoughtworks.xstream.converters.SingleValueConverterWrapper
wrapped-converter   : com.thoughtworks.xstream.converters.basic.StringConverter
line number         : 6
class[1]            : org.geoserver.catalog.impl.CoverageStoreInfoImpl
required-type[1]    : org.geoserver.catalog.impl.CoverageStoreInfoImpl
converter-type[1]   : org.geoserver.config.util.XStreamPersister$StoreInfoConverter
version             : 2.16-SNAPSHOT
-------------------------------

And this is the GeoServer log:

Caused by: com.thoughtworks.xstream.io.StreamException:
at com.thoughtworks.xstream.io.xml.XppReader.pullNextEvent(XppReader.java:124)
at com.thoughtworks.xstream.io.xml.AbstractPullReader.readRealEvent(AbstractPullReader.java:148)
at com.thoughtworks.xstream.io.xml.AbstractPullReader.readEvent(AbstractPullReader.java:135)
at com.thoughtworks.xstream.io.xml.AbstractPullReader.getValue(AbstractPullReader.java:180)
at com.thoughtworks.xstream.converters.SingleValueConverterWrapper.unmarshal(SingleValueConverterWrapper.java:49)
at com.thoughtworks.xstream.core.TreeUnmarshaller.convert(TreeUnmarshaller.java:72)
... 116 more
Caused by: org.xmlpull.v1.XmlPullParserException: entity reference name can not contain character ='
(position: START_TAG seen ...c/test/geoserver/eu_dem_v11_E00N20_cog.tif?useAnon=true&awsRegion=... @6:106)
at org.xmlpull.mxp1.MXParser.parseEntityRef(MXParser.java:2222)
at org.xmlpull.mxp1.MXParser.nextImpl(MXParser.java:1275)
at org.xmlpull.mxp1.MXParser.next(MXParser.java:1093)
at com.thoughtworks.xstream.io.xml.XppReader.pullNextEvent(XppReader.java:109)
... 121 more

Any idea about this error?
We have tried to load the same S3 path directly from GeoServer UI Interface and in this case, works.
Any idea?


texts: ST_Intersects throws an error on mixed SRIDs that && permits

texts:
I found that there is a difference in result between using A && B operator or one of the geo functions like ST_Intersects(A,B)

While the ST_Intersects functions return query results without errors the operator version SOMETIMES results in errors over the SRID:

  ERROR:  Operation on two geometries with different SRID's.

  ERROR: Operation on two geometries with different SRIDs
  SQL state: XX000

All of the records seem have the correct SRID, none missing.

The query is generated by MAPSERV WMS server, so what is there to do?

Sample query:

Operator &&

select encode(AsBinary(force_collection(force_2d("geometry")),'NDR'),'hex') as geom, "id"
from (select * from soils.vw_zones) as vw
where geometry && GeomFromText('POLYGON((75722.1945223652   437331.342330005,75722.1945223652 503295.946741193,169438.897876772 503295.946741193,169438.897876772 437331.342330005,75722.1945223652 437331.342330005))',find_srid('','soils.vw_zones','geometry'))

ST_Intersects

select encode(AsBinary(force_collection(force_2d("geometry")),'NDR'),'hex') as geom, "id"
from (select * from soils.vw_zones) as vw
where st_intersects(geometry,GeomFromText('POLYGON((75722.1945223652 437331.342330005,75722.1945223652 503295.946741193,169438.897876772 503295.946741193,169438.897876772 437331.342330005,75722.1945223652 437331.342330005))',find_srid('','soils.vw_zones','geometry')))


texts: How to combine the content of several mapsets into one new mapset in GRASS.?

texts: I have some mapsets and i want to assemble them into one new mapset and have everything into it. How do i do it?


texts: Setting Merge Rules in Spatial Join for 90+ fields or changing default merge rule?

texts:
I am trying to perform a spatial join that will aggregate about 90 fields using a SUM function.

In the Field Map I get the option to right click on each field and set the merge rule to SUM.

This would be extremely tedious for me to do, especially since I have to perform such a spatial join for 30+ layers.

Is there a way to change the "default" merge rule option to SUM instead of FIRST? How about change the field merge rules in bulk?

Screenshot to illustrate the point


texts: When I import polygon kml layer, how do I change projection and offset to correct for GMT wrap?

texts:
I created a polygon in Google Earth and saved to kml. After importing as a vector layer, where the polygon crosses the prime meridian, the polygon wraps around the world.

So far, I have figured out that I need to change the projection [and offset?], but I cannot figure out how to do that.

I am new to QGIS (and GIS in general). I have already spent many hours trying to figure this out. Can somebody please help?


texts: How to determine irrigable acres from a single source/point?

texts:
I am trying to graphically depict the exact amount of irrigable acres from a single point, in this case a headgate from a drainage ditch. Basically, I have the original water right decrees (i.e. confusing legal documents) from the late 1800's and am interested in accurately illustrating the legally allowable area to which they can irrigate.

Example: From the stream a drainage ditch runs west for 1 mile (on contour) before it opens into a small pasture, irrigating 50 acres below its outlet.

Is there a way to show this acreage below this specific outlet point? Perhaps a sort of reverse Snap Pour Point operation? Some function that fills in a specified area below a certain point, following topography lines?

:


texts: How to create a polygon from Excel coordinate file

texts:
Currently I have an Excel file with 2 columns.

area number.
coordinates.

The coordinate column contains several X and Y points those are seperated with a comma.
Next an example how the points are noted:

5.72567034 -3.61719203,5.72567034 -3.61719203,5.72565126 -3.61720133,5.72562551 -3.61721635,5.72560167 -3.61723351,5.72557211 -3.61725235, etc.

point 1: 5.72567034 (longitude) & -3.61719203 (longitude) {,}

point 2: 5.72567034 -3.61719203,

point 3: 5.72565126 -3.61720133,

I want to extract this data to ARCMap 10.6 and make polygon shapefiles out of it.
The idea is to extract a polygon from points (points to polygon). but I am not able to create points out of these file. Because all the points of one area are in the same cell.
When I try to put every point in one cell by doing in Excel "text to column" it exceeds the width of the worksheet (meaning >256 points per area). And thus by doing this I would lose some data.
Then there is this problem that there are over 1000 areas of which this process should be done so manually is not an option.

Does anyone have a solution?


texts: ArcGIS Server Error: "Feature service requires a registered database"

texts:
I'm attempting to publish a feature service to ArcGIS server. I'm running SQL Server Express and ArcSDE. When I try to publish the service I get the error "00090: Feature service requires a registered database".

When I validate the database registration it appears to succeed, but I still can't publish the service. The features I am publishing are located inside of the geodatabase that I am registering. This geodatabase is located on the same server that ArcGIS Server resides, but I am attempting to publish from another machine. Am I missing some crucial step?

Here is a bigger version of the screenshot.


texts: Restrict pan and zoom in ArcGIS Javascript API?

texts:
I'm creating a web map of a university campus.

I have the map functioning, but I don't want to be able to pan off of the basemap, into whitespace.

Do you know of a function that will auto-pan the extent back to the basemap?


texts: How to write values from arcpy statistical output to table?

texts:
The attached script is my stalled attempt to perform the following workflow:

Create a table and create fields
Loop through a folder of point shapefiles and perform average
nearest neighbor analysis
Create a row and write "nn_values[]" to the corresponding fields
move to the next shapefile and repeat

The script successfully produces a table, runs the aNN analysis for a shapefile and adds an empty row to the table. As an end-product, I would like a table with all of the values from the statistical analysis. How can I write the associated aNN statistical output (i.e. "nn_values[]") to the corresponding fields for each iteration in the for loop?  Many thanks for the guidance.

# Import system modules
import arcpy, os
from arcpy import env

# Set over write
env.overwriteOutput = 1

# Local variables...
env.workspace = r"C:\temp"
Dir = env.workspace
out_name = "NNtable.dbf"

### Execute CreateTable
arcpy.CreateTable_management(Dir, out_name)

### Add fields
arcpy.AddField_management(Dir + "\\" + out_name, "index", "FLOAT")
arcpy.AddField_management(Dir + "\\" + out_name, "z_score", "FLOAT")
arcpy.AddField_management(Dir + "\\" + out_name, "p_value", "FLOAT")
arcpy.AddField_management(Dir + "\\" + out_name, "expected", "FLOAT")
arcpy.AddField_management(Dir + "\\" + out_name, "observed", "FLOAT")
arcpy.AddField_management(Dir + "\\" + out_name, "pathway", "TEXT")

# List FCs

fclist = arcpy.ListFeatureClasses()

for fc in fclist:
    try:

        # Obtain Nearest Neighbor Ratio and z-score
        # Process: Average Nearest Neighbor...
        nn_output = arcpy.AverageNearestNeighbor_stats(fc, "EUCLIDEAN_DISTANCE", "NO_REPORT", "#")

        # Create list of Average Nearest Neighbor output values by splitting the result object
        nn_values = nn_output.split(";")
        print "The nearest neighbor index is: " + nn_values[0]
        print "The z-score of the nearest neighbor index is: " + nn_values[1]
        print "The p-value of the nearest neighbor index is: " + nn_values[2]
        print "The expected mean distance is: " + nn_values[3]
        print "The observed mean distance is: " + nn_values[4]
        print "The path of the HTML report: " + nn_values[5]

    except:
        # If an error occurred when running the tool, print out the error message.
        print arcpy.GetMessages()

    # Create insert cursor for table
    #
    rows = arcpy.InsertCursor(Dir + "\\" + out_name)
    row = rows.newRow()
    rows.insertRow(row)

    # Delete cursor and row objects to remove locks on the data
    del row
    del rows


texts: Calculate the annual mean NDVI for a region using the Maximum Value Composite

texts:
I first used the Maximum Value Composite to generate a annual NDVI for a watershed, and then calculated the mena NDVI  for the entire basin. I don't know if this code is correct. In addition, I want to get the annual mean NDVI  of the basin in 2013-2018 (after Maximum Value Composite), I don't know how to modify my code.

var roi = ee.FeatureCollection("users/qq2385973028/shape/Chabagou_catchment");
var dataset = ee.ImageCollection('MODIS/006/MOD13Q1')

function generateImage() {
  Map.centerObject(roi, 10);
  var img = dataset.filterDate("2017-1-1", "2017-12-31")
              .filterBounds(roi)
              .select("NDVI")
              .mosaic()
              .clip(roi);
  var visParam = {
   min: 0.0,
  max: 8000.0,
    palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
      '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
  };
  Map.addLayer(img, visParam, "NDVI");
  return img;
}
var img = generateImage();

 //image collection
  var l8Col = dataset.filterDate("2015-1-1", "2015-12-31")
              .filterBounds(roi)
              .select("NDVI")
  var image = l8Col.reduce(ee.Reducer.max());
  print(image);

  var visParam = {
     min: 0.0,
  max: 8000.0,
    palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
      '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
  };
  Map.addLayer(image.clip(roi), visParam, "NDVI_min1");

 var image2 =l8Col.max()
 Map.addLayer(image2.clip(roi), visParam, "NDVI_min2");
  //image
  var dict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: roi,
    scale: 250,
    maxPixels: 1e13
  });
  print("dict", dict);


texts:
In my C# Console Application, Visual Studio 2010 highlights new FeatureClassNameClass() in my below line of code and states "ESRI.ArcGIS.Geodatabase.FeatureClassNameClass has No constructors defined" and then states cannot be embedded. Use aplicable interface instead.  How else can I define pOutFeatClassName as type IFeatureClassName?

IFeatureClassName pOutFeatClassName = new FeatureClassNameClass();

If I change the code to below, then Visual Studio does not throw an immediate error but once I run my application Visual Studio immediately bombs at this line and only states "error vshost32.exe".  What is going on?  Many Thanks for any suggestions..

IFeatureClassName pOutFeatClassName = (IFeatureClassName)new FeatureClassName();


texts: I'm have a huge shapefile and I'm trying to delete large portions of it. I can't seems to do it in the edit mode. Can someone please walk me through this process? I'm using QGIS 1.8.


texts: Programmatically creating ArcGIS Image Service and adding it to map using C# / .NET?

texts:
How can I programmatically create ArcGIS Image Service and add that service to my MXD file programmatically?

I am using C# .Net.


texts: Suitable algorithms for simplifying isochrone polygons

texts:
I'm looking for an algorithm for simplifying isochrones which were generated using OSM data.

These are generally characterised as having lots of small "holes" mostly due to missing OSM data but also frequently having small "islands" due to the use of transport like subways. I have tried a number of tools for this but all are too aggressive at shrinking space - most critically, deleting the islands. Reachability to these public transport islands is very important.

So rather than focusing on tools, I wanted to step back a bit and ask if any particular algorithms are well or not well suited to such a task? i.e. being aggressive at removing holes but conservative about not removing islands.

The isochrones start life in GeoJSON format. They eventually will be imported to PostGIS.


texts:
How to check if a file is a valid raster image? Valid in the sense of not being corrupt, the file record is complete. The GIS corollary to "zip --test somefile.zip" to test if an archive's integrity is intact.

An example of where this would be useful: have just batch processed 500 geotiff files, losslessly compressing and adding overviews, and would like to verify resultant files are not broken before removing originals.

I'm partial to gdal in this context, because it works better on the command line, but Arcgis solutions are welcome too.


texts: Splitting raster into smaller chunks using GDAL?

texts: I have a raster (USGS DEM actually) and I need to split it up into smaller chunks like the image below shows. That was accomplished in ArcGIS 10.0 using the Split Raster tool. I would like a FOSS method to do this. I've looked at GDAL, thinking surely it would do it (somehow with gdal_translate), but can't find anything. Ultimately, I'd like to be able to take the raster and say how large (4KM by 4KM chunks) I would like it split up into.


texts:
I'm developing a geospatial db for a company that operates at a number of large facilities. Its primary purpose is to help staff locate themselves and equipment within each facilities' multi-storey buildings. I'm using ArcGIS's Building Interior Space Data Model (BISDM) as a schema template and implementing it in PostGIS. The table structure looks something like this:

facility
  building
    floor
      boundary
      interior_space
        equipment

  transport
    road
    pedestrian_path

  perimeter
    fence
    vehicle_gate
    pedestrian_gate

Each parent table has a one-to-many relationship with it's children, so the schema defines a simple hierarchy. Each table has short_name, long_name and geom columns (as well as others, depending on the table).

My problem is this: I want to build a reverse-geocoder API that takes lat/long as input and returns a string representing the address of the matching locations. I've gotten stuck on how to construct this string.

My idea was to create a query that takes as its starting point the row containing the geom closest to the lat/long point and then traverse up the table hierarchy, collecting the value of short_name for each parent record until it reaches facility. So the assembled address in this case might be 'Room 323A, Floor 4, Building 12, Facility Q' (interior_space->floor->building->facility). However, because I don't know which table my query will start from (since it could be a point anywhere within the facility) I don't know what the parent table names will be in advance, and so I can't construct any JOIN queries.

I'm familiar with approaches like adjacency lists and closure tables to model hierarchal systems, but these all require the data to be in a single table. The PostGIS requirement (well, strong recommendation) for different geometry types to be kept in separate tables means I can't use these approaches.

Can anyone suggest how I can construct this string without having to create a bunch of hard-coded SQL queries?


texts: How best to share maps created in MapInfo Online

texts:
Does anyone know of a good free or low cost option for sharing Maps produced in mapinfo online?

I would want up to 5 users to access the maps with a password and be able to turn layers / points on and off and zoom in and out.

Basically I am looking to upload the Mapinfo files or Shapefiles and create thematic maps (ideally multivariate) online.

I have looked at Mapinfo Stratus but costs are prohibitive for a small company (<5 people)!

Many thanks for any help

Kev


texts: Extract information from a layer underneath of a buffer

texts:
I have a line layer (bus route) and I created an one-mile buffer around it.
Now I want to extract information (population/employment) from the census layer underneath of the buffer area only.
How can I do that?


texts: Is there a size restriction on KML layer in esri map?

texts: My web site utilizes ESRI's map and it contains several KML layers. Most of them show up fine except for the ones bigger than ~4MB. Are there any restrictions on KML file size?


texts: Arcobjects- IAnnotateLayerPropertiesCollection2.QueryItem in Python

texts:
Anyone tried access label properties trough ArcObjects with Python?

I try to obtain reference to label object with IAnnotateLayerPropertiesCollection2.QueryItem, but when I do:

pA = pGLyr.AnnotationProperties
new = NewObj(esriCarto.LabelEngineLayerProperties,esriCarto.IAnnotateLayerProperties)
pA.QueryItem(0,new)

I get:
Runtime error : call takes exactly 2 arguments (3 given)

I'm not sure where is the error comming from.

Any ideas?


texts: Creating forest mask using PALSAR Forest/Non-Forest Map for Sentinel 2 image collection

texts:
I trying mask forest areas in Sentinel- 2 image collection over my study area using the 25m resolution global PALSAR Forest/Non-Forest Map. I followed the masking example given in GEE guides that uses Hansen et al. forest change dataset. It shows fnf.eq is not a function as an error when I am defining the mask layer. The code is attached below. How can I solve the problem?
var geometry:https://code.earthengine.google.com/?asset=users/monish_d/MP
var s2= ee.ImageCollection("COPERNICUS/S2_SR").filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10));
var vis = {bands: ['B4', 'B3', 'B2'], max: 2000, gamma: 1.5};
var area = ee.FeatureCollection(geometry); var opacity = 0.5; Map.centerObject(area, 6.5);
var ImCol = ee.ImageCollection(s2.filterDate('2020-02-15','2020-02-28') //// Filter by dates.
.filterBounds(area));
//// Creat a Cloud Mask.
function maskS2sr(image) {
 var cloudBitMask = ee.Number(2).pow(10).int(); ////cloud band
 var cirrusBitMask = ee.Number(2).pow(11).int();//// cirrus band
 var qa = image.select('QA60'); //// Get the pixel QA band.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0)  ////  Flags set to zero for clear conditions.
.and(qa.bitwiseAnd(cirrusBitMask).eq(0)); ////  Flags set to zero for clear conditions. return image.updateMask(mask) //// Return the masked image, scaled to TOA reflectance, without the QA bands.
  .copyProperties(image, ["system:time_start"]);}
//// Apply cloud mask to image collections
var CM_ImCol = ImCol.map(maskS2sr);
////creat forest mask
var dataset = ee.ImageCollection('JAXA/ALOS/PALSAR/YEARLY/FNF')
              .filterDate('2017-01-01', '2017-12-31'); ////global 25m resolution PALSAR-2/PALSAR SAR mosaic
var fnf = dataset.select('fnf');
var Fmask = fnf.eq(1); ////1:forest, 2:nonforest, 3:water
////Define forest Mask function
var maskforest = function(image) {return image.updateMask(image.select(Fmask));};
/////Apply forest mask to image collections
var CFM_ImCol = CM_ImCol.map(maskforest);
////Mosaic and clip images
var mos = CFM_ImCol.median().clip(area);
Map.addLayer(mos, vis,'Masked image');


texts: I am running ArcGIS Server 10, with ArcSDE 10. I have a geodatabase that I would like my colleague to connect to over the LAN at our office. How do I go about doing this? I would prefer a list of simple steps than a link to the ESRI Web Help URL.


texts: Extracting raster values and creating new raster in QGIS?

texts:
I want to extract certain values of a raster image in order to create a new raster which then should only contain the extracted/chosen values.

How to do this in QGIS?


texts: Difference between ImageCollection and FeatureCollection

texts: I wanted to know the difference between FeatureCollection and ImageCollection in Google Earth Engine. The documentation does not mention the difference, it just says that many methods such as maps, reduce can be applied to both and that both return a list


texts: Reducing number of decimal places after X/Y coordinates in a FC

texts: I'm using arcpy.AddXY_management to add coordinates, and I want to reduce each X/Y field by x number. I've tried to reduce it by string slicing and that obviously didn't work. I also tried setting the precision but that doesn't work on FC's. Anyone have a better idea?


texts: Obtaining real world long/lat from GeoTIFF file using GDAL and C#?

texts:
I am trying to obtain and convert latitude and longitude values into British National Grid coordinates from GeoTIFF using GDAL C# file, I managed to read the coordinate system and corner coordinates as shown in the below picture.
The thing is, since that GDAL is more of Python documentation and I can't go through it, I am not sure what is the next step or how to do it using SpatialReference() and CoordinateTransformation().
This is the code I tried so far but still getting the wrong numbers, I think its because they are  projected coordinate system and needed to be converted in British coordinate system.
How can I do this?

GdalConfiguration.ConfigureGdal();
GdalConfiguration.ConfigureOgr();
Gdal.AllRegister();

string fileName = @"path/DSM_D01.tif";
Console.WriteLine(DateTime.Now.ToString());
List<int> resolutionId = new List<int>() { 9 };
Dataset ds = Gdal.Open(fileName, Access.GA_ReadOnly);

double[] gt = new double[6];
ds.GetGeoTransform(gt);

int Rows = ds.RasterYSize;
int Cols = ds.RasterXSize;
Band band = ds.GetRasterBand(1);
double startX = gt[0];
double startY = gt[3];
double interval = gt[1];
double x, y;

for (int k = 0; k < Rows; k++)
{
    y = startY - k * interval; //Current lat
    int[] buf = new int[Cols];
    //ReadRaster parameters are StartCol, StartRow, ColumnsToRead, RowsToRead, BufferToStoreInto, BufferColumns, BufferRows, 0, 0
    band.ReadRaster(0, k, Cols, 1, buf, Cols, 1, 0, 0);
    //iterate each item in one line
    for (int r = 0; r < Cols; r++)
    {
        if (buf[r] != -32768)   //if pixel value is not NoData value
        {
            x = startX + r * interval;  //current lon
            string BritishNationalGrid =

            //Console.WriteLine("Lat: " + y + " Long: " + x);
        }
    }
}
Console.ReadLine();


texts: Adding/updating field to polygon feature class that lists bordering (neighbor) polygons?

texts:
I have a dataset of zipcodes.

I want to add a new field "NEIGHBORS" that includes bordering polygons in a comma delimited manner such as attached image:


texts: Can you determine the row,col of an un-orthorectified image corresponding to a lat,lon?

texts:
Is there are way to determine the row and column of a pixel in an un-orthorectified image that corresponds to a latitude and longitude given a particular DEM? That is, convert a spatial coordinate to a row,col coordinate for a specific RPC file and DEM source?

Use case: I have a stack of satellite images in nitf format and lat,lon coordinates for points of interest. I want to crop the images around each of the POIs and circumvent any pixel distortions due to orthorectification. So rather than orthorectify the images and crop them around a lat,lon, I want to convert the lat,lon to row,col and crop the un-orthorectified image.

Current approach: The only way I can think to do this for a particular image is to:

Create a new 2-band image that is the same size as the original image where the pixel values are the row and column.
Orthorectify this new image
Sample the new image at a lat,lon to get a row,col coordinate.

For my processing I'm using Python and GDAL bindings.

Problem: This approach seems very inefficient with time and memory, and I don't fully understand the accuracy and precision. So if it's possible, I would prefer to essentially reverse-orthorectify a lat,lon coordinate.


texts: Shift between OpenLayers and vector layer QGIS 2.20

texts:
Can you please help me solve this issue with OpenLayers Plugin. I see some shift with a vector layer of Laos water boundaries when the background is Google Physical Layer.  I am using CRS EPSG:32647 (UTM 47N). QGIS version is 2.2.0.

Please see the picture attached.


texts:
I want to draw only one polygon, modify it and then save the polygon. I already do the saving. I'm not able to modify the polygon as it snaps.

Here is what I need but I'm not able to make it work. https://openlayers.org/en/latest/examples/draw-and-modify-features.html

This is my map:

function map() {
  map = new ol.Map({
  controls: ol.control.defaults({
  attributionOptions:  ({
  collapsible: true
  })
}).extend([
  scaleLineControl, new escalaControl(), new quitarContorno()
]),
logo: false,
layers: layers,
target: 'map',
view: new ol.View({
  center:ol.proj.fromLonLat([-61.530766,-34.865575]),
  zoom: 7
  })
 });
}

This is the function that draws the polygon and then puts an end to it:

function addInteraction() {
 var value = "polygon";
 noDraw();
 draw = new ol.interaction.Draw({
      source: source,
      type: "Polygon",
    });

    map.addInteraction(draw);

    draw.on('drawend', function(evt){
      var feature = evt.feature;
      var p = feature.getGeometry();
      var format = new ol.format.WKT();

      poligonoWKT  = format.writeGeometry(p.transform(ol.proj.get('EPSG:3857'),ol.proj.get('EPSG:4326')));
      p.transform(ol.proj.get('EPSG:4326'),ol.proj.get('EPSG:3857'));
      console.log(poligonoWKT);

      noDraw();

      $("#listgrouparea").show();
      getArea();

    });

}


texts: Generate delaunay triangulation from 3D(x,y,z) points

texts:
I have 3D Points array.each 3D point have x,y,z data. I want to generate delaunay triangulation.I am new in delaunay triangulation. My array as below.
-256.24 85.35 5.0
-262.24 84.35 15.0
-263.24 87.35 25.0
................
I am working in c++ with Visual studio IDE.
I want to generate indexes of triangulation. As well as array of triangles points.


texts: SaaS for high performance geographic intersections?

texts:
I have a fairly simple set of requirements, but I'm struggling to find anything that will handle them "out of the box" at scale. This could be my ignorance with GIS terms, or that there is nothing.

Essentially, we would like to load a system with thousands, or even millions of geographic "areas" (polygons / multi-polygons). A user will define an area/ line/point. We then need to see

a) what stored areas does the user geometry intersect with and

b) get the geometry defined by that intersection. The application is very much read heavy, with the stored "areas" changing infrequently and not in real time.

The SQL Server geographic data type seems to have this functionality (as does PostgreSQL) and would be a good solution for storing the data, but I worry about scale in terms of processing intersections with masses of data and many concurrent users (public facing website).

Cosmos and Mongo DBs seem to allow intersection detection, but not geometry generation.

I was expecting there would be some sort of a "server" (with many GPUs probably) that I could deploy across AWS and handle trillions of these sorts of operations per second, but so far my searches have come up empty. This is surprising because what I'm doing seems like a common use case.


texts: Converting coordinate system of feature layers and datasets programmatically

texts:
I am trying to convert coordinate system of featurelayer and datasets using ESRI.ArcGIS.DataManagementTools.Project . The tool works fine in ArcCatalog .

I have doubts about the parameters i am giving. here is the code :

private void button5_Click(object sender, EventArgs e)
    {
        try
        {

            IWorkspaceFactory pFactory = new FileGDBWorkspaceFactory();

            string pPath = folderBrowserDialog1.SelectedPath;

            string pFileName = System.IO.Path.GetFileName(pPath);

            IWorkspace workspace = pFactory.OpenFromFile(pPath, 0);

            ESRI.ArcGIS.Geodatabase.IFeatureWorkspace featureWorkspace = (ESRI.ArcGIS.Geodatabase.IFeatureWorkspace)workspace;

            IFeatureLayer layer = new FeatureLayerClass();
            layer.FeatureClass = featureWorkspace.OpenFeatureClass(comboBox2.Text);

            ILayer layer2 = (ILayer)layer;
            IFeatureDataset fdataset = featureWorkspace.OpenFeatureDataset("test");

            ESRI.ArcGIS.Geoprocessor.Geoprocessor gp = new ESRI.ArcGIS.Geoprocessor.Geoprocessor();
            // Set the parameters
            ESRI.ArcGIS.DataManagementTools.Project proj = new ESRI.ArcGIS.DataManagementTools.Project();

            proj.in_dataset = layer;
            proj.out_dataset = OutputFiletxt.Text;
            proj.out_coor_system = layer2.AreaOfInterest.SpatialReference;

            RunTool(gp, proj, null);
        }
        catch (Exception ex) { throw ex; }

    }

private void RunTool (Geoprocessor geoprocessor, IGPProcess process, ITrackCancel TC)
            {
            // Set the overwrite output option to true
            geoprocessor.OverwriteOutput = true;

            try
            {
            geoprocessor.Execute (process, null);
            ReturnMessages (geoprocessor); // (error Message = Error HRESULT E_FAIL has been returned from a call to a COM component. )

            }
            catch (Exception err)
            {
            Console.WriteLine (err.Message);
            ReturnMessages (geoprocessor);
            }
            }

The error is pretty vague and I dont know if there is an other way to trace this problem.
COMPlus exception code (-532462766)


texts: How to handle a geotiff map together with points with OpenLayers (or other libraries)?

texts:
I create a lot of printed maps using various applications (ArcGIS, QGIS, Microstation, Illustrator etc). Sometimes I'd like to turn this map somewhat interactive on the web.

What I'm looking for is a simple way to turn a huge geotiff into tiles and present them in OpenLayers (or similar) together with some points or polygons on top. If it wasn't for the points and polygons I guess I could have sticked to OpenLayers using zoomify tiles.

I want to stay in the local projection and because of that I guess it has to be OpenLayers since Modest Maps etc only handles the spherical Mercator projection?

Is there a tool for slicing up my tiff into tiles and then presented so that OpenLayers can handle it together with other point features in the same projection?

Do I need a TMS/WMS? A tileserver?


texts: Accessing the internal data structure storing Postgis spatial index (PostGres GiST)

texts:
I am curious about how the internal data structures and storage mechanism of the spatial GiST index in Postgis. In Oracle, I used to like to show how the SDO indexes were "just another table", and how the levels were attributes of the indexes, and the bboxes of the indexed geometries were stored as blobs (that could be extracted). HOw does postgis store the indexes, and where?

The system tables http://www.postgresql.org/docs/8.3/static/catalogs.html can be used to identify indexes, but I do not know how to access the actual storage ( if they are indeed tables)


texts: Problem with WMS/WMTS Layer Image Quality in QGIS

texts:
I have been using an aerial provided by the state. I add the layer as a WMS/WMTS. Typically the quality is on-par with the resolution depicted in the online hosted mapping program (GeoWeb). However, after use the image suddenly turns pixelated. There is another person in my office using the same version of QGIS (most recent update) and the aerial as a WMS/WMTS with no resolution issues. Once the pixelation begins it does not seem to go away, even if I remove and reload the layer, and the resolution issues are carried through when printed or made a PDF. What could be causing the resolution issues and what could be done to restore it back to the original quality?

Please let me know if additional information is needed, I am just beginning to work with QGIS and may have left information out due to inexperience.

A screenshot of my resolution issues and a link to the original aerial are attached. Please note the "zig-zag" appearance of the building outlines. They usually start as more crisp, particularly at this scale.

http://geodata.state.nj.us/imagerywms/Natural2015?


texts: Determining when layer has been loaded in ArcObjects?

texts:
Is there an event in ArcObjects to determine whether a layer (specifically a shapefile) has been loaded into a map?

I'm using ArcObjects 10.4.1 and C#


texts: How find Nearest points using STBuffer and STContains in SQLServer using given point from point table?

texts:
This is the Query I am using now..

 SELECT top 1 coloumname FROM tablename WHERE ROUND
    ((geometry::Point(74.8427116215336,12.852857058953,4326).STDistance(the_geom))*111.325,2) <2
    ORDER BY ROUND((geometry::Point(74.8427116215336,12.852857058953,4326).STDistance(the_geom))*111.325,2)

How can I change this to using STBuffer?


texts: How do you split a circle polygon into equal quadrants in QGIS?

texts: How do you split a circle polygon into equal sectors that are 5 degrees in width in QGIS? I have shapefile from a boundary of a field and want to break it into sectors.


texts:
I'm trying to import some data from a SHP to GFT (Google Fusion Table).
I get a 401 error. I'm using an Ubuntu 14.04 LTS Linux 64 bits with GDAL 1.11.0, released 2014/04/16. GFT is available (tested with ogrinfo --formats |grep GFT)

The steps at the moment are :

Go to https://console.developers.google.com/project
Create a project
Activate Google Fusion as one of the available API
Create credentials with OAUTH for Installed applications with Other to get a client ID and a client secret
Export variables in system with

export GOA2_CLIENT_ID=client ID>

export GOA2_CLIENT_SECRET=<client secret>

Both override GDAL_CLIENT_ID and GDAL_CLIENT_SECRET default parameters
Run python /usr/bin/gdal_auth.py to open a browser to accept connexion where you get an authorization token (from the browser after connexion). In the console, you will get this message below

Please login and authorize access in web browser...
Enter authorization token:

and will need to paste the token. It will generate a refresh token with message in the console like below:

Consider setting a configuration option like:
GFT_REFRESH_TOKEN=<refresh token>

export with export GFT_REFRESH_TOKEN=<refresh token>
Try ogr2ogr -f GFT "GFT:access=<authorization token>" ~/countries.shp
Get an error like followed

ERROR 1: HTTP error code : 401
ERROR 1: HTML error page returned by server
GFT driver failed to create GFT:access=<authorization token>

I also tried with GFT_REFRESH_TOKEN=<refresh token> ogr2ogr -f GFT "GFT:access=<authorization token>" ~/countries.shp but go the same error

According to some researches, 401 errors in GFT are related to authorisation or scope but I should miss the point here.

So, any clue is welcome!


texts:
I have this code to plot a choropleth from a GeoPandas DataFrame.

def create_gpd_dataframe(self
    d = {'Angle of point': yard_list, 'Latitude': yard_lat_list,
    'Longitude': yard_lon_list}
    df = pd.DataFrame(d)
    self.geopd = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))
    self.optimal_gdf = (self.geopd[(self.geopd['Angle of point'] >= 88) & (self.geopd['Angle of point'] <= 100)])
    self.optimal_gdf.crs = {'init': 'epsg:4326'}
    a = self.optimal_gdf['Latitude'].tolist()
    b = self.optimal_gdf['Longitude'].tolist()

    self.polygon_optimal_list = list(zip(a, b))

    lat = (self.geopd[(self.geopd['Angle of point'] >= 49) & (self.geopd['Angle of point'] <= 88)])
    lon = (self.geopd[(self.geopd['Angle of point'] >= 100) & (self.geopd['Angle of point'] <= 150)])

    self.polygon_plot = list(zip(lat, lon))

def plotter(self,
            midlat=midlat,
            midlon=midlon,
            gdt1=gdt1,
            gdt2=gdt2,
            target=target
            ):

    f_map = folium.Map(
        location=[midlat, midlon],
        tiles='Stamen Terrain',
        zoom_start=10
    )

    folium.Marker(
        location=[gdt1[0], gdt1[1]],
        popup='GDT1'.format(gdt1),
        icon=folium.Icon(color='red')
    ).add_to(f_map)

    folium.Marker(
        location=[gdt2[0], gdt2[1]],
        popup='GDT2'.format(gdt2),
        icon=folium.Icon(color='red')
    ).add_to(f_map)

    folium.Marker(
        location=[target[0], target[1]],
        popup='Target Area'.format(target),
        icon=folium.Icon(color='red')
    ).add_to(f_map)

    folium.Choropleth(
        self.optimal_gdf['geometry'],
        key_on='feature.id',
        fill_color='YlGnBu',
        fill_opacity=1,
        line_opacity=0.2,
        bins=[0, 20, 30, 40, 50, 60]
    ).add_to(f_map)

    f_map.save('in5.html')

    new = 2
    url = "/home/yovel/PycharmProjects/TriangulationCalc/in5.html"
    webbrowser.open(url, new=new)

get = PlotOnMap()
get.create_gpd_dataframe()
get.plotter()

And for some reason my output is this:


texts:
ST_Envelope and ST_OrientedEnvelope have visual examples in the docs:

How would a meaningful example for ST_Extent look like?


texts:
Got this code from a teacher and I don't know what I am doing wrong! I have searched around but I can't find the error. Visual studio shows no error.

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">

<title>OpenLayers map</title>


texts:
I am trying to Reclassify raster using

pygeoprocessing.geoprocessing.reclassify_raster(
      base_raster_path_band, value_map,
      target_raster_path,gdal.GDT_Int16,
      target_nodata,values_required=True)

in this my input raster is DEM with min and max value of -1024 to 6608 in pixel which i want to reclassify in 3 class (-1024 to 0),(0 to 3000) and (3000 to 6608)
and to replace each pixel value with 1,2 and 3
how do i make value_map dictionary for this.


texts: QGIS labels change their position as map is dragged

texts:
I spent a while fixing all the labels across my map. Unfortunately, they tend to change their position, as the map is dragged.
I need a bigger map than displayed in the QGIS map canvas.
See the problem below:

As the map has been dragged south, the labels concerned move too. It applies to all labels across the map.
Is there any solution for it?


texts: "GDALGetActualBlockSize" error in QGIS while loading any GDAL/OGR related functions

texts:
I'm getting an error in QGIS 2.18.10 in openSuse Leap 42.2 which I have could not solve so far. When launching the program the following error message arises: "Unable to load GdalTools plugin. The required "osgeo [python-gdal]" module is missing. Install it and try again". After that message I get another one about the "Processing" plugin as follows:

Couldn't load plugin 'processing'

ImportError: /usr/lib64/python2.7/site-packages/osgeo/_gdal.so: undefined symbol: GDALGetActualBlockSize
Traceback (most recent call last):
  File "/usr/share/qgis/python/qgis/utils.py", line 308, in loadPlugin
    __import__(packageName)
  File "/usr/share/qgis/python/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/__init__.py", line 28, in
    from processing.tools.dataobjects import *          # NOQA
  File "/usr/share/qgis/python/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/tools/dataobjects.py", line 37, in
    from processing.algs.gdal.GdalUtils import GdalUtils
  File "/usr/share/qgis/python/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/algs/gdal/GdalUtils.py", line 32, in
    from osgeo import gdal
  File "/usr/share/qgis/python/qgis/utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/lib64/python2.7/site-packages/osgeo/__init__.py", line 21, in
    _gdal = swig_import_helper()
  File "/usr/lib64/python2.7/site-packages/osgeo/__init__.py", line 17, in swig_import_helper
    _mod = imp.load_module('_gdal', fp, pathname, description)
ImportError: /usr/lib64/python2.7/site-packages/osgeo/_gdal.so: undefined symbol: GDALGetActualBlockSize

Python version: 2.7.13 (default, Mar 22 2017, 12:31:17) [GCC]
QGIS version: 2.18.10 Las Palmas, exported

Python Path:
/usr/share/qgis/python
/home/XXXX/.qgis2/python
/home/XXXX/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python27.zip
/usr/lib64/python2.7
/usr/lib64/python2.7/plat-linux2
/usr/lib64/python2.7/lib-tk
/usr/lib64/python2.7/lib-old
/usr/lib64/python2.7/lib-dynload
/usr/lib64/python2.7/site-packages
/usr/lib64/python2.7/site-packages/PIL
/usr/local/lib64/python2.7/site-packages
/usr/local/lib/python2.7/site-packages
/usr/lib64/python2.7/site-packages/gtk-2.0
/usr/lib/python2.7/site-packages
/usr/lib64/python2.7/site-packages/wx-3.0-gtk2
/home/XXXX/.qgis2//python

Same happens with any GDAL-related plugin like  Semi-Automatic Classification:

The python-gdal package has been installed in the system since long ago. I'm not fully sure but I think the problem started after installing an additional instance of GDAL 2.1.3 from the sources (make, make install...) because it was the only accepted version by Geonode (installed in a python virtual environment). What makes me suspicious about this option is that the "About" section shows that QGIS is running against version 2.1.3 of GDAL/OGR (see picture below) but the only approach I've had with that version was the previously mentioned Geonode installation and currently there is no trace of it in my system (or at least I can't detect it).

I've tried everything I know but nothing works so far, like uninstalling and cleaning every trace of GDAL libraries (both 2.2.1 and 2.1.3 versions), reinstalling python, qgis, grass and all the GIS-software, reinstalling GDAL and the python bindings either from Yast or pip (both are version 2.2.1), changing vendor for QGIS and all GDAL dependencies either to Science or Geo repositories, performing several distro upgrades, deleting QGIS configuration folders and files, etc., but the problem remains.


texts: Concatenating XY values in attribute table without messing up numbers in ArcGIS for Desktop?

texts:
Here's the procedure I use and the resulting problem I'm having in ArcGIS 10.1.

First I use the Add XY coordinates tool to add POINT_X and POINT_Y fields to the attribute table of a point layer. I then add a new text field so that I can concatenate the XY values using the field calculator and the following query: [POINT_X] & "," & [POINT_Y]. For some reason, most of the entries in the concatenated field go to many more decimal places than I want.

For example: POINT_X = -97.831771, POINT_Y = 29.782709, and the concatenated field (XY) = -97.8317708001,29.7827091999. I just want the XY field to equal -97.831771,29.782709! What is causing this error? Somehow I've been able to work around this problem occasionally but I can never reproduce it when I want to.


texts: Is there image editing software that will not alter a GeoTIFF's embedded georeferencing?

texts:
I know there are a couple workarounds for this but I was hoping a direct solution exists.

I am looking for an image editor (like Photoshop but hopefully free software) which can be used to edit GeoTIFFs without removing their georeferencing information upon saving/exporting the edited file. I know that I can create a .TFW file to avoid this problem, but I'd like to avoid creating an extra file and would prefer to keep them as GeoTIFFs (.TIFF images with the georeferencing information embedded in them).

I've tried the GIMP, PhotoFiltre, LazPaint, and Paint.NET, none of them save the TIFF's georeferencing data when saving the edited image.


texts: Is there an equivalent of ArcGIS Hyperlink script editor in QGIS?

texts:
I am trying to find similar functionality in QGIS to ArcGIS Hyperlink script editor.

In ArcGIS I can parse VBScript from a selected field in an attribute table. The script would execute on mouse click.

ArcGIS 10.1 Hyperlinks help

Can I achieve the same task in QGIS using VBScript?


texts: Determining if collection of coordinates (polygon) is ellipse?

texts:
I'm attempting to determine a geometry type based on a collection of coordinates and have come across a situation where I'd like to differentiate between what is a "polygon" and what is an "ellipse". This question is specific to ESRI's ArcGIS Runtime for WPF 10.1.1 SDK, but I imagine this is generic enough to have bearing in any GIS. I'll use ESRI's API for the examples in this sample.

Given the following code sample to generate an ellipse's point collection:

double slice = 2 * Math.PI / 360;
double radiusX = 50;
double radiusY = 20;
        ESRI.ArcGIS.Client.GeometryMapPoint center = new ESRI.ArcGIS.Client.Geometry.MapPoint(0,0);
        ESRI.ArcGIS.Client.Geometry.PointCollection pointCollection = new ESRI.ArcGIS.Client.Geometry.PointCollection();
        for (int angle = 0; angle <= 360; angle += 6)
        {
            double rad = slice * angle;
            double px = center.X + radiusX * Math.Cos(rad);
            double py = center.Y + radiusY * Math.Sin(rad);
            pointCollection.Add(new ESRI.ArcGIS.Client.Geometry.MapPoint(px, py));
        }

And then, given this sample to generate a polygon's point collection (obviously the polygon could be much more complex than this):

ESRI.ArcGIS.Client.Geometry.PointCollection pointCollection = new ESRI.ArcGIS.Client.Geometry.PointCollection();
        pointCollection.Add(new MapPoint(0, 0));
        pointCollection.Add(new MapPoint(0, 10));
        pointCollection.Add(new MapPoint(10, 10));
        pointCollection.Add(new MapPoint(10, 0));
//Close the polygon
pointCollection.Add(pointCollection[0]);

Is there an effective, efficient, and generic way to determine which of these two shapes is an ellipse and which is not? The impetus behind this is that the ESRI WPF API does not differentiate between Polygons and Ellipses.

Perhaps I can make this question more clear, what I would like to discern is whether the given points constitute what could be considered a 2-Dimensional ellipse (perhaps already making this too subjective). The ellipse could have any number of radial points comprising it, and I would like to ideally determine if it meets some test of "roundness" or curvature. The x-radius and y-radius of the sample ellipse could also be varied. I've edited my sample to include this.


texts: Can ESRI JS API show PolyLine as Bezier curve?

texts:
Can ESRI JS API somehow show PolyLine as Bezier Curve?

If yes, how to do it?


texts:
I have a problem with the import of data with FME 2012 and ArcSDE 10.
After import, my data are not visibles in ArcMap or ArcCatalog.
My layers are registered in ArcSDE.

I found a solution:
After import, i write this command line:
sdelayer -o delete -l .....
and
sdelayer -o register

My feature are visible in ArcMap.
Have you a solution to resolve this problem ?


texts: ExecuteError: ERROR 000539 while running expression rcexec()

texts:
from arcpy import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'D:\work\USA'

#pathway to all rasters in workspace directory
rasters = arcpy.ListRasters('*.tif')

# Local variables:
Output_raster = "test_script.tif"

# Local variables:
aa = "flam_length_new.tif"
bb = "burn_prob.tif"
cc = "us_drought_prob.tif"
dd = "100m_fire_prec.tif"

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("(\"flam_length_new.tif")*(Lookup(\"burn_prob.tif\",'Probabilit')*0.6)+(Lookup(\"us_drought_prob.tif\",'Ptobabilit')*0.2)+(Lookup(\"us_fire_1000m_proba.tif\",'prob')*0.3)", Output_raster)

The error is:

  Warning (from warnings module):   File "", line 1
  SyntaxWarning: import * only allowed at module level

  Traceback (most recent call last):   File
  "D:\work\scripts\risk_score.py", line 30, in
      arcpy.gp.RasterCalculator_sa("(\"flam_length_new.tif\")*(Lookup(\"burn_prob.tif\",'Probabilit')*0.6)+(Lookup(\"us_drought_prob.tif\",'Ptobabilit')*0.2)+(Lookup(\"us_fire_1000m_proba.tif\",'prob')*0.3)",
  Output_raster)    File "C:\Program Files
  (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\geoprocessing_base.py", line
  510, in
       return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 000539: Error running expression: rcexec()  Traceback (most
  recent call last):   File "", line 1, in    File
  "", line 8, in rcexec RuntimeError: ERROR 999998: Unexpected
  Error.

  Failed to execute (RasterCalculator).


texts: Check the minimum pixel-value in a Raster using PostGIS

texts:
I'd like to check some basic statistic of a raster in my database. How can I determine the smallest raster value found in all raster tiles (enroaching all srid's)? My code looks like:

SELECT min(ST_Value(rast)) as min, max(ST_Value(rast)) as max FROM <rastername>

Throws out: function st_value(raster) does not exist.


texts:
where is the Python library for QGIS 2.8 on a mac (10.10)?
I cannot find it...
sys.path in the QGIS python console gave me some path, but where the interesting part is happening it replaces the folder with dots.

'/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing' and '/Applications/QGIS.app/Contents/MacOS/../Resources/python

I need to set my python path to use the QGis.core library which requires the Path.
I got the instruction from this site http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html but it seems not be up to date anymore...


texts: How to increase a PostGIS query efficiency when querying OSM information

texts:
I have this SQL query to get roads from my PostGIS server.

WITH boxed_roads AS
(SELECT * FROM "planet_osm_line"
WHERE ST_Within(way, ST_MakeEnvelope(xmax ,xmin, ymax, ymin, 4326))
)

SELECT jsonb_build_object('type', 'FeatureCollection',
'features', json_agg(ST_AsGeoJSON(json_res)::jsonb)) FROM
(SELECT * FROM boxed_roads
WHERE (highway = 'motorway' or highway = 'trunk' or
highway = 'secondary' or highway = 'primary' or
highway = 'unclassified' or highway = 'residential' or highway = 'tertiary'
or highway = 'track' or highway = 'service')
ORDER BY ST_Length(way) Desc LIMIT 200) as json_res;

I'm trying to increase the efficiency of this postGIS query as it is taking on average 350ms and that seems too much.
Anything I should change query-wise?

I've been reading the PostGIS Documentation and found talk about indexes but haven't found a way to manually configure to use them. I guess it happens within geospatial functions like ST_Within.

Also, this is OSM data and each row is filled with None value properties.

Did anyone have the same issue? I found questions here about unions and joins but I only want to query. The information in this DB is static and will probably never change until I update it at some point. No handlers will make any changes to it.

This is my first GIS related question so if you need me to add more information, let me know.
This DB is PostgreSQL 12 with PostGIS 3+
I'm running this query with psycopg2 on python 3.8.1.

Here's the EXPLAIN ANALYZE output:

"Aggregate  (cost=54725.32..54725.33 rows=1 width=32) (actual time=359.728..359.730 rows=1 loops=1)"
"  ->  Subquery Scan on json_res  (cost=54686.90..54699.56 rows=100 width=1131) (actual time=236.965..252.877 rows=100 loops=1)"
"        ->  Limit  (cost=54686.90..54698.56 rows=100 width=1115) (actual time=236.934..250.602 rows=100 loops=1)"
"              ->  Gather Merge  (cost=54686.90..60576.64 rows=50480 width=1115) (actual time=236.926..250.370 rows=100 loops=1)"
"                    Workers Planned: 2"
"                    Workers Launched: 2"
"                    ->  Sort  (cost=53686.87..53749.97 rows=25240 width=1115) (actual time=226.442..226.520 rows=57 loops=3)"
"                          Sort Key: (st_length(planet_osm_line.way)) DESC"
"                          Sort Method: top-N heapsort  Memory: 197kB"
"                          Worker 0:  Sort Method: top-N heapsort  Memory: 272kB"
"                          Worker 1:  Sort Method: top-N heapsort  Memory: 235kB"
"                          ->  Parallel Bitmap Heap Scan on planet_osm_line  (cost=252.45..52722.22 rows=25240 width=1115) (actual time=20.177..192.332 rows=30106 loops=3)"
"                                Filter: ((highway = ANY ('{motorway,trunk,secondary,primary,unclassified,residential,tertiary,track,service}'::text[])) AND st_within(way, '0103000020E61000000100000005000000E14830254D515DC013CFFFB8E3465DC0E14830254D515DC0EE5807D4D26C4040FCD56DF3B5674040EE5807D4D26C4040FCD56DF3B567404013CFFFB8E3465DC0E14830254D515DC013CFFFB8E3465DC0'::geometry))"
"                                Rows Removed by Filter: 6318"
"                                Heap Blocks: exact=1117"
"                                ->  Bitmap Index Scan on planet_osm_line_index  (cost=0.00..237.31 rows=3319 width=0) (actual time=28.226..28.226 rows=109272 loops=1)"
"                                      Index Cond: (way @ '0103000020E61000000100000005000000E14830254D515DC013CFFFB8E3465DC0E14830254D515DC0EE5807D4D26C4040FCD56DF3B5674040EE5807D4D26C4040FCD56DF3B567404013CFFFB8E3465DC0E14830254D515DC013CFFFB8E3465DC0'::geometry)"
"Planning Time: 1.707 ms"
"Execution Time: 360.202 ms"


texts:
This is my first time with QGIS and maybe it is an easy question but I can't figure out how to cut part of a layer to add it to another one.

In the image you can find (red circles) parts of a layer totally embedded in other layers. I would like to "cut" that parts and join/add to the surrounding layer.

I have tried to freehand selection tool but after defining the selected area the whole layer is selected. I just want to select the embedded piece, then create a new layer and add to the surrounding one.

Shapefiles have only one attribute, so I can't (do not want) to select by attribute but for region.

Shapefiles available in this file (Dropbox link)

I'm using QGIS 2.18.3 under Ubuntu 16.04.

As a novice any idea or advice is welcome


texts: Comparing two point shapefiles that contain lat/longs and compute distance using ArcGIS Desktop?

texts:
I have 2 sets of point shapefiles in ArcGIS Desktop (group A and B), each set has around 1000 entries that contain lat/long.

What I want to do is compare Group A against Group B and figure out for every point in Group A what point in Group B is the closest and what that distance is.

What would be the easiest way to do this in ArcGIS Desktop?


texts: Why is Yellowstone so empty on OpenStreetMap?

texts:
I stumbled upon the map of Yellowstone National Park on OpenStreetMap:

Is there a particular reason for this map to be so strange?

the borders of the park cut off the "forest" marking, right through the borders
the inside of the national park is almost empty, the roads are not marked, the village is missing, etc.

What is particularly enigmatic to me is that this is a very touristic area and I would expect it to be very well updated.

I checked for instance Zion National Park or the Grand Canyon and they are rich in information. I also had a look at Yellowstone on Google Maps and there are plenty of details.


texts: How to align/georeference distorted vector layers in QGIS?

texts:
I've got a shapefile which has been digitized from a strongly distorted historical map, and would like to align it as nicely as possible to another precise shapefile layer.

The transformation must not be affine : since the map is historical, it has both cartographic and measurements errors.

The best I've achieved for now is to use :

ogr2ogr -gcp 01_a_x 01_a_y 01_b_x 01_b_y -gcp 02_a_x 02_a_y 02_b_x 02_b_y ... -tps

The results is not bad at all (see images below), but the deformation is not "precise" enough.
I mean that the points I defined do not align perfectly in the resulting geometry.

I would like to have those point align perfectly. Of course, points that are not set using -gcp would not have the same precision since they are only interpolated.

Is there a way to achieve that ?

Thanks a lot !!

Olivier

Step 1 : misaligned layer in red, reference layer in gray

Step 2 : points that I defined using -gcp

Step 3 : resulting layer. Not bad, but as you see, the points I defined are not perfectly aligned...


texts: How to highlight kml features in openlayers

texts:
I'm trying to make the features from a kml layer highlight on hover, but it doesn't seem to work for me. On the other hand, if i create a vector layer with a polygon feature it works just fine. I made an example with kml and polygon feature just to try to figure out what the difference is. Does anyone know what I am missing?
This is the .js code:

var map;

function init(){
        map = new OpenLayers.Map('map');

        var kml = new OpenLayers.Layer.Vector ('Australia Boundaries',{isBaseLayer:true, protocol: new OpenLayers.Protocol.HTTP ({url: 'kml/australia.kml', format: new OpenLayers.Format.KML({ extractStyles: true, extractAttributes: true})}),projection: new OpenLayers.Projection("EPSG:4326"), strategies: [new OpenLayers.Strategy.Fixed()]});
        map.addLayers([kml]);

    //create vector layer with a polygon feature

        var vector = new OpenLayers.Layer.Vector('Vector', {projection: new OpenLayers.Projection("EPSG:4326") });

        var point = new OpenLayers.Geometry.Point(90.752725,10.475584);
        var pointList = [];
        for(var p=0; p<6; ++p) {
            var a = p * (2 * Math.PI) / 7;
            var r = Math.random(1) + 20;
            var newPoint = new OpenLayers.Geometry.Point(point.x + (r * Math.cos(a)),
                                                         point.y + (r * Math.sin(a)));
            pointList.push(newPoint);
        }
        pointList.push(pointList[0]);

        var linearRing = new OpenLayers.Geometry.LinearRing(pointList);
        var polygonFeature = new OpenLayers.Feature.Vector(
            new OpenLayers.Geometry.Polygon([linearRing]));
        vector.addFeatures([polygonFeature]);
        map.addLayer(vector);

   //create highlight and select controls for the vector and kml layers

        var highlightCtrlVector = new OpenLayers.Control.SelectFeature(vector, {
            hover: true,
            highlightOnly: true,
            renderIntent: 'temporary'});

        var selectCtrlVector = new OpenLayers.Control.SelectFeature(vector,
            {clickout: true}
        );
        var highlightCtrlKml = new OpenLayers.Control.SelectFeature(kml, {
            hover: true,
            highlightOnly: true,
            renderIntent:'temporary'});

        var selectCtrlKml = new OpenLayers.Control.SelectFeature(kml,
            {clickout: true}
        );
        map.addControl(highlightCtrlVector);
        map.addControl(selectCtrlVector);
        map.addControl(highlightCtrlKml);
        map.addControl(selectCtrlKml);
        highlightCtrlVector.activate();
        highlightCtrlKml.activate();
        selectCtrlVector.activate();
        selectCtrlKml.activate();

        map.addControl(new OpenLayers.Control.LayerSwitcher());
        map.zoomToMaxExtent();
    }

And here is the example: http://geoteste.comli.com/app/indexTest.html


texts: How does OpenLayers 3 handle coordinates?

texts:
Ol3 always uses the flatCoordinates because obviously it needs 2d-coordinates as the screen is only 2-dimensional. But how does OL3 calculate the flat representation?

For example if my view’s projection was EPSG:25832 (UTM Zone 32) it could calculate the transversal Mercator Projection (with help from proj4.js) to obtain the flat coordinates. How is this processed with, for example, EPSG:4326 (WGS84), which is using geographical coordinates? Which projection does OL3 apply here?


texts: Load a layer from a PostgreSQL using Graphical Modeler in QGIS after filtering with multiple input variables

texts:
I have a table in a PostgreSQL database. I want to pull in features from it as a layer that match the Date range and Lat/Long specified as an input while running a QGIS Graphical Modeler. Some other functions would then be performed on this layer within the model.
I tried to use the 'PostgreSQL execute and load SQL' algorithm within the model, but I can only input one variable at the end of a query by using the 'Pre-calculate Value' option.
How can this be undertaken?


texts: Presentation of terrain and sub-terrain data in 3D using Qgis

texts:
I have DEM of karst and .shp files of caves in that area exported from "speleoGIS" Therion
http://therion.speleo.sk/

Qgis2threejs is amazing in making visualisation of 3D sourface, which is easy to share and present.
But how can I add there models of caves, when Qgis2threejs allows me to use only polyline file?
And it is possible to add there such format as .dxf?


texts: I have large point set which I want to overlay on open-Street map. That is to create vector tiles of the point set based on zoom levels of the underlying Open-Street map. For example if the current zoom level of open-Street map that I'm viewing is 5, I want to create vector tile of the point set for the same zoom level, and then overlay this. With this approach, I'm trying to cluster the points and serve as tiles.Is there some one who have used and succeeded implementing  this approach? or does this approach make sense?