texts: Can I change the color of selected query in QGIS2.18.9 (Las Palmas)?

texts:
It is a shapefile of Gujarat roadways, and I tried to select some features which appear selected in the attribute table but not in the map. Shapefile is perfect.
Selected features are related to specific SH from Ahmdabad to Modasa.


texts:
I'm loading a WMS layer which comes in WGS84 into Google Maps. From what I've read so far, I should reproject this layer to Google Pseudo Mercator (either 3857/900913). However, when I save the layer in either one of those two projections, it is not displayed in Google Maps, but I can view it in Data Preview of Geoserver. Others with WGS84 are displayed correctly but with a position shift.

Do I really need to reproject the layer? I'm not sure whether it is not just a matter of some inaccuracy.


texts:
I'm trying to convert coordinates to DMS, I try to use field calculator, but I still have no result.

I also try this script:

(CASE WHEN $x < 0 THEN '-' ELSE '' END) ||
floor (abs($x)) || '°' ||
floor(((abs($x)) - floor (abs($x))) * 60) ||'\'' ||
substr(
  (tostring(
    (((abs($x)) - floor (abs($x))) * 60) - floor(((abs($x)) - floor (abs($x))) * 60)) * 60),
  1,
  5) || '" E'"

Any proposition?


texts: Cannot get specific WMS to display in Leaflet or OpenLayers 3

texts:
I am trying to create a map in both Leaflet and OpenLayers 3 which displays the  layer EMAP8 from the WMS http://wms.nlsc.gov.tw/wms. I can't get it to work. I can open the layer in QGIS, but when I try to use it in Leaflet like this:

var overlay_taiwanemapen0 = L.tileLayer.wms('http://wms.nlsc.gov.tw/wms', {
    layers: 'EMAP8',
    format: 'image/png',
    transparent: true,
    continuousWorld : true,
});

or in OpenLayers 3 like this:

var lyr_taiwanemapen = new ol.layer.Tile({
    source: new ol.source.TileWMS(({
        url: "http://wms.nlsc.gov.tw/wms",
        params: {"LAYERS": "EMAP8", "TILED": "true"},
    })),
    title: "taiwane-map(en)",
});

I get nothing displayed. I can see layers being requested, and data being transferred, but nothing on the map.

What am I doing wrong?


texts: How to calculate flight line intersections?

texts: I am new to GIS, and I have flight lines and I am interested in determining the crossover or intersection of those lines. Can I use QGIS for that and what would be the procedures for doing it?


texts: Bilinear resampling yields inaccurate/small values (Earth Engine)

texts:
I am trying to upscale / resample some GSMaP observations from 0.1 x 0.1 degrees to 30m. I want the resolution to match Landsat 8 so that I can perform some modeling.

But in Google Earth Engine when I resample using bilinear interpolation, the GSMaP values don't geographically overlap the original values well, and they are also significantly smaller. Because the original cells are so large, I thought at least some of the new resampled cells would have the same values as the original image, given that values are based on the 4 nearest neighbors. Instead, the values are significantly smaller. E.g. the red GSMap cells have a precip rate of 42 mm/hr, the resampled cells have a rate of 4 mm/hr.

Any ideas? I'm not sure if I am resampling incorrectly in Earth Engine, or misunderstanding the technique.

https://code.earthengine.google.com/73e2becea4a5ff86594f8bb40a46183a

PS. The resampled layer is transparent in Earth Engine but changing the opacity does not make it opaque. Would like to know if there is a way to fix that as well.


texts:
In LASTools you can reproject LiDAR files using LAS2LAS.exe by selecting the the target CRS under target projection. But I can't seem to find an option in the UI to convert from NAD83 to WGS_1984_UTM_Zone_11. I can only find UTM and a drop down menu for Zones.
How can I do this transformation in LASTools?


texts: Proper ground-to-grid correction and projection from old survey maps

texts:
I asked this question over on the ESRI geonet site, but haven't had any luck there,  sorry for the duplication...

I'm working on digitizing a group of maps from the 1950's that show locations of easements along a pipeline. The easement boundaries are recorded with quadrant bearings and distances in feet. (I'm not sure if it's international feet or survey feet, and according to a few sources it shouldn't make much difference? Would love to have someone confirm that though.)

I'm using the Traverse tool to draw the easement parcels, and trying to figure out if I need to do a ground-to-grid correction. The original maps indicate that the basis of bearings is "State Plane Coordinates (Lambert System) California Grid Zone V," with no specific indication of declination or any other corrections. I want the file that I'm creating to be in State Plane CA Grid Zone V (Lambert) as well. So, my first question is: Do I even need to bother with a ground-to-grid correction? And if so, can anyone shed any light on how to find the right correction angle?

Along with this, I'm wondering about keeping everything accurate since I'll need to go from NAD27 (on the original maps) to NAD83 (to play nice with current files). My plan of attack right now is to:
1) create the easement feature class and set the coordinate system to NAD27 in the appropriate SP zone
2) Draw easement polygons
3) Project into NAD83.
Would this workflow keep the features accurate? (Or, any suggestions on other/better ways to do it?)

Hope this is clear; dealing with datum/projection details isn't my forte.


texts: Trying to use Python to clean up null/bad data in Lat/Lon Load

texts:
I'm using ArcMap 10.4.  Trying to load a list of Lat/Lons and connect the surface hole and bottom holes with a line connector (Toolbox XY to Line).  I have two issues:

1) Values which have BHL and no SHL.  I took care of those with a definition query Surf_Lat >0 in the data table which just doesn't plot the wells which don't have a surface hole.

2) Values which are from vertical wells and have a BHL Lon = 0 and BHL Lat = 0.  What I need to do is make the BHL = SHL if the values are zero.  I've not been successful in trying to do a definition query because I can't figure out how to do the statement from (1) above, and also do an IF statement.

I need to fix this somehow as you can see in the background I have lines going crazy on the map trying to "connect" to this 0 Latitude value. I've taken some Python classes but not sure about assigning the name of the value here in my equation.

Code

def mycalc(Surf_Lon, Surf_Lat):
    if (Surf_Lon=0):
      return mycalc=Surf_Lat
    else:
     return mycalc=Surf_Lon


texts: What is the category limit for category widgets in CARTO BUILDER?

texts:
Is there a way to set the numbers of items displayed in the category widget? Is 5 the current limit?

I can not find anything in the documentation:
https://carto.com/learn/guides/widgets/search-and-filter-category-values

I can use "Search in categories" but then the items are not selectable anymore.


texts: Is there any way to extract coordinates of a road segment using QGIS?

texts:
I add my own road vector data on top of Google Maps using the OpenLayers plugin with QGIS. The vector layer is generated using a PostgreSQL DB.

Is there any way to extract coordinates of a road segment?


texts:
I'm trying to import some CAD files from microstation (dgn) v7 to ArcGis 9.3.1 and i need to access the User Data Linkage from it but i have not the Data Interoperability License.

Is there any alternative way to do this?


texts:
I am new to GIS. I would like to share a ArcMap project that I initiated with someone else.

Is there a way to share directly this project with the data, to allow the person to keep working on it?


texts: split polyline with point and add vertex on that point as well as reflect in attribute table?

texts:
I am new to ArcGIS.

I am using the ArcObjects SDK and I couldn't understand how to:

Add vertex to that clicked location as well
Add attributes in that attribute table.

Can someone suggest some books so that I can get better knowledge of ArcGIS.


texts: Selection of only monthly data on Google Earth Engine?

texts:
Can you suggest an approach to select January Landsat 5 data between 2000 to 2010 and calculate its mean?
In other words, the mean of the all January images from 2000 to 2010.


texts: How to delete char in table attribute in QGIS 2.6.1

texts:
I need to label many items but in the table attribute field there are some char (#F) for every item that I have to delete.

The length of the description are not always the same.

I create a new field but I am not able to set the logic operation to do that.


texts: IDW: Default Maximum Distance and Coincident Points

texts:
The ArcGIS Help file says the Default Maximum Distance, when using the Variable Search Radius in IDW, is the extent's diagonal, however, no value is automatically filled in the corresponding box after I define the extent, so is ArcMap really using this default value if the box is empty?

Regarding coincident points (different measurements at the same location) when running IDW, it is possible to set the Environment "Geostatistical Analysis" to Min, Max, Mean, etc., but is this selection really honored when running IDW using Spatial Analyst, or does it only apply to tools from the Geostatistical Analyst extension? I am confused because the ArcGIS help says: "Tools that honor the Coincident Points environment define how coincident data is treated in Geostatistical Analyst".
If the Coicident Point selection is not honored when using Spatial Analyst, how do I know what IDW is doing about them?
Thanks a lot for your help!


texts: Round gdal_translate AAIGrid output to 4 decimal places

texts: How can I use gdal_translate for AAIGrid so that the raster values are not more than 4 decimal places.


texts:
I'm using  Postgres 9.4.4 with Postgis 2.2.1., my issues is that i have about 3 million of points i want to fetch and display on a map with leaflet.
The problem is that those are too many points, what is the best aproach to tackle the problem since fetching all the rows from the database is slow?

i have considered taking a bunch of random points and group the most closest one to those points but it might not be best.

Other solution  i thougth is paginate all the points and fetching them by groups, since this would be faster that query all of them in one time, but not sure how to do this with Geoserver.

Any hints/ideas about this?


texts: Publishing a vector layer from GeoServer to use in Mapbox

texts:
I have published a vector file in GeoServer 2.15.1. How can I add this layer to Mapbox?
The example in doc.geoserver.com works fine but when I use my layer in the provided format it keeps showing 404 error for the request.
After that, I used the code below but still can't get any progress.

mapboxgl.accessToken = 'pk.eyJ1IjoibW9oYW1tYWRob3NzZWluOCIsImEiOiJjazZrYmM5YmEwMmQ2M2RwNmppMzBvbHd3In0.Ki5iFg47J1PJB3W0Z6wlgQ';
    var map = new mapboxgl.Map({
        container: 'map',
        style: 'mapbox://styles/mapbox/streets-v11',
        center: [80, -50],
        zoom: 4
    });

    map.on('load', function() {
        var layers = map.getStyle().layers;
        // Find the index of the first symbol layer in the map style
        var firstSymbolId;
        for (var i = 0; i < layers.length; i++) {
            if (layers[i].type === 'symbol') {
                firstSymbolId = layers[i].id;
                break;
            }
        }
        map.addSource(
        'large_area_tight_modelling', {
            'type': 'vector',
            "tiles": [
                "http://localhost:8080/geoserver/gwc/service/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&LAYER=test:large_area_tight_modelling&STYLE=&TILEMATRIX=EPSG:900913:{z}&TILEMATRIXSET=EPSG:900913&FORMAT=application/vnd.mapbox-vector-tile&TILECOL={x}&TILEROW={y}"

            ],
            "minZoom": 0,
            "maxZoom": 14
        });

       map.addLayer(
        {
            'id': 'large_area_tight_modelling',
            'type': 'fill',
            'source': 'test:large_area_tight_modelling',
            'layout': {},
            'paint': {
                'fill-color': '#f08',
                'fill-opacity': 0.4
            }

        },
        firstSymbolId
    );
});


texts: Replicating"Intersection" function in QGIS using PostGIS queries

texts:
I have two multipolygon shapefiles that I've imported into my PostGIS database - one of congressional districts (cds) and one of US counties (ctys). Below shows the counties in blue and the districts overlain with orange border, no fill.

INPUT DATA

I would like to perform a simple intersection using a PostGIS/PostgreSQL query that would produce the same output as the "Intersection" tool in QGIS, which creates a multipolygon output that consists of polygons for all intersected areas between the two files, as well as all of the attributes of each file. This desired output is shown here:

DESIRED OUTPUT (QGIS/Intersection)

I tried to achieve this in PostGIS using ST_Intersection() along with ST_Intersect() like this:

CREATE TABLE test AS
SELECT ST_Intersection(a.geom, b.geom) AS geom
FROM cds AS a, ctys AS b
WHERE ST_Intersects(a.geom, b.geom);

but the resulting output is a messy, multiple geometry-type output with polygons, lines, points and gaps in seemingly random places.

BAD OUTPUT (PostGIS/ST_Intersection())

The geometries of each input layer (checked with the QGIS/Check Validity function) appear to be valid.

What am I doing wrong?


texts: Creating a TopoJSON file with countries that also includes sovereign countries

texts:
I am looking for a TopoJSON file that includes all countries including what I ended up after looking as 'sovereign countries' found in the following link:

https://en.wikipedia.org/wiki/List_of_sovereign_states

Given the fact that the TopoJSON will be used in D3.js, it is stated that the TopoJSON files used in D3 are using vectors from 'Natural Earth' found here:
http://www.naturalearthdata.com/downloads/110m-cultural-vectors/

However, the list of sovereign countries found in the 'Natural Earth' data is not the same - eg while Monaco is listed in wikipedia as a sovereign country it is not listed in the sovereign countries of natural earth. It is listed in the 'populated places' file though.

I am looking for a file that includes all countries including the sovereign countries - wikipedia version. I don't mind creating what is needed from scratch, so here I am asking for efficient ways to achieve what I want.

Googling around I found out that you can merge multiple GeoJSON files and then convert to TopoJSON, however the processes aren't very clear.

Looking for any pointers, tips, tools - getting a ready file would be an added bonus, but I would like to learn how to do it on myself.

In other words, is it easier to find a list, strip it and merge it?

Is there a way to manually add vectors and additional information on a TopoJSON directly?

Do I need to extract information from it, edit and re-build it somehow?


texts: Why curl returns "userName%" by GeoServer REST API create new user

texts:
When I insert this below link like GeoServer 2.17.1 RESTful API I got only: userName%
curl -u admin:@dmin9901 -X POST http://localhost:8080/geoserver/rest/security/usergroup/users/ -H  "accept: application/json" -H  "content-type: application/json" -d "{  \"userName\": \"test\",  \"password\": \"test1234\",  \"enabled\": true}"

It can't create new user in my GeoServer, it would be bug?


texts: How can I render custom google maps tiles with mapnik?

texts:
I'm currently developing my own little GIS service, that displays custom computed information on a google map.

As you can see, I'm loading a custom tile overlay, where I renderd the information in a color scale. The problem is with the Google Maps texts for places like Berlin. Those are now partially covered by my tiles and hard to read. I would like to have the text displayed above my tiles.

Can I render the Google Maps texts for places and streets using mapnik? Or can I put my custom tiles between the text and the map image (but I dont think those two information are splitted into separated images, so it might not be possible)?


texts: Delete duplicate geometries from ORACLE table

texts:
I have searched a while now, but there seems to be no answer including an example on how to find and delete duplicate geometries from an Oracle DB table.
I am facing a situation where I have a table with 3million+ entries with duplicate point geometries (exact same location) and I need to remove them. The other attributes are not helpful, they are allowed to have duplicates. So it doesnt't matters which record is deleted.
So in my example I would expect to end up with the entries with pnt_id 1 and 3.

CREATE TABLE spatial_table (
    pnt_id number,
    pnt SDO_GEOMETRY);

INSERT INTO spatial_table VALUES (
    1,
   SDO_GEOMETRY(
      2001,
      NULL,
      SDO_POINT_TYPE(12, 14, NULL),
      NULL,
      NULL));
)

INSERT INTO spatial_table VALUES (
    2,
   SDO_GEOMETRY(
      2001,
      NULL,
      SDO_POINT_TYPE(12, 14, NULL),
      NULL,
      NULL));
)

INSERT INTO spatial_table VALUES (
    3,
   SDO_GEOMETRY(
      2001,
      NULL,
      SDO_POINT_TYPE(20, 18, NULL),
      NULL,
      NULL));
)


texts:
I am exporting tables from Postgresql 9.2 to shapefiles using postgis 2.0. I have the same system on other machines and it works correctly. I have tried running the command from a prompt and get an error about the table not existing, with a truncated version of the name.

"path_to_exe/pgsql2shp" -f "..\shapefiles\W06_Exporters_150814_141145_with_geom" -h localhost -u postgres -P
    postgres -p 5432 db_name "select * from analysis.\"W06_Exporters_150814_141145_with_geom\""


texts: From the thousands of the overlaid polygons, get clusters of density

texts:
I have a few thousands of the polygons, over the given region. Many of them overlay. What I want is to detect areas, where density (number of overlaying polygons at the point) is substantially higher than around.

My first attempt is to create a grid of points and for each point count the number of polygons that it intersects. From that, I can get an overall heat map and spot the areas of the highest density.

However, as density varies a lot, I imagine some minor (local) spots remain hidden, as total numbers are not that hight. So I am looking for a technique that will allow me to detect spatially segregated clusters that have the density that is way higher than their surrounding.

It sounds like a contour detection in CV, so I really hope to find the code for the problem


texts:
I have a table that contains a column polygon shape. I need a query that I Want to pass a point(coordinates) if point contains in a polygon_shape of an instance the match will be returns .

select * from "DistrictAreas" where ST_Contains(polygon_shape, ST_GeomFromText('POINT(51.3525284 35.6960617)', 4326));

this is a error comes from postgresql

  Reason: SQL Error [XX000]: ERROR: Operation on mixed SRID geometries

  I find that this is the problem :
  org.jkiss.dbeaver.model.exec.DBCException: SQL Error: Error parsing
  SQL data:java.lang.IllegalArgumentException: Points of LinearRing do
  not form a closed linestring


texts: Can I make a multi layers pdf from a mapnik xml?

texts:
I am playing with mapnik and trying to export multiple vector files as one pdf with multiple layers.

For the moment, my vectors get merged when rendering. Is there a way to keep each vector in a separate layer of my rendered pdf?


texts: Remove islands and completely surrounded polygons after polygonization with QGIS

texts:
I have a line map of buildings in a city. I was able to correctly polygonize the buildings in QGIS with the Polygonize processing algorithm. All the holes/islands became features as expected. Now I need to select and remove all features that are completely surrounded by other (one or more) features. Examples shown in red on the image:

These are not part of buildings, but inner yards.


texts: Draw and save graphics to IndexedDB / Local Storage using ArcGIS API for JavaScript?

texts:
I need to add graphics to a map with class.

Is it any sample or way to save it to IndexedDB or Local Storage ?

I tried to use offline-editor-js but this library don´t offer possibility to save graphics.


texts: Exporting table from ArcGIS Desktop to Excel file?

texts:
Is there a way to export a table with selected records to an Excel spreadsheet?

I was trying to export just to .txt file but ArcGIS gives me an error message.

  An error occurred exporting the table

The other problem I have is with coordinates. I want them only with one locus after the comma, I mean like the format: xx xx xx,x not like I have now xx xx xx,xxx.

The version of my software is 9.3.


texts:
I create geodatabase "training"

sudo su
su postgres
createdb -E UTF8 template_postgis2
createlang -d template_postgis2 plpgsql
psql -d postgres -c "UPDATE pg_database SET datistemplate='true' WHERE datname='template_postgis2'"
psql -d template_postgis2 -f /usr/share/postgresql/9.1/contrib/postgis-2.1/postgis.sql
psql -d template_postgis2 -f /usr/share/postgresql/9.1/contrib/postgis-2.1/rtpostgis.sql
psql -d template_postgis2 -c "GRANT ALL ON geometry_columns TO PUBLIC;"
psql -d template_postgis2 -c "GRANT ALL ON geography_columns TO PUBLIC;"
psql -d template_postgis2 -c "GRANT ALL ON spatial_ref_sys TO PUBLIC;"
createdb training -T template_postgis2

geodatabase training http://savepic.su/3233935.jpg

Then I changed "default.style"

sudo cp "./default.style" "/usr/share/osm2pgsql/default.style"

Then I load shapes

sudo apt-get install subversion
mkdir "/home/oleg/mapnik"
svn co "http://svn.openstreetmap.org/applications/rendering/mapnik" "/home/oleg/mapnik"
bash /home/oleg/mapnik/get-coastlines.sh

Now I want create style

python /home/oleg/mapnik/generate_xml.py osm.xml my_osm.xml --dbname training --user postgres --password postgres --accept-none

link http://savepic.su/3233934.jpg

Please, tell me how create style.

Update: Now I have this problem.

oleg@oleg-pkvb:~/mapnik$ ./generate_xml.py osm.xml osm_new.xml --dbname training --user osm --world_boundaries /home/oleg/world_boundaries/ --accept-none
Traceback (most recent call last):
  File "./generate_xml.py", line 204, in <module>
serialize(template_xml,options)
  File "./generate_xml.py", line 80, in serialize
mapnik.load_map(m,xml,True)
RuntimeError: Postgis Plugin: FATAL:  Peer authentication failed for user "osm"

Connection string: ' dbname=training user=osm connect_timeout=4'
  encountered during parsing of layer 'landcover' in Layer at line 341 of 'osm.xml'

Update 2:

oleg@oleg-pkvb:~/mapnik$ ./generate_xml.py osm.xml osm_new.xml --dbname training --host localhost --port 5432 --user osm --password postgres --world_boundaries /home/oleg/world_boundaries/ --accept-none
Traceback (most recent call last):
  File "./generate_xml.py", line 204, in <module>
    serialize(template_xml,options)
  File "./generate_xml.py", line 80, in serialize
    mapnik.load_map(m,xml,True)
RuntimeError: ERROR:  relation "planet_osm_polygon" does not exist

Update 3

I try download data in my database "training"

oleg@oleg-pkvb:~$ sudo osm2pgsql -U osm -d training -W postgres -H localhost -P 5432 Moscow.osm.gz
osm2pgsql SVN version 0.70.5

Password:
Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point" does not exist, skipping
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
PREPARE get_way (int4) AS SELECT AsText(way) FROM planet_osm_point WHERE osm_id = $1;
 failed: ERROR:  function astext(geometry) does not exist
LINE 1: PREPARE get_way (int4) AS SELECT AsText(way) FROM planet_osm...
                                     ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Error occurred, cleaning up

Update 4:

Now I use osm2pgsql 0.83, postgis 2.1, postgres 9.1, mapnik 300000. I create user 'gisuser' and database 'gis'. But I have the same problem ((

postgres=# ALTER USER gisuser WITH PASSWORD 'postgres';
ALTER ROLE
postgres=# \q
oleg@oleg-pkvb:~$ sudo osm2pgsql -U gisuser -W postgres -d gis -H localhost -P 5432 Moscow.osm.gz
osm2pgsql SVN version 0.83.0 (64bit id space)

Password:
Warning: program compiled against libxml 209 using older 207
Using projection SRS 900913 (Spherical Mercator)
Setting up table: planet_osm_point
NOTICE:  table "planet_osm_point" does not exist, skipping
NOTICE:  table "planet_osm_point_tmp" does not exist, skipping
SELECT AddGeometryColumn('planet_osm_point', 'way', 900913, 'POINT', 2 );
 failed: ERROR:  AddGeometryColumn() - invalid SRID
CONTEXT:  SQL statement "SELECT AddGeometryColumn('','',$1,$2,$3,$4,$5, $6)"
PL/pgSQL function "addgeometrycolumn" line 5 at SQL statement


texts:
I want to remove existing vector layers on the map and

map.getLayers().getSource();

returns me ol.source.Source object as it stated in the API.

I cannot find the way how to jump from ol.source.Sourceto ol.source.Vector object. Everything I find in answers was to use something like this

map.getLayers().getSource().getFeatures()

which assumes that getSource() returns ol.layer.Vector I believe (generates error).

How to get to ol.layer.Vector from map?


texts:
I found data about offshore Windparks in Germany I'd like to import into QGIS.
According to this tutorial it should be quite easy:

http://www.igismap.com/convert-kml-shapefile-qgis/

But all I get are empty layers. The attribute table is empty too.

The dataset is downloaded from here:

http://www.geoportal.de/DE/Geoportal/Karten/karten.html?lang=de&wmcid=111

I download the Windfarm data under ''Speichern'' and I get a small 5,5kb file. Importing it into QGIS works technicially but the layer is empty.


texts:
What formula could I use to calculate the size in Kilometres of a bounding box based on a given South-West latitude/longitude and a North-East latitude/longitude points?

The bounding box format is defined as:

bounds = sw_latitude, sw_longitude, ne_latitude, ne_longitude


texts: Using GeoJSON Services in Arcade Expression

texts:
I am trying to use an Arcade expression with GeoJSON services in a map using the ArcGIS Javascript API. The intent of the Arcade expression is to return a count of lakes that intersect a county that the user clicks on. Both the lakes and county polygons being GeoJSON services. Swapping out the lakes GeoJSON service with a ArcGIS Feature Layer service results in the expression functioning properly, returning the count of lakes that intersect the county GeoJSON service.
  var template = {
  title: "{LABEL}",
  content: "{expression/countyInfo}",
  expressionInfos: [
  {
    name: "countyInfo",
    title: "Marker Info",
    expression: `

    var lake = FeatureSetByName($map, 'Lakes');

    var lakeIntersect = Intersects(lake, $feature);

    var cnt = Count(lakeIntersect);

    return \`
    Number of lakes: \${Text(cnt, "####")}
    \`
    `
  }
  ]
};

I've tried changing the $map global variable in the FeatureSetByName function to the other global variables options ($datastore, $feature, $layer) but have had no luck.
Has anyone had any success creating a FeatureSet from a GeoJSON service, how did you declare it?
EDIT
Code can be found here


texts:
When I try to add a WFS-layer in QGIS, it crashed automatically. It generates a .dmp file, which I opened in Notepad++, but I don't know where to look for the actual cause of the crash.

Tried various WFS-links, all crash. Tried various options (change CRS, zooming to certain area, "only request features overlapping the view extent", using different titles of the same WFS-server... etc...

Some of the WFS links I've used:

https://services.geodan.nl/public/data/org/gws/MBIZ4280LEEF/MBIZ4280LEEF/wfs?

//geodata.nationaalgeoregister.nl/rwsgeluidskaarten/wfs?request=GetCapabilities

//geodata.nationaalgeoregister.nl/bag/wfs?request=GetCapabilities

Any ideas on what the problem might be?

Windows 10, 64-bit
QGIS Desktop 2.18.8


texts:
I'm relatively new to PyQGIS and I'm building a plugin to perform some calculations between raster layers. I was seeking to use the QgsRasterCalculator for some operations, however, the function that I wrote return an invalid raster layer. I think the problem might be in the output path (I already printed the other parameters to make sure, and everything seems correct). The layer is not created there, and therefore the variable result does not pick it up. Do you know if there is any way to solve this?

def calculator(self, layers, rates100, rates5000, entries):
    formula = str(rates100[0]) + '*' + entries[0].ref + '+' + str(rates100[1]) + '*' + entries[1].ref + '+' + str(rates100[2]) + '*' + entries[2].ref + '+' + str(rates100[3]) + '*' + entries[3].ref + '+' + str(rates100[4]) + '*' + entries[4].ref + '+' + str(rates100[5]) + '*' + entries[5].ref + '+' + str(rates100[6]) + '*' + entries[6].ref + '+' + str(rates100[7]) + '*' + entries[7].ref + '+' + str(rates100[8]) + '*' + entries[8].ref
    calc = QgsRasterCalculator(formula,
                               'D:/Teste/teste.tif',
                               'GTiff',
                               layers[5].extent(),
                               layers[5].width(),
                               layers[5].height(),
                               entries)
    calc.processCalculation()
    result = QgsRasterLayer('D:/Teste/teste.tif','Vulnerabilidade')
    return result


texts: Getting the distance in km between two coordinates within Hungary

texts:
I have a table called properties (generated with EF Core) which has a column called GeoPoint. In the database, this column is a geomerty(Point) type. All properties are located within Hungary. I am trying to create a procedure for search properties within a distance of a specific coordinate.
Distance should be in meters.
Properties table has two text columns, longitude, latitude, which I have used to populate this new extra column called GeoPoint and I have used the following script to do this.
UPDATE "Properties" SET "GeoPoint" = (ST_MakePoint(CAST("Longitude" AS double precision), CAST("Latitude" AS double precision)))

After this I had all my GeoPoint columns populated, but not sure if this was the correct way.
To search I tried to use the following script
SELECT p."Id", p."Longitude", p."Latitude", p."GeoPoint"
FROM "Properties" p
WHERE ST_DWithin(
        p."GeoPoint",
        ST_Transform(
            ST_GeomFromText(
                'POINT(19.081750000000 47.54770000000)',
                4326),
            2100),
        5000)

and I get the following error;

ERROR:  LWGEOM_dwithin: Operation on mixed SRID geometries (Point, 0) != (Point, 2100)

I am completely new to PostGIS, but it seems like the script I have used to populate the GeoPoint columns is incorrect.
UPDATE:
I have changed my update script to
UPDATE "Properties" SET "GeoPoint" = ST_SetSRID((ST_MakePoint(CAST("Longitude" AS double precision), CAST("Latitude" AS double precision))), 23700)

and my select statement to
SELECT p."Id", p."Longitude", p."Latitude", p."GeoPoint"
FROM "Properties" p
WHERE ST_DWithin(
        p."GeoPoint",
        ST_Transform(
            ST_GeomFromText(
                'POINT(19.081750000000 47.54770000000)',
                4326),
            23700),
        500000)

Running this select statement returns no error, but it does not return any results either. Distance is set to 500000, which should be 500 km
and the table contains multiple sample record within 3 to 20 km distance.
If I change the distance to 5 000 000, then it returns all rows even some located 17+ km away
UPDATE AGAIN:
I have googled it and SRID 23700 supposed to be the correct one Hungary and distance should be in meters. I still do not get it. If I run the following:
select ST_Distance(
'SRID=23700;POINT(19.067931856332 47.511673637045)'::geometry,
'SRID=23700;POINT(19.081750000000 47.54770000000)'::geometry)

It returns 0.03858548849237025 which is clearly not in meters. What am I doing wrong?


texts:
I've got some issues with OpenLayers 3 when displaying a Vector layer with opacity. The more I pan or zoom in the map, and the less transparency there is on the features. I think this comes from the Vector Source that loads several times the same features, which makes some kind of overlap on the same data. Here is my code for WFS layer, where myStyle is a ol.style.Style with opacity=0.3 :

var vectorSource = new ol.source.Vector({
    loader: function(extent) {
        var url = gsWFS + '?service=WFS&version=1.1.0&request=GetFeature&typename=' + lyrName + '&outputFormat=application/json&srsname=EPSG:3857&bbox=' + extent.join(',') + ',EPSG:3857';
        $.ajax({
            url: url,
            dataType: 'json',
            jsonp: false,
            method: 'POST',
            success: function (response) {
                vectorSource.addFeatures(geojsonFormat.readFeatures(response));
            }
        });
    },
    strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({ maxZoom: 20 }))
});

var vecLyr = new ol.layer.Vector({
    source: vectorSource,
    style: myStyle,
    visible: true
});

Here are some iamges of the map, on rendering (opacity=0.3) and after a few moves in the map (opacity should still be the same, it is in the object, but it's not on the map).

Do you have any ideas of how to deal with this opacity problem ?


texts: How to rename field names in a shapefile attribute table without table manager using python

texts: By using the table manager we can rename the field names but i want to rename of fields using python. How can i do that?


texts: Free 3rd party library for showing geospatial vector information

texts:
Up to now, I have used GDAL to extract raster data from geotiffs and then, OpenGL to show that raster data.
But when it comes to extract vector data with the help of OGR, what library would you suggest to show these vector data in visual-studio picture controls that is more appropriate than OpenGL?

I've found out from your comments that you have thought I'm searching for a way to read vector data. This is the subject of my another question. Here, I'm just asking suppose that we have extracted all the geometries. Which environment is the best to show these geometries?


texts: Average frequency of zip code change for a given address over the course of year

texts:
This is a difficult question to answer with a high degree of accuracy, but could someone estimate to within an order of magnitude for me the chances of a zip code for a given address changing over the course of a year?

For example, would you say the chances of this are:

a) ~10%

b) ~1%

c) ~.1%

d) ~.01%

Just trying to get a rough sense of the constancy of these things for some code I'm writing. If the answer were ~10% I would proceed very differently than if it were ~.01%.


texts: Clipping large IUCN range map vector file to shapefile using ArcGIS Desktop?

texts:
I'm having trouble clipping my large IUCN species range maps vector file to a southeast asia boundary shapefile. The single vector file consists of many individual species polygons. I've tried eliminating most of the polygons that are non SEA, made sure both shapefiles are in the same projection, and then clip it to SEA boundary shapefile, but i get an error message saying I've generated an empty output. After reading some answers online, I tried repairing geometry but the same issue arise, empty output. I've tried QGIS but that doesn't work either,software crashes.

What could be the problem and how can I solve this?


texts: Changing Area Output (somewhere to 8th decimal point) each time ArcPy script run?

texts:
I currently have shapefile of US PUMAS and another feature shapefile (call it feature X). I am using ArcPy to project both of these files to the following projection coordinate system (using Project_management):

http://spatialreference.org/ref/sr-org/nad_1983_contiguous_usa_albers/html/

After projecting, I use intersect_analysis and then CalculateAreas_stat to compute the intersection areas. Finally I output a .csv file, which I then use to compute area weighted means.

The issue is that every time I run my simple script, the intersection areas vary by a small amount (somewhere to the 8th decimal point). These small changes are impacting my analysis.

If this is a GIS precision issue or if I making a mistake, how can I go about fixing it?


texts: Passing a file into ModelBuilder based on user input string

texts:
What I need to do is pass a string into my model to tell it to look in a certain file for a State's road network kind of like this:

C:\Users\User1\Documents\ArcGIS\USA\%State%\Shape\GU_IncorporatedPlace.shp

Where "%State%" is a user input string. I've tried the Parse Path tool but I'm not sure how to load the file into the Model.


texts: Geocoder of Google Maps on OpenLayers map

texts:
I'm trying to integrate the geocoder of Google Maps into my ol-map. There are a lot of similar questions asked. I tried out this: http://dev.openlayers.org/sandbox/topp/almanac/examples/geocoder.html example but it doesn't work. Looking through: How to use Google Maps geocoding with Openlayers? also didn't help me.

(As I'm new here I cannot comment on peoples answers so I thought it might be OK to open a new question - tell me if I'm wrong) See code below for my attempt. Do you have any tips or know any working example?

<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Google styled in OL mit kml</title>
<link rel="stylesheet" href="OpenLayers-2.12/theme/default/style.css" type="text/css">
<link rel="stylesheet" href="OpenLayers-2.12/theme/default/google.css" type="text/css">

OSM + Google Maps + KML Reprojection

 <input id="address" type="textbox" value="Sydney, NSW">
 <input type="button" value="Geocode" onclick="search()">


texts: Stopping classbreakvalues from changing symbology of point layer (remove strokes) using ArcPy?

texts:
I'm using ArcMap 10.4.1.

I am trying to update a point layer symbology with a premade lyr file, then adjusting it's class break values in arcpy. It's working fine, only it changes the symbols back to default (while the colours stay correct). I don't want my points to have a stoke/outline and can't change that with arcpy, so thats a problem.

Here's my code:

updateLayer = arcpy.mapping.ListLayers(mxd, "FC", df)[0]
FClyrpath = "FCsymology.lyr"
sourceLayer = arcpy.mapping.Layer(FClyrpath)
arcpy.mapping.UpdateLayer(df,updateLayer,sourceLayer, symbology_only = True)

lyr = arcpy.mapping.ListLayers(mxd, "FC", df)[0]

if lyr.symbologyType == "GRADUATED_COLORS":
    lyr.symbology.valueField = "depth"
    lyr.symbology.classBreakValues = [3, 5, 8, 11, 12, 15, 18, 21, 24, 72]

arcpy.RefreshTOC()

This is what my premade ".lyr" symbology looks like:

This is what happens when changing the classbreakValues: The classification is good, the strokes around the symbol not.

I also tried changing the default symbology in AdvancedArcMapSettings, but no avail.


texts: How to disallow edits in layers on a QGIS project?

texts:
So I'm adding layers to QGIS from Postgresql, is there a way to "lock" certain layers and make sure that users don't accidentally move stuff around or delete them?

Ideally something done on the project level would be ideal but if it can't be done then how to prevent it from the db?


texts: How do I perform spatial join of raster data in QGIS? I have raster soil data (tiff), I'd like to join this to an attribute table in excel? I'm not able to view the attribute table either? Is this possible in QGIS, specifically working with raster data?


texts: Is there a way to rasterize polygon with hatch pattern with GDAL CLI or in a Python script?

texts:
I am trying to automate some work previously done in QGIS. Target platform is AWS Lambda so there is no way to use QGIS or any other standalone GIS-specific application.
QGIS provide a way to draw polygon on GeoTIFF, fill polygon with hatch pattern, and then export result as a PNG picture. I am trying to re-create this behavior with Python and GDAL only.
Polygon is stored as feature with single polygon in GeoJSON format. Currently I am able to render polygon with gdal_rasterize and then combine result with original GeoTIFF with gdal_merge.
What's the problem?
I couldn't find a way to programmatically fill (originally stored in GeoJSON format, rasterized, saved in GeoTIFF format) polygon with hatch pattern.
What I tried?
Not much.
First of all, there is a Feature Style Specification page that defines the way feature style information (i.e. colors, line width, symbols, etc.) should be handled at the various levels in GDAL’s vector drivers (OGR). And there is no GeoJSON among vector drivers that provides support of feature style information.
On the other hand, on the same page I found few other vector drivers, that claimed to support feature style information and support georeferencing at the same time: LIBKML, DGN, PDF, JML, MITAB, CAD. None of them is familiar to me.


texts: Field calculatur in ArcMap 10: How to specify field row?

texts:
I'm new to arcmap, so I'm not entirely sure about terminology which is probably why I can't seem to find the answer to this seemingly trivial question.

In the attribute table I want to be able to create a new field which takes the value from another field and adds the value from the previous row to it. My idea was simply to make a some kind of a for loop but to do that it seems I need to be able to specify a row number of a field, and I don't know how to do that.

Any help would be appreciated, ranging from a better solution to the answer to this specific question :).


texts: Using CARTO in a different language (multilanguage)?

texts: Is it possible to configure CARTO, CARTO Builder, etc to be available in different languages e.g. Spanish, French, Portuguese?


texts: OpenLayers - Click on a panel triggers a map click but it should not

texts:
I added a panel to a map that is contained inside an IFrame. Then I attach a click handler to the div of this panel with JQuery to toggle a control.

However, I notice that a click on the panel also trigger a map click event when it should not because another control handle the map click and both the panel and the map controls are triggered together.

var olmap = iframe.contentWindow.map // openlayers map inside the iframe
var i = olmap.getControlsByClass("Korem.Control.Info")[0] //map click control
var layerControl = olmap.getControlsByClass("Korem.Control.LightLayerControl")[0]
var customPanel = new iframe.contentWindow.OpenLayers.Control.Panel({
  id: "customPanel",
  autoActivate: true
});
olmap.addControl(customPanel)
iframeBody = $(iframe).contents().find("body");
var panel= iframeBody.find("#customPanel")
panel.css({"position":"absolute", "top":"0", "right":"0", "width":"auto", "height":"auto", "z-index":"1033"})
panel.append('
Afficher les écoles
');
$(iframe).contents().find("#clickLink").click(function(){
    i.deactivate()
    if (layerControl.active == true) {
      layerControl.deactivate();
    } else {
      layerControl.activate();
    }
});
layerControl.events.on({deactivate : function(){
    i.activate();
  }
});

Is it normal that a click on a panel trigger a map click also? How could I prevent it?


texts: Width and distance when using line patterns in QGIS

texts: I am trying to make a cross hatch polygon fill (example link attached) using the values (defined by law) of 0.1mm for line width and 0.8mm for distance between the lines. But I don't get/see any hatch with such a small values, so I'd like to know if the values less than 1mm are supported anyway or I have to solve this differently?


texts: ID unknown projection for shapefile in ArcMap?

texts:
I am using ArcMap 10.3 and I have a data file with xy coordinates (New York City/State) that I'm trying to project properly.  Having trouble figuring it out and don't have a way to reach whoever created it.  I know there are a million posts on this, but most rely on testing out a bunch of different projected coord. systems after narrowing them down based on the extent.  Here's the extent for this data:

TOP: 5058236.102283 | BOTTOM: 4940052.540569 | LEFT: -8264805.988673 | RIGHT: -8180005.401691

From the above, I guess it'd be a UTM projection (?), but I'm not sure what to do about those negative numbers.

I also have partial address info and geocoded one to try to troubleshoot:

XY geocoded: 40.7218342,-73.9787655

UTM based on that lat/long (easting/northing): 586248.1, 4508380.2

What I have in the data for that point (x/y): -8235276.4, 4971409.736

Nothing I've tried comes even remotely close.

Any ideas?

EDIT: Have an idea what these numbers are (if not a solution).  I gave up on the shapefile above with the wacky XY coordinates and switched to a different shapefile, which had most of the same points but with a projection specified (NAD_1983_UTM_Zone_18N) and added it to a map that was projected using WGS_1984_Web_Mercator_Auxiliary_Sphere.  Later, I was using the "Move to" tool on a point and noted that the pre-populated values in the tool window were in the same wacky coordinates as the other shapefile.  E.g. for the above example, it gave me: -8235272.518, 4971391.237.  My idea is that these are the transformed XY values resulting from the on-the-fly projection, and that the original data somehow had these also.  Maybe it's decimal degrees converted to meters, but relative to the meridian/parallel of the data frame coordinate sys.  That said, I'm not sure what they represent or how to transform them back for use in other maps.


texts: Can an external application start ArcCatalog.exe and immediately navigate to a specific folder?

texts:
I'm looking something similar to explorer's command line switches:
explorer.exe /select,c:\examples\condensed\need\obtains\dispute\under\billion\floor

or

explorer.exe /e,c:\examples\condensed\need\obtains\dispute\under\billion\floor

Can an external application create a folder c:\examples\condensed\need\obtains\dispute\under\billion\floor and then start ArcCatalog.exe in that location?


texts: ArcPy SignInToPortal() Error "Failed to sign out"

texts:
Does anyone know how to sign out of portal using arcpy?
The code below let me sign into portal. Running the script anytime after that causes an error "Failed to sign out". That tells me that it is probably signed into portal but I haven't found a way to sign out. The error traceback is below.
Things I tried:

SignOutFromPortal_server() gives the error. "module 'arcpy' has no
attribute 'SignOutFromPortal_server'". Was this ported to Pro?
Catching the error works but it's a problem if my portal name or
credentials are different than the one that's currently signed in.
Using ProPy.bat gives the same error.
Adding an 'r' before the portal string gives the same error.
Closing everything and logging out of the server changes nothing.
Here is an ESRI Community post with a very similar title but the
answer provided was not helpful.

This is running from in IDLE Python 3.7.9 [MSC v.1922 64 bit (AMD64)] on win32
ArcGIS Server 10.8.1
import arcpy

# This returns "https://gisproweb01p.houtx.lcl/cohportal/"
ptl = arcpy.GetActivePortalURL()

arcpy.SignInToPortal("https://gisproweb01p.houtx.lcl/cohgis", "domain\user", "xxxxxx")

more code...

Traceback (most recent call last):
  File "M:\Workspace\Locators\test_batch\BARC_batch\BARC_Geocode.py", line 25, in <module>
    arcpy.SignInToPortal(r"https://gisproweb01p.houtx.lcl/cohgis", "houtx\arcgismanger", "@arcg15manager")
  File "M:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py", line 2450, in SignInToPortal
    return _SignInToPortal(*args, **kwargs)
ValueError: Failed to sign out.


texts: how to use GPS traffic data existing in MSSQL format, on toQGIS?

texts:
I have GPS traffic data in MSSQL format, consisting of ID, LAT, LONG, SPEED, DIRECTION,TRACK TIME.

How do I use this in QGIS? I tried with MYSQL, but connecting to server was an issue.

As my data is not spatial and does not have any geometry information I cannot use SpatiaLite either.


texts: Import .tif file to server and display it in PuTTy (grass gis)

texts:
I am quite new in working on server and have one problem with importing file.

I have .TIFF files in my computer and would like to send them to server which I work with using PuTTy and GRASS. I tried to use WinSCP to copy files to my folder in server, but I am not sure if it is enough to display it in GRASS via PuTTy (anyways, does not work when I try).

I tried to use: r.in.gdal and r.import. As well just to display file using d.mon and then d.rast filename.

It's basic problem, but please, tell me step by step how to do that to be able to display raster.


texts: Plot an error ellipse at the tip of a vector arrow

texts:
I am using QGIS 3.12 and I would like to plot an error ellipse at the tip of a vector arrow.

I have plotted GNSS velocity vectors in a map and now I would like to plot the error associated to the vector. How could I do it ?

I have tried to add an ellipse in the symbology, but not sure how to locate it at the tip of the vector arrow. The idea is to produce a similar map as the one below. The arrows indicate the displacement direction and velocity (vector length) of Iberia respect a fixed Africa from GNSS measurements and the ellipse the associated error.

The origin of the vector is the location of the GNSS station and the arrow tip location depends on the vector length. Is at that arrow tip that the ellipse would be located.


texts: How to clip an xarray to a smaller extent given the lat/lon coordinates?

texts:
I have a netcdf file showing elevation over Europe derived from GTopo30. Turns out that the extent of this raster file is too big to combine it with other layers, so I would like to crop it to a given extent specified in lat/lon coordinates. I have seen some examples here and here trying to make this operation I am describing, but it seems they are not working in this case.

The raster file description looks like this:

<xarray.DataArray (band: 14, y: 5520, x: 8400)>
[649152000 values with dtype=float64]
Coordinates:
  * band     (band) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14
  * y        (y) float64 71.0 70.99 70.98 70.97 70.96 ... 25.03 25.02 25.01 25.0
  * x        (x) float64 -25.0 -24.99 -24.98 -24.97 ... 44.97 44.98 44.99 45.0
Attributes:
    transform:     (0.0083333333, 0.0, -25.000139509, 0.0, -0.0083333333, 70....
    crs:           +init=epsg:4326
    res:           (0.0083333333, 0.0083333333)
    is_tiled:      0
    nodatavals:    (-1.7e+308, -1.7e+308, -1.7e+308, -1.7e+308, -1.7e+308, -1...
    scales:        (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1....
    offsets:       (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0....
    descriptions:  ('tpi', 'slope', 'aspect', 'Pvec', 'Qvec', 'alt', 'dist2co...

So, as you can see, the longitude and longitudes are stored under the "x" and "y" coordinates and they come in WGS84 (EPSG:4326) format.

The extent is this:

min_lon = -24.995
min_lat = 25.05
max_lon = 45.50
max_lat = 71.55

I tried to use xarray's filters and subsetting, but it seems I am doing something wrong, since I can't get the required slice:

# Example 1
sel1 = band.sel(x=(band.x < max_lon) | (band.x > min_lon))
sel2 = sel1.sel(y=(sel1.y < max_lat) | (sel1.y > min_lat))

# Example 2
sel1 = band.where((band.x < max_lon) | (band.x > min_lon))
sel2 = sel1.where((sel1.x < max_lat) | (sel1.x > min_lat))

Both of them produce: KeyError: "not all values found in index 'x'"

What is the correct xarray way of cropping a large raster to a given lat/lon window?
Thanks!


texts: Landsat 8 bands of 432 looks too dark in ArcGIS

texts:
I have created a composite map using the combination of Landsat 8 bands of 4, 3 and 2 (I read this combination gives true natural color) to see the scene in true color in ArcGIS, but it looks too dark. I couldn't see things clearly from this map.

but I wanted to create this kind of composite map, which is easy to recognise the forest, river etc. How to create this kind of map?

What is wrong with it?


texts: Buffer Problems While Using Model Builder

texts:
I've created a model that is supposed to allow a user to select an area of road and will then highlight the attributes within a user-specified buffer distance from the selected area. This is what I've made:

However, it doesn't show the actual buffer part of it. The final product is supposed to have the area of road and the attributes within the buffer distance show up (which it is doing), but it doesn't show the buffer:

How can I fix this?


texts: Displaying multiple fields as labels if one field is null in QGIS?

texts:
I am having some issues trying to show multiple fields for a single label in QGIS. In QGIS2.6, attempting to display labels for a feature consisting of multiple fields results in an empty label if one of the fields is null. So for example, in the label expression dialog:

concat("Temp Site",'\n',"Notes",'\n',"Function")

will function normally unless one of the fields is null, then it results in no label being shown. Using the '||' operator instead in the expression dialog also results in the same behaviour:

"Temp Site" || '\n'  ||  "Notes" || '\n'  ||  "Function"

Is there a way to display the labels even if one of the fields is null?


texts:
Is it possible to install QGIS Plugins with Python code?
I would like to be able to provide users who have a basic QGIS image served up to them, to be able to quickly install all their required plugins with a tailored Python script to save them manually installing each via the usual menu route, one at a time.
I cannot see any code that runs in logs after a successful plugin install and no clues online in my searches.


texts: Google Maps and Displaying Coordinates while drawing?

texts:
Currently we are doing this on a project:

google.maps.event.addListener(map, 'mousemove', function (event) {
    displayCoordinates(event.latLng);
});

function displayCoordinates(pnt) {
    var coordsLabel = document.getElementById("tdCursor");
    var lat = pnt.lat();
    lat = lat.toFixed(4);
    var lng = pnt.lng();
    lng = lng.toFixed(4);
    coordsLabel.innerHTML = "Latitude: " + lat + "  Longitude: " + lng;
}

However, when we go to draw a polygon this no longer works.  We had to attach another event to the polygon mouseover once it is drawn to call displayCoordinates, but we are having a hard time finding what to call how to implement this with Google Maps.  ESRI and OpenLayers makes this kind of stuff easier IME.


texts: ECW/MrSid plugin for GDAL from Kyngchaos not working on Mountain Lion?

texts:
Subject says it all. I have installed all of Kyngchaos packages, but

gdalinfo --formats

does not list ECW or MrSid formats.
Could someone else with Mountain Lion try if they can get any JP2 driver installed into GDAL?


texts: I have collected data (spectral bands) from Landsat 4-5 and Landsat 8. I was calculating NDVI when I noticed that the two Landsats have different band numbers designated as Red and Near Infra Red (needed for NDVI calculation). When I checked the wavelengths associated with these bands here Landsat Band Description I noticed that the wavelengths also differ a bit. For example I used bands 3 and 4 (Red & Near Infra Red) from Landsat 4-5 and bands 4 and 5 (Red and Near Infra Red) from Landsat 8 for NDVI calculation, but the wavelength range of band 3 in Landsat 4-5 is not exactly the same as of band 4 in Landsat 8. Is there any way to convert things so that they are same and comparable. I got big variances between NDVI values from the two satellites that's why the confusion.


texts: Finding diameter of pipes on either side of valves using ArcGIS Desktop?

texts:
I'm using ArcGIS desktop 10.2 and I need to find the diameter of pipes which are either side of a particular valve. There are over 2000 of these valves on my map so I'll need some sort of tool.

The valves are represented as points/dots and the pipes are lines. The attribute table of the pipe layers have a Diameter field, showing the diameter of each section of pipe.


texts: How to adjust logo with CSS Viewer in QGIS Cloud Pro

texts: I am new to creating web maps and am using QGIS Cloud Pro. I uploaded our company logo onto the web map, but most of the logo is cut off. Using the Viewer tool and CSS there I have tried to adjust it so that all of it shows. I can adjust most other things on the web map using CSS with trial and error, but have not figured out how to get our whole logo to show.


texts: How do I set color of map dots in Carto Javascript API?

texts:
I'm using Carto JS API.

        function main() {
            // create leaflet map and define some properties
            var map=    L.map(
                "map",
                {
                    zoomControl: true,
                    center: [40,-95],
                    zoom: 4
                }
            );

            // add a base layer to map
            var Hydda_Full = L.tileLayer('http://{s}.tile.openstreetmap.se/hydda/full/{z}/{x}/{y}.png', {
                attribution: 'Tiles courtesy of OpenStreetMap Sweden &mdash; Map data &copy; OpenStreetMap'
            });
            Hydda_Full.addTo(map);

            // viz.json of map
            var layerURL=           "URL TO MY CARTO MAP GOES HERE";

            cartodb.createLayer(map, layerURL, {infowindow: false, tooltip: false})
                .addTo(map)
                .done(function(layer){
                    layer.setInteraction(true);
                }
       }

This gives me my map with dots. The default color is orange. I want to change their colors in the JS API. How do I do this?


texts: Merging or joining two tables with some matching columns in QGIS to calclulate unique values

texts:
I have two tables (A and B) that display the presence or absence of a species within a polygon and within each polygons associated buffer layer. The purpose of the buffer layer is to account for a margin of error for species recorded at low spatial accuracy? By that I mean the species presence points are outside the target polygon, but close enough to potentially be within the polygon for example some species were recorded using a 6 figure grid reference which has a 200m margin of error. As you can see within the tables, species may be present in both tables (eg Spp1 in ID 1), however they may also be only present in only one (eg spp6 is present in ID 2 in table B but not A).

I wish to know the number unique species recorded at low spatial within the buffer in addition to the number of unique species recorded at higher spatial accuracy within the polygons but the two must be combined in the same operation to avoid counting the same species as two separate ones. So far I have used a spatial join on the datasets to produce a presence absence matrix (Tables A and B) but this does not work as it duplicates all of the species layers, eg creating spp1A and spp1B rather than just one column for spp1, spp2, spp3 etc. Is there anyway to achieve this?

Table A (eg polygon)

ID    spp1    spp3    spp4    spp6
1     1       1               1
2             1       1
3             1               1
4                     1
5     1

Table B (eg polygon buffer)

ID    spp1    spp2    spp3    spp4    spp5    spp6
1     1       1                               1
2                     1                       1
3                                             1
4                             1
5     1


texts: How to set image as background for a shapefile polygon using QGIS

texts: I would like to change the background of a polygon to an image using QGIS. Under Symbology, there is an option for Raster Fill but the .jpg image is so zoomed in you cannot see it. Is there a way to make the polygon background an image and move it around so it looks good? Any tips on file types, resolutions, or anything?


texts: GeoServer 2.2 and KML with multiple layers

texts:
I would like to be able to create a single URL that allows the download of three layers in a single KML file. Here's two ways to do that:

WMS KML Plain
GeoServer Application Link

Both fail. The first complains of un-matched tags in the KML. The other doesn't load any points from one of the layers. The broken layer in the KML works in OpenLayers fine...

Any ideas on how to fix this, what's going wrong and why?

Upgrading GeoServer may happen later but isn't a fix.


texts: How does the tolerance influence QgsGraphBuilder and tieing of points?

texts:
I am developing a script which uses OpenStreetMap data to calculate the distance (over the road network) between each building and the nearest hospital.

To do this, I mainly follow the official tutorial http://www.qgis.org/en/docs/pyqgis_developer_cookbook/network_analysis.html

My problem comes when calling QgsGraphBuilder. Because my buildings obviously don't lie exactly on the road network, I want to specify a tolerance value when creating the builder.

However, with

builder = QgsGraphBuilder(crs, True, 200.0, "WGS84")

it doesn't "snap" my buildings to the network... Is there anything wrong with the way I call the function QgsGraphBuilder? I haven't seen any example which actually implements this tolerance parameter...

EDIT

I actually also tie all my buildings'centroids to the network:

builder = QgsGraphBuilder(crs, True, 200.0, "WGS84")

all_points = [] # A list that will contain the coordinates of all objects to be tied to the network
all_points.append(hospital) # I add the coordinates of the hospital

for building in coord_buildings.values():
    all_points.append(building) # I add the coordinates of all the buildings'centroids

    tiedPoints = director.makeGraph( builder, all_points ) # I tie the buildings'centroids to the graph

    graph = builder.graph() # I build the graph

    startId = graph.findVertex(hospital) # I take the id of the hospital in the graph

    print startId # Displays -1

If I understand correctly, this means that my hospital doesn't match to anything in the network because it hasn't been tied properly... maybe because of the tolerance value?

By the way, my data is in a PostGIS database, unprojected (EPSG 4326), in case it makes any difference

EDIT 2

I just saw that some buildings are actually tied to the network when I use only one parameter

builder = QgsGraphBuilder(crs)

or

builder = QgsGraphBuilder(crs, True, 0.0)

But rougly 1000 are still not tied. So probably really to do with the tolerance factor?

When I specify the other parameters

builder = QgsGraphBuilder(crs, True, 200.0, "WGS84")

the dijkstra process fails (the tree only contains -1 values)


texts: How do I make a path cross the most possible classes in a raster?

texts: I want to create a function and possibly use Arcmap tools to create something similar to Cost Path. Instead of just choosing the LCP, it must cross through the most different integers. Let's say I have a raster with values 1-5. I want this path to use all five values but still use the shortest path in relation to number of cells. The raster values do not represent any sort of cost.


texts: How to execute GDALRasterize() function wtih options?

texts:
I am using GDAL 2.1.0, in both Windows and Linux environments.
I want to programmatically rasterize a shape file in my C++ app.  Thanks for any clues you can provide.
gdal_rasterize
I can successfully execute the stand-alone gdal_rasterize app from the command line.  For example, the following works to rasterize a shape file that contains areal features about farms.
gdal_rasterize  -a 'crop_id'  -l 'farms'  -init 0  -ts 1024 1024  -ot Byte  farms.shp  farms.tif

GDALRasterize()
The following is my attempt to perform the same rasterization within my C++ app.
GDALAllRegister();

...

//------------------------------
// Paths of the source and destination files, relative to current working directory.
string srcPath = "farms.shp";
string dstPath = "farms.tif";

//------------------------------
// Open the vector dataset.
unsigned int openFlags = GDAL_OF_VECTOR | GDAL_OF_READONLY;
GDALDataset * pSrcDataset = static_cast<GDALDataset *>( GDALOpenEx(
  srcPath.c_str(),
  openFlags,
  NULL, NULL, NULL ) );

// The one layer.
OGRLayer * pSrcLayer = pSrcDataset->GetLayer( 0 );

//------------------------------
// Rasterize options.
// TODO: Are these intended to be just like the command line arguments
// I would use if executing gdal_rasterize?
// https://gdal.org/programs/gdal_rasterize.html
char ** argv = NULL;

// -a <attribute_name>
// An attribute field on the features to be used for a burn-in value.
argv = CSLSetNameValue( argv,  "-a", "crop_id" );

// -l <layer_name>
// The layer from the datasource that will be used for input features.
argv = CSLSetNameValue( argv,  "-l", pSrcLayer->GetName() );

// -init <value>
// Pre-initialize the output image bands with these values.
argv = CSLSetNameValue( argv,  "-init", "0" );

// -ts <width> <height>
// The output file size in pixels and lines.
argv = CSLSetNameValue( argv,  "-ts", "1024 1024" );

// -ot <type>
// Force the output bands to be of the indicated data type.
argv = CSLSetNameValue( argv,  "-ot", "Byte" );

GDALRasterizeOptions * pOptions = GDALRasterizeOptionsNew( argv, NULL );

//------------------------------
// Perform rasterization.
int usageError;
GDALDataset * pDstDataset = static_cast<GDALDataset *>( GDALRasterize(
  dstPath.c_str(), NULL,
  pSrcDataset,
  pOptions,
  &usageError ) );

//------------------------------
// Cleanup.
GDALRasterizeOptionsFree( pOptions );
CSLDestroy( argv );
GDALClose( pSrcDataset );
GDALClose( pDstDataset );

Problem: Invalid arguments
The first problem is that I am providing the rasterization arguments incorrectly.   The following error printed to standard output lets me know that I am doing it wrong.
GDAL code 6: Unknown option name '-a=crop_id'

Providing the arguments incorrectly is probably the cause of my other problems.
Problem: Projection and SRS unknown
The following errors printed to standard output.
GDAL code 10: Pointer 'hDS' is NULL in 'GDALGetProjectionRef'.

GDAL code 1: The input vector layer has a SRS, but the output raster dataset SRS is unknown.
Ensure output raster dataset has the same SRS, otherwise results might be incorrect.

Problem: Crash from uninitialized pointer
I eventually crash because of a bad pointer.
First-chance exception at 0x10050ee0 in make_wcf.exe: 0xC0000005: Access violation reading location 0x00000040.
Unhandled exception at 0x10050ee0 in make_wcf.exe: 0xC0000005: Access violation reading location 0x00000040.


texts: Merge shapes from the same layer by attributes

texts:
I have a single layer with shapes (raster cells) that represent all countries on earth (the figure below shows france). Each shape has specific attributes.

Is it possible to merge the shapes by their attributes? E.g. if attributes x,y and z in shape 1 are the same as attributes x,y and z in shape 2 ... then create a shape that covers the area of shape 1 and shape 2 and contains the attributes x,y,z?

Example: If I have the attribute table displayed below, I would like to merge the shapes with the ID 2 and 7 to one shape.

I have tens of thousands of shapes with hundreds if not thousands of different attributes. I'm afraid doing anything manually is not an option.


texts: 'QVariant' is not defined error in QGIS plugin

texts:
This is the error I am getting when trying to run my plugin.

"NameError: name 'QVariant' is not defined"
This same code will work in the python console in QGIS. So, I don't have an idea as to why it won't work in the plugin that I am making.

from PyQt5.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt5.QtGui import QIcon
from PyQt5.QtWidgets import QAction, QFileDialog
from qgis.core import *
from qgis.core import QgsField

def Create_fields(self, layer):

    list_fields = []

    for field in layer.fields():
        field_id = field.name()
        list_fields.append(field_id)

    # Adding fields to shapefile

    if len(list_fields) == 1:
        res = layer.dataProvider().addAttributes([
            QgsField("Avg_NDVI", QVariant.Double),
            QgsField("Date", QVariant.Date),
            QgsField("Study_Type", QVariant.Int),
            QgsField("Trial_Num", QVariant.Int),
            QgsField("Plot_Num", QVariant.Int),
            QgsField("Trt", QVariant.Int),
            QgsField("Rep", QVariant.Int)])
    layer.commitChanges()

edit

The fix was changing the qVersion to QVariant in the import statement area.


texts: Creating automated maps for attributes in QGIS

texts:
I understand that automated maps can be created using information in the attributes. I want to stop manually creating each map for land owned/managed (this is coded in the attribute tables).

How do I do this?

I am using QGIS 3.8.2


texts: Getting real north direction or real heading?

texts:
I try to get the north direction in a Cesium globe but it seems to always point to the viewer perspective.

Follow the compass instrument on top of the map:

This image I'm pointing (correctly?) to north:

Now, not soo good. I think the north is not there!

Even more strange:

This is how I'm getting the heading:

var viewer = new Cesium.Viewer('cesiumContainer',{
    terrainProvider : terrainProvider,
    timeline: false,
    animation: false,
    baseLayerPicker: false,
    skyAtmosphere: false,
    fullscreenButton : false,
    geocoder : false,
    homeButton : true,
    infoBox : false,
    sceneModePicker : true,
    selectionIndicator : false,
    timeline : false,
    navigationHelpButton : false,
});

var camera = viewer.camera;

var heading = Cesium.Math.toDegrees( camera.heading );

I'm very new on Cesium and took this code from an example.

I need to always point to the real magnetic north.

Seems to get the "false" north or map north I need to:

var heading = 360 - Cesium.Math.toDegrees( camera.heading );

Now I need to point to magnetic north (without to calculate the Magnetic declination )


texts:
Does anyone know of a method to split vector data based on a corresponding raster layer.
For example, I have a classification with two classes and have an overlaid vector layer as seen in the image below. The classification has identified 2 classes in one of the polygons.

What I'm trying to do is split the polygon based on the detected split in the classification.

I've looked into polygonising the classification but when implementing this on a large scale this would talk a huge amount of time which I'm trying to avoid. I've also looked into clustering but couldn't seem to get any good of it.
By any chance has anyone else looked into something similar?


texts: How to edit a layer from a personal geodatabase in QGIS

texts:
I want to open and edit layers from a personal geodatabase in QGIS. I have read this post and I follow the ODBC way to open a personal geodatabase in QGIS with success.

However, I cant edit layers in QGIS.

The editing toolbar is grayed out.

Can I edit the layer from a personal geodatabase in QGIS? Is it possible ?


texts: QField expression 'current_value()' not working

texts:
QField v1.8.4 on a Samsung Galaxy Tab S4 running Android 10.
Project prepared in QGIS 3.14
I have a point layer (trees_pt) - a shapefile. The layer has attributes like inter alia, spreadN, spreadE, spreadS, spreadW.
I would like to set spreadE, spreadS, spreadW to default to the current edited value of spreadN. I have set the default value for each of spreadE, spreadS, spreadW' to current_value('SpreadN'). This seems to have no effect, either in QGIS or QField.
Any ideas what I am doing wrong?


texts:
I have two layers … assume A and B … one is landuse and one is limit of that.

The main reason to union these two layers is to calculate street.

But my problem is after I union these two the area increase.

Logically it should not be more than the limit, because the landuse layer is in limit.


texts: Determining if file in projected or geographic coordinates using ArcGIS Desktop?

texts:
I have download the European grid DEM raster data from Copernicus website, which I would like to use to calculate solar radiation for smaller study site.

The ArcGIS:Area solar radiation tool requires that raster will be in projected coordinate system. So, how can I check if my current coordinate system is projected, or geographic?

Opening the raster file, I see in properties that Current coordinate system is ETRS89_ETRS_LAEA. But, is is projected or geographic? The EPSG code is missing. Googling ETRS89_ETRS_LAEA revelas EPSG: 3035, but the projection/geographic information is still missing: https://epsg.io/3035

Searching for properties, I again find both projected and geographic CS information:

I have found a python scripting solution here: Determining if spatial reference Geographic or Projected using ArcObjects?

but maybe there is a simpler solution to answer this question without need to run a script? Something like a golden rule, without extra typing, googling, just simply included by ArcGIS engineers?


texts:
I am having trouble with the ArcGIS 10.1 and ERDAS Imagine 9.2 license conflict.
From what I found out it isn't that rare, but I do not understand why there such a conflict in the first place.

Therefore I am looking for some alternatives, heard about PCI's Geomatica but never tried it, any open source solutions?


texts: How to calculate average number from points within a polygon

texts:
I have multiple data points that represent the center of a survey plot with a radius of 5, 10, or 15 feet.  These points also contain percent coverage of different plants.

These points are within different polygons.

Using area of the survey site, I would like to average percentage of this coverage, and then using that estimate the total coverage of each plant across the polygon.

Which tools would I use to calculate this?  I have attempted to use the field calculator but to no avail.


texts: Remove QGIS watermark on QGIS Server tiles

texts:
I am trying to set up QGIS Server in a leaflet client. Everything works fine, except that I get a QGIS logo watermark on each tile, which obviously is not so great.

Here is my WMS request:

http://0.0.0.0/cgi-bin/qgis_mapserv.fcgi?
map=map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetMap&
BBOX=37,-122,38,-121.7&
CRS=EPSG:4326&
WIDTH=722&
HEIGHT=352&
LAYERS=map&
STYLES=default&
FORMAT=image/png&
DPI=96&
TRANSPARENT=true

I have searched everywhere but I can't seem to find where to remove it.

Did anyone face the same problem?


texts: Extracting x,y coordinates from `QgsGeometry: Point(x, y)` using PyQGIS

texts:
I'm trying to extract the coordinates of a feature vertex. When I write on my editor
print(vertex)

what I get is
<QgsGeometry: Point(58746.43000004297, -256665.088999123)>

But when I try to extract the coordinates,
print(vertex.x)

I get

AttributeError: 'QgsGeometry' object has no attribute 'x'

How do I extract the coordinates x,y?
If it is not possible with <QgsGeometry: Point(x,y)>, then how do I simply convert to a form that allows me to extract the coordinates?


texts: I have a grid and some lines within a grid and I was wondering how to measure the length of the lines inside the grid. I can use the measure tool for a straight line, however in the future I will not have straight lines. Is there still a way to find the length of the line? I've tried using the tabulate intersection tool only to receive a "the geometry is not M aware error".


texts: Coordinate values are not displayed correctly in OpenLayers

texts:
I'm trying to display coordinates of map click event. I do not think the coordinates are displayed correctly. The two numbers shown to me are (3889586.1245487547, 6296722.285457358).
Where do you think the problem comes from?
here is my code:
map.on('click',function(e){
        map.forEachFeatureAtPixel(e.pixel,function(feature , layer){
            let clickedCoordinate = e.coordinate;
            let clickedFeatureCountryName = feature.get('NAME_0');
            let clickedFeatureCityName = feature.get('NAME_1');
            let clickedFeatureAreaName = feature.get('NAME_2');
            let clickedLatitudeCoordinate = clickedCoordinate[1];
            let clickedLongitudeCoordinate = clickedCoordinate[0];
            overlayLayer.setPosition(clickedCoordinate);
            overlayFeatureCountryName.innerHTML = clickedFeatureCountryName;
            overlayFeatureCityName.innerHTML = clickedFeatureCityName;
            overlayFeatureAreaName.innerHTML = clickedFeatureAreaName;
            overlayLatitudeCoordinate.innerHTML=clickedLatitudeCoordinate;
            overlayLongitudeCoordinate.innerHTML=clickedLongitudeCoordinate;
        })
    })
}


texts: Generating distance from each address point on map to nearest fire hydrant using ArcGIS Network Analyst?

texts:
I am attempting to measure the distance between one set of points (Address point data) to the closest point from another set of points (Fire hydrants).

When I perform Closest Facility Analysis only a few routes between points show up (and it is not even the closest fire hydrant to that address point).

What am I doing wrong that is preventing me from generating a distance from each address point on my map to the nearest fire hydrant?

Is this a problem with adding the data to the Incidents/ Facilities - what can I do to mitigate/ solve this problem?

A majority of the time it says that it was unable to locate the facility

I am using ArcMap 10.2.2 and have a feature dataset containing all of the feature classes I am using. I can convert to shapefiles if needed.


texts: Missing Contours when extracted from DEM in QGIS?

texts:
I used QGIS Raster->Extraction->Contour tool to extract contours at 1m interval from a virtual raster. The virtual raster was built using around 100 tiles in ASC format that represented a 1m DEM.

When zoomed out the contours looked as I expected conforming to the shape of the terrain very well, even picking up minor detail such as dish drains on the side of roads.

However when I zoomed in I noticed that there was numerous missing contours.  For example, in the image below showing a ridge, the contours step down at 1m increments from about 101m (bottom of image below) to 96m. However then there is a significant jump to the next contour - some examples are 72m in the top right of the image, to 85m on the left of the image.  Looking at the attribute table there also appears to be gaps in data (for instance if I filter display to show contours in the map, there are gaps)

I checked the virtual raster for this area and it looked fine.

I seek advice on how to obtain all the contours as the extracted contours are probably not useful in this state.

Below is a copy of the raster with the contours overlaid. Apologies for the lack of graduation, but you can see different shades of gray showing the terrain.  Due to the size of the raster I had alot of trouble changing the styling.

Below is a copy of the settings I used to extract the contours also.


texts: Move Zoom In & Out Button in ArcGIS JavaScript Web AppBuilder (Developer Edition Download)

texts:
After creating a web app in ArcGIS Online, I downloaded the app using the Web AppBuilder for ArcGIS (Developer Edition). Several folders, .json, .html, .js, and .css files are contained within the download. I have been able to modify the header and website color scheme, as well as move the Home button and various other widgets. However, I cannot find any reference to the Zoom Slider. I have found several posts on how to hide the slider, but not within the downloaded AppBuilder developer edition.

How would I move the slider down the page to accommodate a larger header? Where would this code be inserted (in the config.json, index.html, or another file included within the download)?


texts: Width and value problems with rasterstats

texts:
I am trying to extract data from a raster with a gpkg file as a vector.
The issue is that, the data can't be extracted because the internal Bounding Box is negative. There I am not sure, whether this a bug, or my files are ill formatted.
In detail: I want to extract mean and the medium for the Polygons of the vector for a given raster file.
The CRS of the raster file is: EPSG:25833
the affine transformation is: Affine(1.0, 0.0, 379999.5, 0.0, 1.0, 5795999.5)
I am a computer scientist but I think this transformation is odd? Isn't e  of the affine transformation a negative number?
Example Polygon of the vector has the bounding box: (380670.907, 5796016.979, 380686.087, 5796057.668)
In it's main.py rasterstats calls: rasterio's read method. This one calculates the Bounding Box of the raster file to: ((58, 18), (671, 687)). Which has a negative width and therefore an Error is raised.
I found out when I manipulate the Affine transformation and set affine.e to -affine.e, then this error is not raised, but no statics is computed. This is due to the fact, that when the sub raster is extracted by:
fsrc = rast.read(bounds=geom_bounds)

Calling rasterios read method my the main.py of rsaterstats, this returns a array with only nodata values set.
Is this a bug, or do I have to fix my files somehow? How do I do it?
Full example code:
import numpy as np
from affine import Affine
from rasterstats import zonal_stats
from rioxarray import open_rasterio

raster = open_rasterio(r"path to elevation.xyz")
raster = raster.rio.write_crs(25833)

affine = raster.rio.transform()
affine2 = Affine(a=affine.a, b=affine.b, c=affine.c, d=affine.d, e=-affine.e, f=affine.f)
stats = zonal_stats(r'Path to levelo f detail.gpkg', np.squeeze(slope),
                    affine=affine2, layer="roofsurface", stats=['min', 'max', 'mean', 'median'])


texts: Inserting value to column which is not existing in input using ogr2ogr

texts:
I am inserting geopackage to a table in the PostgreSQL if a certain value is not null with the following comment
ogr2ogr -append -nln pgtablename  -f PostgreSQL "PG:host=xx user=xx dbname=xx password=xx"  E:\workspace\DATA_WORK\data\edit\70319.gpkg Bauwerkslinie -sql "SELECT * from Bauwerkslinie where begehung_von IS NOT NULL"

It works as expected. In the table pgtablename, there is a column that I do not have in gpkg transactionid and I would like to insert a value, that will be always the same for the affected objects in the command above.
Any idea how can I do it?


texts: Dynamic scale dependency from attribute query using ArcGIS for Desktop?

texts:
If I have a feature class called roads with an attribute for type e.g. highway etc.

Is there a way to setup a dynamic scale dependency on the layer based on the type of road?

Like do a SQL query and apply a specific scale dependency range on the  records that match?


texts: Reproject CSV file to different projecton using ogr2ogr

texts:
I have a CSV file "visivisi.csv" with two columns "x,y" which is in my country local projection. The proj4 line is

+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I am trying to use ogr2ogr to convert this CSV to a new CSV file in wgs84 projection, is this even possible? This doesn't work:

ogr2ogr -f csv -s_srs "+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" -t_srs "EPSG:4326" test.csv visivisi.csv

EDIT: Thanks to below answer, this is the solution (too long for comment section):

if you have a csv file with fields "x,y" called "visivisi.csv", make such visivisi.vrt file:

<OGRVRTDataSource>
  <OGRVRTLayer name="visivisi">
    <SrcDataSource>visivisi.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="x" y="y"/>
    <Field name="latitude" src="x" />
    <Field name="longitude" src="y" />
  </OGRVRTLayer>
</OGRVRTDataSource>

Then run this command:

ogr2ogr -overwrite -f CSV -lco GEOMETRY=AS_XY -t_srs EPSG:4326 test.csv visivisi.vrt

test.csv can not exist beforehand. the SRS in the above example is latvian LKS92


texts: How to add loading bar in QGIS Plugin development?

texts:
Every time i use tool or plugin in QGIS, i saw a kind of bar like this for indicating loading progress :

The different between that bar and Progress Bar is that progress bar will increase until our defined iteration (run progress), but that bar only showing up reversibly.
I need to replace it with this code (statusBar) in my main script as alternate of my "please wait.." message :

self.iface.mainWindow().statusBar().showMessage("Please Wait..")
open_set.add(base_coordinate)
....

How to add that bar in my main plugin script (.py)? Is that kind of bar is similar with QProgressBar?


texts: What must be learned to build 3D web map?

texts:
I wish to build a 3D web map based on Silverlight or JavaScript. I cannot build 3D in Silverlight.

What I should do then? Should I learn "Balder" or OpenGL or what?  I must do it before 09/2014. Can you give some ideas?


texts:
I'm making a map package to use it later in a WPF application, I shared the layers as a map package and there is no problem but when I tried to open the map package file again I found that all the layers have the red exclamation mark which means that the link to the data is broken, I fixed it and set the data source then  re shared again but the problem persist.

Here are the pics of the situation:

First before sharing the map package everything is good

Opening the map package file I found the link is broken

Here when I tried to share it again before fixing the data source to see what error will result


texts: Implementing standard deviation classification in QGIS symbology?

texts:
Is it possible to implement Standard deviation classification in QGIS symbology of a layer like in ArcGIS Desktop?
http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Standard_deviation
I need to show values above and below the mean.


texts: Splicing a string and adding a charcter with update cursors issue

texts:
I am new to python, so I am trying to figure these things out on the fly so please bare with me if this is an easy task.  I am using the name of my Feature Classes to populate an attribute field with update cursors, but I am trying to add a dash 2 spaces from the left.  For example the feature class is name OB619 and I need to populate the ID field as OB-619.  I was able to get the update cursor to take the name of the Feature Class, but I haven't been able to figure outnhow to use splicing in the cursor to add the dash.  Any help would be greatly appreciated.

arcpy.env.workspace = r'\\\Public\OB_Test_3.gdb\Test1'

def makeFullPath(path, name):
    return path + "\\" + name

file_name_field = "ID"
fieldname2 = "FLOATABLES"
fieldname3 = "PATHOGENS"
fieldname4 = "NITROGEN"
fieldname5 = "OWNERSHIP"

feature_classes = arcpy.ListFeatureClasses()

fieldsList = [file_name_field, fieldname2, fieldname3, fieldname4, fieldname5 ]

for fc in feature_classes:
    existing_fields = [f.name for f in arcpy.ListFields(fc)]
    if fieldsList not in existing_fields:
       arcpy.management.AddField(fc, file_name_field, 'TEXT', field_length=200)
       arcpy.management.AddField(fc, fieldname2, 'TEXT', field_length=50)
       arcpy.management.AddField(fc, fieldname3, 'TEXT', field_length=50)
       arcpy.management.AddField(fc, fieldname4, 'TEXT', field_length=50)
       arcpy.management.AddField(fc, fieldname5, 'TEXT', field_length=50)
       # write the file name into each row of the file name filed

    with arcpy.da.UpdateCursor(fc, [file_name_field]) as uc:
        for row in uc:
            uc.updateRow([str(fc)])
    del row, uc


texts: Installed QGis on macbook pro. If I needed to reinstall or modify a plugin, where is the plugins folder located?


texts:
What is "self" and how can I use it error-free in a Q script?

Example:

mc = self.iface.mapCanvas()

But when I try to run this in the Console I get the following error:

NameError: name 'self' is not defined

Maybe I need to import a specific module of some kind to use 'self'?


texts: QgsMapRenderer to export current extent using PyQGIS

texts:
I am trying to accomplish PyQGIS simple map rendering via the QGIS 2.0 Python console, but end up with only a blank image. I've followed instructions of this link. My objective is to export a simple image of the current extent of layers of interest to a png. Where can I be going wrong?

#identify crs
crs = QgsCoordinateReferenceSystem(3297, QgsCoordinateReferenceSystem.EpsgCrsId)

# create image
img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)

# set image's background color
color = QColor(255,255,255)
img.fill(color.rgb())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

render = QgsMapRenderer()

#set crs
render.setDestinationCrs(crs)

# set layer set
layers = qgis.utils.iface.legendInterface.layers()
lst = []
for layer in layers:
    lst.append(layer.())

render.setLayerSet(lst)

# set to CURRENT extent**map of layer's full extent is successfully rendered when QgsRectangle(render.fullExtent() is used)
rect = QgsRectangle(render.extent())
rect.scale(1.1)
render.setExtent(rect)

# set output size
render.setOutputSize(img.size(), img.logicalDpiX())

# do the rendering
render.render(p)

p.end()

# save image
img.save("render.png","png")


texts: Showing thousands separator in row count at bottom of ArcMap tables?

texts:
Using Arcmap 10.4.1 is it possible to add thousands separators (such as commas) to the numbers (selected and total) of table records displayed at the bottom of each table display?  I think that readability would be greatly enhanced!

I know how to add thousands separators to the numeric values within the table by right-clicking the field name and modifying its properties, but here I'm interested in the totals displayed at the bottom of each table.

I've looked in the ArcMap Options dropdown, but... nothing.  I would like to make this the default feature.  Any thoughts?


texts:
I want to edit some complex polylines, which I imported from my GPS. These polylines have many vertices and are selfintersecting many times.
When I try to split these lines using the 'Split Features'-command the line is split on every selfintersect and not only where I draw the cuting line. What am I doing wrong? Any suggestions?
I'm using QGIS 1.8 on Windows8 (64).


texts:
I am writing a Python script which uses a vector polygon outline to construct a raster map of a hypothetical 'mound'. Finally, the script will output the volume of this mound using r.volume.

The problem is that the r.volume calculation is resulting in an error. I can't get much info on the error but when I run the volume command in terminal I get no output other than:

:~> Floating point exception: 8

If I run the volume command using python I get:

grass.run_command("r.volume", input=rounded, clump='clumped')
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 395, in run_command
    return handle_errors(returncode, returncode, args, kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 313, in handle_errors
    returncode=returncode)
CalledModuleError: Module run None ['r.volume', 'input=ditchpoly_area_comp_int_round', 'clump=clumped'] ended with error
Process ended with non-zero return code -8. See errors in the (error) output.

Below is the annotated python script that I used. If you want to reproduce the situation then I suggest drawing your own 'mound' shape and then inputting the name of this layer as the first python object and running the script.

import grass.script as grass
import os
import numpy as np

# name of your mound polygon shape layer
geooutline = "ditchpoly_area"

############ mound particulars ############

# distance from edge of polygon to foot of mound
erosion = float(-1.5)

# horizontal distance from mound base to mound peak
baseoffset = float(-3.00)

# height of mound
depth = float(3.30)

# rounding value for mound raster
rastrounding = float(0.1)

# names of mound objects
erooutline = geooutline + "_ero"
erooutlinerast = geooutline + "_ero_rast"
baseoutline = geooutline + "_base"
baseoutlinerast = geooutline + "_base_rast"
surface = geooutline + "_surface"
ditchbase = geooutline + "_ditchbase"
comp = geooutline + "_comp"
comp_int = geooutline + "_comp_int"
rounded = geooutline + "_comp_int_round"
ditchonly = geooutline + "_ditchonly"

############ grass commands ############

# set region
grass.run_command("g.region", flags="a", vectory=geooutline, nsres=1, ewres=1)

# create buffer outlines
grass.run_command("v.buffer", overwrite=True, input=geooutline, output=erooutline, distance=erosion)
grass.run_command("v.buffer", overwrite=True, input=erooutline, output=baseoutline, distance=baseoffset)

# prepare surface raster
grass.run_command("v.to.rast", overwrite=True, input=erooutline, output=erooutlinerast, use='cat')
formula = '"%s" = if(isnull("%s"), 0,null())' % (surface, erooutlinerast)
grass.raster.mapcalc(formula, overwrite=True)

# prepare ditch base
grass.run_command("v.to.rast", overwrite=True, input=baseoutline, output=baseoutlinerast, use='cat')
formula = '{0}" = if({1}" == 1,"{2}",null())' .format(ditchbase, baseoutlinerast, depth)
formula = formula.replace('"', '')
grass.raster.mapcalc(formula, overwrite=True)

# patch composite raster
compformula = surface + ',' + ditchbase
grass.run_command("r.patch", overwrite=True, input=compformula, output=comp)

# fill nulls interpolate ditch profile
grass.run_command("r.fillnulls", overwrite=True, input=comp, output=comp_int, method='rst')

# round values
formula = '"{0}" = float(round("{1}","{2}"))' .format(rounded, comp_int, rastrounding)
formula = formula.replace('"', '')
grass.raster.mapcalc(formula, overwrite=True)

# extract ditch
formula = '"{0}" = if("{1}" == 0,null(),"{2}")' .format(ditchonly, comp_int, rounded)
formula = formula.replace('"', '')
grass.raster.mapcalc(formula, overwrite=True)

# calculate volume
grass.run_command("r.volume", input=rounded, clump=erooutlinerast)

############  remove old files ############

#clear old
grass.run_command("g.remove", flags="f", type='all', pattern=geooutline + "_*")


texts: Using GDAL in python to stack Landsat bands

texts:
I am processing a large amount of landsat data. Landsat scenes ship as a collection of tiffs- one file for each spectral band, plus a bunch of masks. I need to stack some of the bands into a single tiff for subsequent processing. I can easily do this in various programs, but would like to automate the process using GDAL in python. I will eventually need to automate the stacking of hundreds or thousands of images, so efficiency will be important.

I have written a script which is simple and perfectly functional but feels extremely hacky- I open the desired tiffs, convert them to numpy arrays, stack the arrays, and save the stacked array to a new geotiff with metadata taken from the input. Here is a simplified version of the code:

from osgeo import gdal, gdal_array
import numpy as np
b1 = gdal.Open("LT50250232011160PAC01_sr_band1.tif")
b2 = gdal.Open("LT50250232011160PAC01_sr_band2.tif")
array1 = b1.ReadAsArray()
array2 = b2.ReadAsArray()
stacked = np.array([array1,array2])
gdal_array.SaveArray(stacked.astype("int"), "b12_stacked.tif", "GTiff", gdal.Open("LT50250232011160PAC01_sr_band1.tif"))

If possible, I'd like to avoid converting the data to numpy arrays, and just stack the gdal.Dataset objects (b1 and b2). This would make the code cleaner and presumably more efficient. It seems like it should be straightforward to do, but I have not been able to find any information on how to do it.

Alternate approaches which I have come across in research would be to use gdal_merge.py or vrts (as suggested here). Neither of these strikes me as ideal; gdal_merge is designed for tasks much more complicated than mine and presumably involves a lot of overhead, while vrts are only usable in the command line. However, if the approach I outlined above is not feasible, would one of these appraoches be more efficient/pythonic than the hack I have already written?

Edit: my code needs to be run on a Windows system.


texts: OpenGeoSuite/ Geoserver error message "404 status code"

texts:
I am trying out OpenGeosuite in  my WIndows 8 64-bit machine.  OpenGeo Jetty and OpenGeoPostgreSQL services are up and running. Working the Qgis Open geo plug-in I'm trying to publish my layers in Geoserver but I get something like this QGIS error message:

Tried to make a GET request to ttp://localhost:8080/geoserver/rest/workspaces/xx/datastores/xxxx_public/featuretypes/viewdivision.xml but got a 404 status code:

Can someone help me?


texts:
I have a shapefile, which contains track points of a gpx file. Every track point has a number in its "number" field that tells the numeric name of the track the point belongs to. So points are clustered in a way. Unfortunately, there is no information about the order of the points (track segment) in a group.
I have an unclear idea of using Python and a nested For-loop but my Python knowledge is super-limited. This VBA-ish method is the way I think it could work (I'm really sorry, I still have a lot to learn)

For n = 1 To 'max number of groups'
    fcount = 'number of features in a group'
    For m = 1 To fcount
        track_seg = m           #field name where result goes
    Next m
Next n

And this would be the preferred output in QGIS:

(Number one: start point; last number: end point).


texts: Creating AutoCAD-like grading lines in ArcMap?

texts:
I often need to show grading for roads and ditches. This is what the map readers require and what the standard symbology for grading is.  These 3D lines can be easily created in AutoCAD civil 3d and similar software but I would like to, if not create the same, then at least emulate these grading lines as 2D lines or symbology. Ideally, a polygon fill but that's unlikely to work effectively.  I've also tried a number of techniques using hash line symbols and offsets but there are always issues with scale, etc.

Ultimately, I would like to be able to easily create such lines between to irregular lines.  A link to a functioning script, addin, or an ArcObjects code snippet is acceptable.


texts: How to find how grass's r.viewshed works programmatically?

texts:
I want the Python source code of Grass's Viewshed analysis. In QGIS software there is only an .exe  file of that algorithm. Although I'm not sure whether it's written also in python or not.

Where can I find a Python (preferably using pyqgis) source code of that analysis to get more familiar with the coding of that operation?