texts: WMS layer is not overlaying exactly when zooming out. How to add WMS layers to google map layer in openlayers?

texts:
My openlayers code is as follows. But there is shift in between both the layers. How to rectify it?

map = new OpenLayers.Map( 'map_element',{} );
        var gsat = new OpenLayers.Layer.Google(
            "Google Satellite",
            {type: G_HYBRID_MAP,
             numZoomLevels: 22,
             isBaseLayer: true,
             'displayInLayerSwitcher':false},
              {sphericalMercator: true}
            );
        var roads = new OpenLayers.Layer.WMS(
        "Roads", "http://localhost:8088/geoserver/cite/wms",
        {layers: 'cite:viz_roads',
                     isBaseLayer: false,
                     transparent: true,
                     opacity: 1.0}
        );

        var org = new OpenLayers.Layer.WMS(
        "points of interest", "http://localhost:8088/geoserver/cite/wms",
        {layers: 'cite:organisations',
                     isBaseLayer: false,
                     transparent: true, opacity: 1.0}
        );
        map.addLayers([gsat,roads,org]);
        map.addControl(new OpenLayers.Control.LayerSwitcher({}));

if(!map.getCenter()){
map.setCenter(new OpenLayers.LonLat(83.324947,17.732909), 14);
}
    }

This is the snapshot of the map:

Observe the shift from position 1 to 2. 1 is the actual position where the feature should be. It is rectified automatically when I turn layer off and on. But how to rectify this error while zooming?


texts: Performing polygon drilldown/overlay in PostGIS?

texts:
I face a challenge with PostGIS that I cannot seem to wrap my head around. I know I can solve this using a programming language (and that is my backup-plan), but I really do like to solve this in PostGIS. I've tried searching, but could not find any answers that matches my problem, this might be because I'm uncertain on my search terms, so please excuse that and point me in the right direction of there indeed is an answer.

My problem is this:

I have a table with mixed Polygons/MultiPolygons
Each (multi)polygon has an attribute that ranks it (priority)
Each polygon also has a value I would like to know
I have a search area (polygon)
For my query area I want to find the area covered by each polygon value

Example:

Say I have the three polygons depicted in red, green, and indigo here:

And that the smaller, blue rectangle is my query polygon

In addition, the attributes are

geom   | rank | value
-------|------|----
red    |  3   | 0.1
green  |  1   | 0.2
indigo |  2   | 0.2

What I want is to select these geometries, so that the highest ranked (green) fills all the area it can (i.e the intersection between my query geom and that geom), then the next highest (indigo) fills the intersection between the query geom and the geom MINUS the already covered) etc.

Something like this:

I've found this question: Using ST_Difference to remove overlapping features? but it does not seem to do what I want.

I can myself figure out how to compute areas and such, so a query that gives me the three geometries as depicted in the second image is fine!

Additional info:
- This is not a large table (~2000 rows)
- there may be zero or multiple overlaps (not just three)
- there may not be any polygons in my query area (or just in parts of it)
- i'm running postgis 2.3 on postgres 9.6.6

My fallback solution is to do a query like this:

SELECT
ST_Intersection(geom, querygeom) as intersection, rank, value
FROM mytable
WHERE ST_Intersects(geom, querygeom)
ORDER by rank asc

And then iteratively "chop off" parts of the geometries in code. But, as I said, I would really like to do this in PostGIS


texts:
I am trying to load all the SRTM data into a PostGIS database.  I tried the process described here: https://postgis.net/docs/using_raster_dataman.html
It created a giant SQL file.  When I tried loading this file, the process used up all my machine's memory and crashed without loading anything.
How can I load all the SRTM data into a single PostGIS table?
Edit: Here is the command I used to create the SQL file:

raster2pgsql -c *.hgt srtm > srtm.sql

And then I loaded it with:
psql -d public -f srtm.sql

Edit: Here is the error I got:
...(many occurrences of the insert messages)...
INSERT 0 1
psql:srtm.sql:8520: ERROR:  invalid byte sequence for encoding "UTF8": 0xd2 0xd2
psql:srtm.sql:8577: invalid command \
psql:srtm.sql:8604: invalid command \
psql:srtm.sql:9043: out of memory


texts: Zonal statistics mean VS Field statistics mean in ArcGIS

texts:
I calculated mean from a raster for a single polygon using zonal statistics as table tool and got the mean as 440,000.

The raster is the exact same shape as the polygon (used extract by mask tool for this). When I open the attribute table of this raster and do field statistics on the value field, I see that the mean is 330,000.

Shouldn't the two means be the same if not exactly the same?


texts: How do you set the tracking interval in ESRI Workforce?

texts: I am using ESRI Workforce (beta) along with Collector for a Public Works project.  I have created a Workforce App and enabled location tracking.  My problem is the location tracking seems to be taking the GPS reading and plotting a point every minute, but at this interval I end up with way too many points and I keep getting the "Layer Can Not Draw Too Many Points" error.  My question is how do I turn the interval to say 15 minutes or 30 minutes.  I have been in the project settings and in the layer settings for each of my layers and I can not find anywhere to change the ping interval.


texts: Setting CRS of raster to EPSG 3035 using R

texts:
I am trying to set the CRS of a raster file with the following code:
library(raster)
r<-raster(xmn=0,xmx=10,ymn=0,ymx=10)
crs(r)<-CRS("+init=epsg:3035")

The warning message is:
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum European_Terrestrial_Reference_System_1989 in CRS definition

How can I succesfully change the CRS to EPSG 3035?
If I set the CRS to WGS84 I don't get any warning.


texts: How can I disable logs being created in standalone python script?

texts:
When I run the following script in QGIS 2.6.1, it produces a folder called 'processing', which contains a log, in the directory it reads files from. And since the script reads files from several locations, there's going to be 'processing' folders scattered around.

Is there a way to disable this?

I've excluded the processes in the script as I believe this is the cause:

from processing.ProcessingPlugin import ProcessingPlugin
plugin = ProcessingPlugin(iface)
from processing.tools import *

Here is the script:

import os, sys, glob

from qgis.core import *
from qgis.gui import QgsMapCanvas
from PyQt4.QtGui import *

from os.path import expanduser
home = expanduser("~")

QgsApplication( [], False, home + "\AppData\Local\Temp" ) # Set a custom config path

QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)
QgsApplication.initQgis()
app = QApplication([])

#Run processes
def run():
...


texts: Rounding coordinates to 5 decimals in GeoPandas?

texts:
I have a shapefile, which I converted into geopandas dataframe. I would like to round the coordinates to 5 decimals.

Is there a method/function in geopandas that can do this?

My shapefile contains multipolygons.


texts: How to convert Shapefile geometries to WKB using OGR?

texts:
I got the geometry of a feature of the shapefile, and i want to save that geometry in the postgis (in WKB formats like it happens when we import shapefiles using shp2pgsql and psql commands). How do I convert then

To get Geometry i have used OSGeo OGR library
eg:

feat = layer.GetFeature(0)
geometry = feat.GetGeometryRef()

and I have got

<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0x0096A2D8> >

so how will convert it into the WKB geometry?
I am using Python for this.


texts:
I'm looking for some package in R to estimate the value of a variable by moving window regression. I have read in Local Models for Spatial Analysis, which is a good method to predict values at unsampled locations. A very similar case is mentioned in the book, PRISM.

I've been reading various R packages, but have not found anything about MWR, GWR only.

Is there another way of naming or can GWR be adapted to function as MWR?

Text from book ->
Local Models for Spatial Analysis


texts: How to sum up all field values in QGIS modeler

texts:
I am building a model to generate two zones around a point, a 250m-buffer, and a ring of 250 to 1000m, where i recalculate values from an other layer according to the intersections share.

The purpose is to examine if values in the 1. buffer are smaller than in the second.
The rest of explanations would take too long and everything goes fine so far.

The final step would now be, to sum up the values for each layer, but somehow I can't find a solution.

There is the "Statistics by categories" function in the modeler which produces a stats-table i could join back to my layers and then dissolve them, but ... for whatever reason ... this only produces values for min, max, mean and stddev, but no sum of values!!

All dissolve functions available, do not offer any summing or whatever actions.

Anyone got any idea?

I'm stuck 1 meter in front of the finish line ...

P.s: there is NO option to do anything manually there. Simply too many points with too many other parameters to check each time. Thats why I made the model ...

UPDATE:
Now I had the ideas to create centroids from all intersected polygons and then do a SAGA Points in Polygon statistic (which includes a SUM !!!!) with the un-intersected buffer a few steps before ... but this time the SAGA function is buggy and doesn't work.


texts: Dimensions of images given their latitude and longitude for Google Earth Engine

texts:
From the point

  long: 16.499050771335575, lat: -30.882

I want to find the corresponding latitude and longitude values that would result in a 200 km by 200 km square in which the aforementioned point would be the bottom left (southwest) corner of the square.

I stumbled upon Understanding terms in Length of Degree formula? that gave the corresponding formulas. So I used this to come up with the points of my square:

  16.499050771335575,-29.07796275949132

  18.584367695114103,-29.07796275949132

  18.584367695114103,-30.882

  16.499050771335575,-30.882

However, when using the following code to download the corresponding image from the Modis NDVI 1km data set through Google Earth Engine, I got a 233 by 201 image.

path = img.getDownloadUrl({
    'scale': 1000,
    'bands': [{'id':'NDVI'}],
    'region': '[[16.499050771335575,-29.07796275949132],[18.584367695114103,-29.07796275949132],[18.584367695114103,-30.882],[16.499050771335575,-30.882],[16.499050771335575,-29.07796275949132]]'
})

I have double checked that the distances between the four points are indeed ~200 km.

I don't understand why the code is producing a 233 by 201 km image instead of the expected 200 by 200 km image?

Furthermore, when I moved the latitude up to the equator, I got an image returned by Google Earth engine that was 201 by 203, even though the distances between the points was ~200 (plus minus 1).

Update: After further investigation I've noticed that the distance between two neighboring pixels of the same latitude is 858 meters and not 1000 meters. This is odd, as I specified in the code the spatial resolution of 1000 meters.


texts: pass an array of geometry to postgresql function

texts:
I am trying to pass an array of geometry to a postgis function using JDBC with JAVA:

            CallableStatement proc = conn.prepareCall("{ ? = call somefunc(?) }");
            proc.registerOutParameter(1, 7);
            PGgeometry[] arr = new PGgeometry[2];
            arr[0] = a;
            arr[1] = b;
            proc.setArray(2, conn.createArrayOf("geometry", arr));
            proc.execute();

a,b are valid PGgeometry types.
I am getting this exception:

Exception in thread "main" java.lang.AbstractMethodError: org.postgresql.jdbc3g.Jdbc3gConnection.createArrayOf(Ljava/lang/String;[Ljava/lang/Object;)Ljava/sql/Array;

Is this even possible to do? I found no good answer for that.

Thanks!


texts:
I am using OGR to read features from a shape file; when the feature is destroyed I am getting an error 0xC0000005: Access violation reading location 0xffffffffffffffff.

This secton of the code looks pointless but is a simplification of a larger section showing where the error is occuring. I have removed lines of code that could interfere with the objects and I am still getting the same error message.

OGRDataSource *hDS;
OGRSFDriver   *Driver;
OGRFeature    *pFeat;

hDS = OGRSFDriverRegistrar::Open(BasePath,FALSE,&Driver);
OGRLayer *Layer = hDS->GetLayerByName(ThisLayerName);
Layer->ResetReading();

while ((pFeat = Layer->GetNextFeature()) != NULL)
{
    // I did intend to do something here
    OGRFeature::DestroyFeature(pFeat); // error occurs here
}
OGRDataSource::DestroyDataSource(hDS);

I am using OGR2OGR source as a reference on how to read an OGR datasource. I can see in the reference that each feature returned is destroyed before grabbing a new one until the last feature and the datasource is destroyed thereafter.

I am using GDAL version 1.10 64bit on Windows 7 64bit; I have tried previous versions 1.9.2 and 1.8 and they all return the same error. I am using Visual C++ if that makes any difference.


texts: Combining multiple mosaics into single one using ArcMap?

texts:
I collected 615 images with Phantom 4dji. My computer couldn't process all of them together on Agisoft Photoscan, so I separated my photos by area, making five orthophotos.
My goal is to combine all of these 5 orthopotos into a single one to make a mosaic of the hole area.
First I was thinking on making it on Agisoft but then I got an error: "SSL support not available". So, I tried in Arcmap following this tutorial, otherwise my mosaic ends up having black gaps between the orthophotos:

How can I solve this issue?


texts: Creating visually standardized RGB images from Sentinel-2 data

texts:
I want to create "visually appealing" RGB images from Sentinel-2 L1C data in a standardized manner. Obviously, this starts with extracting the bands 4, 3, and 2, which correspond to the Red, Green, and Blue channels of the electromagnetic spectrum. The problem, however, is that the histograms of the individual bands in the different images vary significantly in both shape and min/max values. So I am wondering what is the best way to create RGB images that follow some kind of standardized visual appearance.

Is there any fixed conversion factors or thresholds? Or would it be more advisable to use relative thresholds, e.g. by keeping the middle 95% of the gray values within each band?


texts: Adding a shapefile to a feature layer using arcgis js api

texts:
I have the following code and up until a couple of weeks ago (or possibly one or two months because I hadn't used it in some time) everything seemed to be working fine:

function generateFeatureCollection(fileName) {
    require(["esri/geometry/scaleUtils", "esri/request", "dojo/dom", "dojo/on", "dojo/_base/lang"], function (scaleUtils, request, dom, on, lang) {
        var portalUrl = "http://www.arcgis.com";
        app.map.setMapCursor("progress");
        document.body.style.cursor = 'progress';
        //resolving the fakepath issue in the file location url
        var name = fileName.split(".");
        name = name[0].replace("c:\\fakepath\\", "");

        var params = {
            'name': name,
            'targetSR': app.map.spatialReference,
            'maxRecordCount': 12500,
            'enforceInputFileSizeLimit': true,
            'enforceOutputJsonSizeLimit': true
        };

        var extent = scaleUtils.getExtentForScale(app.map, 40000);
        var resolution = extent.getWidth() / app.map.width;
        params.generalize = true;
        params.maxAllowableOffset = resolution;
        params.reducePrecision = true;
        params.numberOfDigitsAfterDecimal = 0;

        var myContent = {
            'filetype': 'shapefile',
            'publishParameters': JSON.stringify(params),
            'f': 'json',
            'callback.html': 'textarea'
        };

        var genRequest = request({
            url: portalUrl + '/sharing/rest/content/features/generate',
            content: myContent,
            form: dom.byId('uploadForm'),
            handleAs: 'json',
            load: lang.hitch(this, function (response) {
                if (response.error) {
                    return;
                }
                addShapefileToMap(response.featureCollection);
            }),
        });
        genRequest.then(
            function (response) {
                //console.log("Success: ", response.layers);
            }, function (error) {
                app.notifyError("Error on upload.");
                $('#IEFill').append("<b>Upload Error: "+error+"</b>");
            });
    });
}

However recently an error is being thrown in the last part of the code where the line is:

app.notifyError("Error on upload.");

The following error is being thrown and I have been through all my code and have no idea what I may be doing wrong especially since this process was working and no changes were made to the best of my knowledge.

SyntaxError: Unexpected token <
at Object.b.fromJson (init.js:230)
at Object.json (init.js:187)
at s (init.js:192)
at c (init.js:76)
at d (init.js:76)
at b.Deferred.resolve.callback (init.js:77)
at init.js:195
at k (init.js:199)
at m (init.js:199)
at f.resolve (init.js:201)

Would anyone know why this error is being thrown and any possible solutions for it?


texts: Editing list of strings using Google Earth Engine?

texts:
I'm trying to edit string elements in a list object. I need to get a list of IDs of a image collection, and then replace the orbit row to another (in this case, 073 to 074). I'm new to GEE and got stucked creating the .map() function. Does anyone know how can I solve this?

var landsat = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
    .filterDate('2008-01-01', '2011-12-31')
    .filter(ee.Filter.calendarRange(6,7,'month'))
    .filter(ee.Filter.eq('WRS_PATH', 222))
    .filter(ee.Filter.eq('WRS_ROW', 73))
    .filterMetadata('CLOUD_COVER', 'equals', 0)

var names = ee.List(landsat.aggregate_array('system:index'));
print(names)

var replace = function(i){
  var size = (i.size().subtract(1)) // get list size
  var imgs = ee.List.sequence(0, 11) // creates list sequence to iterate
  var replace = imgs.map(function(j) {
    var str = ee.String(i.get(j)) // convert each element to a string
    var str2 = str.replace('222073','222074')
    var newlist = []
    return newlist.add(str2)
  })
}
var newnames = names.map(replace)

What I have:

LT05_222073_20080604
LT05_222073_20080706
LT05_222073_20080722
LT05_222073_20100626
LT05_222073_20110715

What I need:

LT05_222074_20080604
LT05_222074_20080706
LT05_222074_20080722
LT05_222074_20100626
LT05_222074_20110715


texts: Formally testing whether building type tends to be close to another building type?

texts:
In arcGIS 10.1 (advanced licence), I have two polygon shapefiles, each representing two different types of building (let's call them A and B).

Each type is actually used for two different purposes. What I would like to accomplish is to formally test if building type B tends to "occur" close to building type A.

What procedure, either within arcGIS or with other programs, could I follow?

As for 'other programs', kindly note that I have some familiarity with SAM, GeoDA, and Passage programs, and that I have no experience with GRASS.


texts: How to find the nearest land mark from a given longitude and latitude?

texts:
I have a table called osm_places with place_name geometry, place_type and many other info related to the city. My database system is PostGIS.

My requirement is:

Given a latitude and longitude (which may or many not be the exact lat and long stored in the database), I need to find the nearest city or place_name  from the table osm_places with the help of the input latitude and longitude. Any help is highly appreciated.

Is there a nearest neighbor implementation in php available for this or SQL to handle this?


texts: Finding accurate distance of a drawing in QGIS

texts:
Attached shows a rough drawn image of a pipeline track. I need to measure the pipeline track from one end to the other, e.g from the red and blue square to the other red and blue square, e.g from the blue small square to the red circle.

The issue I am facing is that I have 100's of pipeline tracks to measure.
I know using the distance to hubtool will give me a length but not follow the shape of the pipeline, rather will give me a distance in a straight line.

I have also used the measure tool. This tool works but it will take forever to measure 100's of pipeline tracks.

Is there another way which can accurately and efficiently measure the distances? Is that this drawing has no coordinates:


texts: QGIS 2.18 Training Manual Exercise 1.2.1 - Query Builder Problem

texts:
For Exercise 1.2.1 in the QGIS Training Manual Training Manual Link I attempted to extract data from multipolygon layer using the query;

building != 'NULL'

I had double checked if there were features with no information in their 'building' field in my data and there were many. Unfortunately the query failed with message shown in image below.

I have tried to exclude spaces, replacing != with NOT, different data, and using " " instead of "NULL"

As I am a bit of a rookie with coding, I have been having no luck.

No luck with Kazuhito Answer, "building" != 'NULL'. This query isn't selecting any features, please see photos for reference.


texts: How to calculate the terrain profile in Profile Tools plugin

texts:
Now, I'm trying to calculate terrain slope(°) using pyqgis using below process.

Along the specific lines which I want to know the terrain profile, I generated the points using 'qgis:pointsalonglines'.

line_lyr_1 = processing.run('native:pointsalonglines',{'INPUT':line_lyr_0,'DISTANCE':10, 'START_OFFSET':0, 'END_OFFSET':0, 'OUTPUT':'memory:'})['OUTPUT']

Using generated points and DEM raster layer, I could get the elevation data for each points by below codes.

elevation = namedtuple('elevation', ['Place_Name','Line_ID','Distance', 'Elevation'])
elevation_datas = []

for feat in line_lyr_1.getFeatures():
    point = []
    tanada_ID = []
    line_ID = []
    distance = []
    elevation_point = []

    point = feat.geometry().asPoint()
    Place_ID = feat.attributes()[0]
    Line_ID =  feat.attributes()[1]
    Distance = feat.attributes()[2]
    Elevation = list(rst_lyr_1.dataProvider().identify(point, QgsRaster.IdentifyFormatValue).results().values())[0]

    elevation_datas.append(elevation(Place_ID, Line_ID, Distance, Elevation))

After the exporting the elevation data, I made a terrain plots using exported data with the equation calculating slope degree atan(Elevation/Distance)*180/pi() in Excel.

However, using the profile tool plugin in QGIS, I got below results.

Why did my result become different between two plots?
Is the calculation method different in profile tool plugin?
I think the result of profile tool plugin is correct.


texts: GeoPandas geometry is different on Windows and Linux

texts:
Building a same application both on Linux and Windows I noticed that the geometry is handled differently after transferring the dataframe from WGS84 projection to OSM.

The results of the program below are different when run on Windows or when run on Linux with the x and y values interchanged:

import geopandas as gpd
from shapely.geometry import Point

EPSG_OSM = 3857
EPSG_WGS84 = 4326

longitude = 45
latitude = 0
my_location = gpd.GeoDataFrame(geometry=[Point(longitude, latitude)])

geom = my_location.loc[0].geometry
print(f'long: {geom.x:15.4f}, lat: {geom.y:15.4f}')

my_location.crs = EPSG_WGS84
my_location = my_location.to_crs(epsg=EPSG_OSM)

geom = my_location.loc[0].geometry
print(f'x   : {geom.x:15.4f}, y  : {geom.y:15.4f}')

Linux answer:

long:          0.0000, lat:        -45.0000
x   :          0.0000, y  :   -5621521.4862

long:         45.0000, lat:          0.0000
x   :    5009377.0857, y  :          0.0000

Windows answer:

long:          0.0000, lat:        -45.0000
x   :   -5009377.0857, y  :          0.0000

long:         45.0000, lat:          0.0000
x   :          0.0000, y  :    5621521.4862

Is this correct behaviour? The Windows results seem incorrect.

I can resolve the matter by switching the latitude and longitude in Windows..., but looks to me there is something not right!


texts: How to write every raster from a list of rasters?

texts:
I have created a list of rasters - "norm" and run a several loops. Now I would like to write every created raster from "norm" (containing 1 layer) into one raster, having its specific name.

# I can get the rasters names
rst.names<-c(1:length(norm))
for(i in 1:length(norm)){
    rst.names1<-names(norm[[i]])
    rst.names[[i]]<-rst.names1
    }

To write and save a rasters I have found a very useful answer on https://stackoverflow.com/questions/14890369/how-to-write-rasters-after-stacking-them

where they use to respond my question (modified):

unst_norm<- unstack(norm)
for(i in 1:lenght(unst_norm)){writeRaster(unst_norm[[i]], file=rst.names[i])}

However if I use "unstack" I have this error:

Error in stats::as.formula(form) :
argument "form" is missing, with no default

and I can´t find any help for that.

I can use  "unlist" instead of "unstack" but it doesn't seem to work...

but I still get an error:

Error in .local(.Object, ...) : Unable to create dataset

Please, any ideas where the problem could be? Or, how to automatically write and save all rasters from raster list?


texts: Keep legend item & feature count even when item does not appear in current Atlas feature

texts:
I would like to see all legend items  irrespective of whether the item appears in the current Atlas feature. I have formatted the legend such that count of each item appears in brackets, next to the item label in the legend. For items that do not appear in the current Atlas feature, I would like to see a zero in brackets, next to the item label in the legend.

I however would like the legend to reflect only information regarding the current Atlas Feature. For example: if for a given Atlas feature, there are 10 Physician offices, 7 Diagnostic Imaging Facilities, 0 Acute care Facilities, 0 Ambulatory care facilities, I would like the legend to look as follows:

Physician Offices (10)
Diagnostic Imaging Facility (7)
Acute Care   Facility (0)
Ambulatory Care Facility (0)


texts: OpenLayers adding/removing layer when checkbox is changed

texts:
I am trying to add and remove layers2 layer when checkbox is changed. How can I do it with jquery?

Following my way  i get: Uncaught ReferenceError: layers is not defined

     var wmsLayer2 = new ol.layer.Image({
          source: wmsSource2
        });

        $("#Buildings_chkbx").change(function(){
   var bool = !bool;
    if(bool)
          layers.push(wmsLayer2);
    else
    //remove layer...
        });

        var map = new ol.Map({
          target: 'map',
          layers: [baseLayer, wmsLayer],
          view: new ol.View({
            center: ol.proj.fromLonLat([33.33386, 35.14710]),
            zoom: 14
          })
        });


texts: Reading geojson properties in openlayers 4?

texts:
My code is very simple. This is my vector-source and this is how I render my geojson :

var vectorSource = new ol.source.Vector({
   format: new ol.format.GeoJSON(),
   url: 'some_geojson_is_here'
});

Then I pass this vectorSource as a source to a vector layer etc etc. Now I'm wondering is there any way I can get every feature this geojson has??

In this link : https://openlayers.org/en/latest/apidoc/ol.format.GeoJSON.html the 3-rd method it's what I think the way to get all features but I couldn't figure out how to use it.

Any ideas?


texts: Adding new Field and writing values for all features

texts:
I need to create new Field in Attribute table of vector layer, and afterwards to calculate values and store them in that Field for all features I have in that layer.
I have solution which works, but thing is that it is too slow. Does anyone has an idea for faster way?

    # Dodavanje nove kolone
    name = "Vrednost"
    provider = trecilayer.dataProvider()
    caps = provider.capabilities()
    # Provera da li kolona vec postoji, vraca "-1" ako ne postoji
    ind = provider.fieldNameIndex(name)
    try:
        if ind == -1:
            if caps & QgsVectorDataProvider.AddAttributes:
                res = provider.addAttributes( [ QgsField(name,QVariant.Double) ] )
    except:
        return False
    #trecilayer.updateFields()

    # Zamisljam da sam koeficijente ucitao iz ulaznog fajla
    #koef=[1.00, 0.85, 0.80, 0.75, 0.70, 0.60, 0.50, 0.40, 0.00]
    lista=open(self.putanja,'r')
    lista.read
    for string in lista:
        x=string.split(",")
    koef=[]
    for i in range(len(x)):
        koef.append(float(x[i]))
    lista.close()

    #priprema za progress bar
    #ukupan broj redova
    trecilayer.selectAll()
    count = trecilayer.selectedFeatureCount()
    i=0
    #racunanje vrednosti
    trecilayer.startEditing()
    for feature in trecilayer.getFeatures():
        QgsMessageLog.logMessage(str(trecilayer.fieldNameIndex('Vrednost')), 'Vrednost_Parcele')
        razred=feature['Procembeni']
        povrsina=feature.geometry().area()
        nova=povrsina*float(koef[int(razred-1)])
        feature['Vrednost'] = nova
        trecilayer.updateFeature(feature)
        #progress bar
        self.pb(i,count)
        i=i+1

    trecilayer.commitChanges()


texts:
I made a script that makes a chart of the NDVI values by dates for different sites.
When I create an APP to share that chart, the only thing visualized is the MAP.

How do I create an APP with only (if possible) a Chart?


texts: Compass error (can't find proper directories)

texts:
Compass::SpriteException in Sessions#new

  Showing /home/cartodb/cartodb/app/views/layouts/front_layout.html.erb
  where line #12 raised:

  No files were found in the load path matching "icons/*.png". Your
  current load paths are: /home/cartodb/cartodb/app/assets/images   (in
  /home/cartodb/cartodb/app/assets/stylesheets/common/header.css.scss)

I've googled and done this: "git submodule foreach git checkout master", and then "rvm use 1.9.2@cartodb --create && bundle install"….
I can get to the next error (and the next) by manually adding various png files but this is not ideal.

Any suggestions?


texts: Travel planning - Finding farthest points of a polygon

texts:
I’d like to organize a tour around the lake with 2 objectives:

maximise the length of shores to view
minimise travel distance

In order to do so I can:

Pick 2 most remote points on the lake in terms of travelling distance
Define ‘left’ and ‘right’ shores of the lake
Travel at equal distance between shores from A to B

What is the best way to find these 2 points?

UPDATE ON ORIGINAL Q:

The approach I am using is to place multiple random points and calculate Euclidean minimum spanning tree.

Iterate trough all pairs of points/nodes in order to find a couple with longest distance. Picture shows distance where destination is point at the North.
It works but a bit cumbersome, so I am looking for ideas. Raster or vector will do.

This is about never ending discussion re polygon length.


texts: Manipulating Azimuthal Equidistant Projections in QGIS

texts:
I have a shapefile of world countries that is projected as Azimuthal Equidistant with Chicago, USA as the center. I would like to change the projection such that Edinburgh, UK is the center. One would think that this should be as simple as changing 2 numbers in the coordinate system parameters...

The problem that I seem to be experiencing is that QGIS makes it difficult if not impossible to view the complete proj4 code of a given shapefile's coordinate system parameters.

This is what I have tried so far:

Added shapefile of world countries projected as Azimuthal Equidistant with Chicago as the center
Opened properties, metadata tab, copy proj4 parameter code
Went to settings, custom projection, paste in proj4 code, changed the lat/lon from that of Chicago to that of Edinburgh
Saved shapefile of world countries, specified my new custom projection, added new file to map
My new file looks like a map of Pangaea, with the continents all smushed together.

I think there must be other parameters aside from those shown in the proj4 code of the metadata tab. The only thing I'm changing in the proj4 code is the center of the map, but the change I'm seeing is much more than that. My observation suggests that more than one independent variable is changing.

I would be happy to send my shapefile to anyone who would like to play with it.

Thanks,

Daniel Wolf
Environmental/Geospatial Enthusiast


texts:
OS : osx 10.9.2
* I just operated a clean osx installation.
Globe have never worked since I started to learn QGIS.
When I tried globe, it terminated QGIS then a message(see below) occurred.

Any Idea?

Thanks.


texts:
Basically I generate random numbers in R language which are indeed times between bus arrivals.
I need to simulate a moving point(a bus) in OpenStreetMap by using this generated numbers.
I have read that R language communicates very good with OSM but not sure how to realize this goal.
Any idea?


texts: How do I add a Corel Draw map as .dxf file to QGIS and how do I georeference it?

texts:
I have a vector data map of Poland, which I exported from Corel Draw. I can display the map in QGIS. At first the boundaries and subdivisions are only displayed as lines, after saving the data as shapefile, I was able to use the "Line to Polygon"-tool to convert the lines. Now I have a shapefile with multiple features.

Now I want to relocate and distort it to make it fit to my basemap. As you can see in the picture, the Poland shapefile is presented close to Africa because the location is not yet defined correctly.

There are now several problems:

I tried to move the Poland shapefile in an edit session but the shapefile contains of multiple polygons and I was not able to move them together. How can I move multiple features at the same time?
I also tried to make a single polygon by merging the multiple features in an edit session, but when doing it with the "Merge Selected Features"-tool, an error occurs and QGIS crashes.
Therefore, I tried it with a third way by creating a new shapefile, drawing a polygon, that covers the whole Poland polygons and then cliping both shapefiles. After Cliping, the new shapefile is created and added to the table of contents, but there is nothing displayed on my map.

As you can see, there may be a mistake in the previous steps but I dont know where, I already repeated every step but there is still no solution found. I suspect the problem to be in the step of converting the Corel Draw .dxf file to a shapefile or the transformation from lines to polygons.

Also if you know a more effective way to get my map to its planned location and how to reference it correctly, please let me know your thoughts.


texts: Connect OpenLayers3 Client with QGIS Server

texts:
at the moment I am developing a map editor on a  web.

After researching at the internet about map clients, I decided to use OpenLayers3 client because it is very simple an has all the tools i need for my project.

So, my question is:

 - It is possible to connect OpenLayers3 Client with QGIS Server?


texts: How to set parameter filters in ToolValidator?

texts:
I am populating parameter[2] using validation code provided by ESRI. How can I modify this code to repeat for parameter[3] what was done for parameter[2]? To elaborate...param[0] is a user supplied  feature class.  param[1] is fields from provided feature class. param[2] provides sorted values from selected field, and i would like to provide and sort these same values again in param[3]...

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
     """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()
    self.fcfield = (None, None)

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parmater
    has been changed."""
    if self.params[0].value and self.params[1].value:
      fc, col = str(self.params[0].value), str(self.params[1].value)
      if self.fcfield != (fc, col):
        self.fcfield = (fc, col)
        self.params[2].filter.list = [str(val) for val in
                                    sorted(
                                      set(
                                        row.getValue(col)
                                        for row in arcpy.SearchCursor(fc, None, None,
                                                                      col)
                                          )
                                       )
                                  ]
    if self.params[2].value not in self.params[2].filter.list:
      self.params[2].value = self.params[2].filter.list[0]

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""


texts: PostGIS ERROR: function x(geometry) does not exist

texts:
I am using the following postgis query for updating a pgrouting table from my custom road data. In the x(ST_StartPoint(part_geom)) as x1, line its returning an error mentioned

  ERROR:  function x(geometry) does not exist

INSERT INTO public.roads_r (
SELECT DISTINCT ON (the_geom)
nextval('public.roads_r_gid_seq') AS gid,
gid AS parent_gid,
name,
no,
type,
NULL as source,
NULL as target,
ST_Length(part_geom, true) as cost,
ST_Length(part_geom, true) as reverse_cost,
x(ST_StartPoint(part_geom)) as x1,
y(ST_StartPoint(part_geom)) as y1,
x(ST_EndPoint(part_geom)) as x2,
y(ST_EndPoint(part_geom)) as y2,
NULL as to_cost,
NULL as rule,
ST_Multi(part_geom) AS the_geom
FROM (
    SELECT
        a.gid,
        a.geom,
        (ST_Dump(ST_Union(b.geom))).geom AS part_geom,
        a.road_name as name,
        a.road_no as no,
        a.road_type as type
    FROM public.roads as a, public.roads as b
    WHERE ST_Intersects(a.geom, b.geom)
    GROUP BY a.gid, a.geom, a.road_name, a.road_no, a.road_type
) AS res
WHERE ST_Contains(geom, part_geom) = true OR ST_Contains(ST_Buffer(geom, 0.00001), part_geom) = true)


texts: Sending data between QGIS and other programs in real time

texts: I would like to be able to grab information about QGIS users' interactions and send things like cursor position and button events to a separate program running concurrently in MacOS. With GRASS tools I could use the Unix Pipe for this. Is it possible to do a similar thing with QGIS? Or is there a better way to get information from QGIS to a concurrent program? (We'll be writing a Python plug-in to capture the data within QGIS.)


texts: How to upgrade PostGIS in OSGEO-Live Xubuntu virtual machine?

texts:
I have downlaoded the OSGEO-Live 6.0 virtual machine running Linux Xubuntu.  I have it running through VirtualBox on a Windows machine.  It currently comes preinstalled with a lot of GIS software, including PostGIS 1.5.  I would like to upgrade PostGIS to 2.x.

Having not used linux much, what is the best method to upgrade PostGIS from 1.5 to 2.x on the Xubuntu virtual machine?  Since I just downloaded the Live-DVD, we can assume there is no data that needs to be retained.


texts:
I am writting an addin using ArcObjects C#. I have setup the following listeners, but they are not being hit.

My code at the top of the addin.Extension class is:

...
private static IActiveViewEvents_ItemAddedEventHandler tocAddHandler;
private static IActiveViewEvents_ItemDeletedEventHandler tocDeleteHandler;
...

protected override void OnStartup()
{
     WireDocumentEvents();
}
....
 private void WireDocumentEvents()
{
ArcMap.Events.OpenDocument +=
                new ESRI.ArcGIS.ArcMapUI.IDocumentEvents_OpenDocumentEventHandler(PostMxDLoad);

}
...
private void PostMxDLoad()
{
... some other stuff

 tocDeleteHandler = new IActiveViewEvents_ItemDeletedEventHandler(SettingsMgr.TableofContents.OnTOCItemDeleted);
((IActiveViewEvents_Event)ArcMap.Document.FocusMap).ItemDeleted += tocDeleteHandler;

tocAddHandler = new IActiveViewEvents_ItemAddedEventHandler(SettingsMgr.TableofContents.OnTOCItemAdded);
((IActiveViewEvents_Event)ArcMap.Document.FocusMap).ItemAdded += tocAddHandler;
}


texts: ModelBuilder overwrites output when using split tool?

texts:
I have one SHP file containing polygons that describe splitlines in rows and cols.
I'd like to use this to split an existing dataset of a few hundred SHPs containing tree positioning data. The extends of the tree SHP files do not match those of the split SHP.
It is therefore likely that contents from various SHPs will (should) be written in one of the split tiles SHPs.

My setup in Modelbuilder is as follows:

Problem now is that the split tool overwrites previously created SHPs when when moving on to new iterations of the tree SHPs.

"Collect Values" doesn't work as the resulting files are named by the split field attribute which I can't pull from the split tool.

Same applies for inline variables as described here.
Since I cannot pull the name of the created SHPs from the split tool.

I also checked Geoprocessing options and overwrite is turned off.

Any ideas how to avoid the split tool to overwrite duplicates?


texts: QGIS expression for label text to display vertex coordinates of polygons

texts:
This question is inspired by Displaying vertex coordinates of a polygon or line without creating a new layer and especially by the answer by @JGH: https://gis.stackexchange.com/a/381130/88814
The basic problem is how to label vertices in QGIS using expressions only (without creating a new layer). Whereas @JGH uses geometry generator to create points on the vertices and sets a Font Marker to it, which in turn is not a character, but an expression to generate in this way sort of a label for the coordinates of every vertex.
This is really a great solution, only that you don't have all the styling options for labels like marker lines, rotation, text buffer etc. (see screenshot below). It would be great to create real labels for each vertex. The only problem seems to be that it is not possible to select the @geometry_part_num Variable (as in @JGH's solution) in the label text field as this variable is not available.
So the question is: is it possible to find an expression that generates the values of the x and y coordinates of each vertex?
Placing a label for each vertex of a geometry (but with the same label text) is possible by setting a label that is placed on a point created by geometry generator with the following expression - maybe it can be adapted somehow for the label text as well:
collect_geometries(
    array_foreach (
        generate_series( 1,  num_points( $geometry)),

        make_point (
            x ( point_n($geometry, @element)),
            y ( point_n($geometry, @element))
        )
    )
)


texts: In ArcGIS Web appbuilder widget, how can i use a variable what is declared in another widget file?

texts:
In Edit widget, i modified the RelatedRecordsEditor.js to automatically upload input fields with data like webMapId.

        var mapid=dijit.byId("dijit_form_ValidationTextBox_1").attr("value");
        console.log(mapid);
        var webMapId = this.map.itemInfo.item.id; //this is not working, only in widget.js
        var update = mapid.replace(mapid, webMapId);
        console.log(res);
        dojo.byId('dijit_form_ValidationTextBox_1').value = update;
      });

But there is no 'map' declared in this file, only in widget.js. How Can I reach that variable from RelatedRecordsEditor.js? Or how can I declare it? I tried:

define(['esri/map'] function(map)
        var map = this.map;
        console.log(map);
--->undefined

Or how can i use widges.js's map object in this js?


texts: Concatenate multiple rows into one where the NULL values are without losing the other information?

texts:
I am wondering if it is possible to concatenate multiple rows in an ArcMap shapefile attribute table into one row based on shape length.

In the table below, I have created the fields "Form B1, Form B2, Form B3, Material B1, Material B2, Material B3". I have populated those values with information from the Zone, Component, Form and Material columns.

All 3 rows have identical geometry and are overlaid on top of each other. I want to combine them into one single row without losing the data.

For example, "Form B1" has been populated with information from:

Zone B
Component 1
Form
Material

"Form B2" with information from:

Zone B
Component 2
Form
Material

So my question is, can I merge/concatenate these 3 rows into one, and have the information fill in where the NULL value current is?


texts: Displaying estimated computational time when running tools in ArcGIS for Desktop?

texts:
I'm almost positive I have set up models to display the estimated runtime/computational time for the models/tools in the past, but I can't seem remember how to turn on this option. I've looked around the web to no avail and would like to determine if this is even possible.

If it is not possible to just turn on the option is there python code I can use to get this information?


texts: Root Mean Square Error of 2 raster layer in QGIS

texts:
Is there a method to calculate the RMSE between 2 raster layers in QGIS?

I have 2 point cloud files, 1 edited and the other not. After triangulation of each layer I have obtained 2 raster files.

I would like to find the RMSE in terms of their vertical height diff..

Is it possible to subtract one raster from the other and form a new raster and from export and calculate the RMSE in excel based on the z value?


texts: Does ITopologicalOperator.Simplify call IGeometry.SnapToSpatialReference?

texts:
Does the ITopologicalOperator.Simplify method call IGeometry.SnapToSpatialReference as part of the simplify process?

I thought I saw some ESRI documentation that indicated that it does, but I cannot find it now.   I would prefer not to make a redundant call if it is not needed.  If you have a link to ESRI documentation indicating that it does, please include it.


texts:
I have a map of streets and buildings that I need to change their angel on the map. It is also necessary that the name labels will change respectively.
This is a screen shot of the map


texts: create diversity map of qualitative features, ArcMap

texts:
I am trying to create a map of forest type diversity with ArcMap. I have a polygon forest layer with the different forest types as attributes (e.g. 0815, 0862, 0813, ..., which stands for 'mixed forest', 'oak forest', and so on.)

Using Patch Analyst I have created a Hexagon map and intersected it with the forest polygon layer. Usually I would now use 'field statistics' to calculate the variance of forest types for each Hexagon field and then join the result table to the Hexagon polygon layer to create a diversity map. Calculations don't work, however, since I have qualitative Features instead of quantitative ones.

Does somebody have an idea how to do this?


texts: Delete duplicate text in Atribute table, qgis

texts:
I'm trying to label correctly a set of points in my map. They refer to localities and water analyses for those places in Spain. But the information in one of the fields is duplicated, it's written two times. I'm using QGIS 2.18

I'm trying to delete the duplicate text from a field in my Attribute table, but so far, I've not been able to find the right expression. The Field is called Municipio.
I've tried with the expression:
replace ( "Municipio", 'Municipio', ' ' )
but it's not working.

Any idea ?


texts: Is there a way to have a location bar (or breadcrumb trail) to show the location path in QGIS 3 browser panel? Most importantly in order to copy-paste a location and go to it directly without having to manually navigate the tree ...


texts: How to zoom to full map extent with mapserver?

texts:
I need to put checkbox with name "Zoom to extents" which will zoom to full map extent (first extent which appears by opening web page). Can someone help me?
In mapfile I have this "EXTENT 563090 5072655 588220 5093325".


texts: subtype & domain from gdb to table (arcgis)

texts:
We have a geodatabase containing an obj class (main table) which has subtypes & these contain domains.

Need a way (e.g. script) to get those subtypes & domains in a "readable" table.

Any help?


texts: Extract geometry from citygml data using citygml4j

texts:
I have a citygml file with some buildings. And I am trying to extract the geometry of the polygon in that building using the citygml4j library. I have looked at the FeatureWalker class but I am unable to get the geometry of the polygon.

Can anyone let me know how this can be done?

    CityGMLContext ctx = new CityGMLContext();
    CityGMLBuilder builder = ctx.createCityGMLBuilder();

    CityGMLInputFactory in = builder.createCityGMLInputFactory();
    CityGMLReader reader = in.createCityGMLReader(new File("/home/vishal/NWW/sampleData/LOD2_Building_v100.gml"));

    while(reader.hasNext())
    {
        CityGML citygml = reader.nextFeature();
        System.out.println("Found class:" + citygml.getCityGMLClass() + "\nVersion"+citygml.getCityGMLModule().getVersion());

        //Counting the no of buildings
        CityModel citymodel = new CityModel();
        if(citygml.getCityGMLClass() == CityGMLClass.CITY_MODEL)
        {
            citymodel = (CityModel)citygml;
            // Counting the no of buildings
            int count=0;
            for(CityObjectMember cityObjectMember : citymodel.getCityObjectMember())
            {
                AbstractCityObject cityobject = cityObjectMember.getCityObject();
                if(cityobject.getCityGMLClass() == CityGMLClass.BUILDING)
                {
                    ++count;
                }
            }
            System.out.println("Building count"+count);
        }

        FeatureWalker walker = new FeatureWalker(){
            public void visit(Building building){
                System.out.println(building.getId());
                //MultiSurface multisurface = boundrysurface.getLod2MultiSurface().getMultiSurface();
                //System.out.println(multisurface.getSurfaceMember().get(0));
                List<BoundarySurfaceProperty> list = building.getBoundedBySurface();
                System.out.println(list);
                System.out.println(list.get(0).getBoundarySurface());
                //HOW TO GET THE POLYGON AND ITS COORDINATES??
            }
        };
        citymodel.accept(walker);


texts: Using AutoCAD to GIS Standards for Easy Import?

texts:
I am trying to setup an AutoCAD Template that will be easy to import into an ArcGIS geodatabase.

Is there a way I can setup a template where the data tables are already attached so a Tech can just select a PVC 150mm BluBrute Pipe and draw it in, then maybe select a 19mm Copper Service Lateral and draw that in?

The only other way I can see doing this is to make each pipe of any given combination of attributes it's own layer. I don't like this because the layer schema to keep track of it starts getting complicated and convoluted.


texts: Uisng Python Parser of Field Calculator with ArcPy?

texts:
I have a point shapefile with two columns i want to process.
The first column is aspect and the second is aspect_m60, there is only one row with the aspect of my point. I have tried to write a script in Python for ArcGIS 10.2.2 because i want to put it in my code but with no success.
The code is this:

# Calculate Field
import arcpy
# Set environment settings
arcpy.env.workspace = "c:/W/Sik"
# Set local variables
inTable = "Point"
fieldName = "aspect_m60"
expression = "getCalc(!aspect!, !aspect_m60!)"
codeblock = """def getCalc(aspect, aspect_m60):
   if (aspect < 60):
      aspect_m60 = (aspect - 60) + 360
   if (aspect = 60):
      aspect_m60 = aspect - 60
   else:
      aspect_m60 = aspect - 60"""


texts: Converting CAD Annotation to geodatabase feature class using ArcPy?

texts:
Is there a way to list through DWG and pick out let's say only the annotations and feed them into a geoprocessing tool?

Similar to below in a feature class...

for ds in datasets:
    #Finds polyline feature classes and appends to sde
    for poly in arcpy.ListFeatureClasses("Polyl*",feature_dataset=ds):
        arcpy.Append_management(poly,  .....

I would like to pick out the annotations and feed them into the Import annotation tool. I can't seem to get it. The link below makes it look so easy but it doesn't work.

http://desktop.arcgis.com/en/arcmap/latest/manage-data/cad/using-python-to-load-cad-data.htm


texts: Additional No Data Value Raster Pyqgis 3

texts:
I would like to improve way to adjust histogram of Sentinel 2 images in Qgis 3. I am using pyqgis to do the first task to it (create histogram rightly). So I have images that 0 is a no data value. To set it I use this code:
Set raster layer additional no data value with PyQgis

See below:

import processing
self=qgis.utils
layers = self.iface.mapCanvas().layers()

path_save = '/home/diogocaribe/sld/'

# Change o backgroud value
for layer in layers:
    layerType = layer.type()
    if layerType == QgsMapLayer.RasterLayer:
        for band in range(1, layer.bandCount()+1):
            provider = layer.dataProvider()
            provider.setNoDataValue(band, 0) #first one is referred to band number
            layer.triggerRepaint()

So, the first thing that I have to do is create histogram. However, if I don't set 0 in additional no data value the histogram is created with different behavior:

When I defined additional no data value as 0 and recompute histogram it works well.

So I would like to fill the blank in proprieties layer using pyqgis in order to put additional no data field as 0 and fix the problem described above.


texts: Importing .shp file into PostGIS with the function shp2pgsql

texts:
I am trying to import a .shp file into PostgreSQL, but so far no success.
After completing the setup, I installed the PostGIS extension and finally went for some geo exploration. But, when I try to import a .shp file, this error appears:
root@e61fec68cbc2:/# shp2pgsql -s SRID br_shp_dir/BR_UF_2019.shp public.geo_br | psql -h localhost -d db-example -U postgres;
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
SET
SET
BEGIN
ERROR:  relation "geo_br" already exists
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ERROR:  current transaction is aborted, commands ignored until end of transaction block
ROLLBACK
ANALYZE

root@e61fec68cbc2:/#

I working Docker, so It's containerized.
I've done it before, but couldn't this time.
What is wrong?


texts: WFS 2.0 with bbox filter and attribute filter

texts:
I am trying to get the following WFS GET request to work:

wfs = "service=WFS&version=2.0.0&request=GetFeature",
name = "&typeName=storeName:layerName",
bboxFilter = '&bbox=' +  bounds._southWest.lng + ',' + bounds._southWest.lat + ',' + bounds._northEast.lng + ',' + bounds._northEast.lat,
attributeFilter = '&name=test',
format = "&outputFormat=text/javascript&format_options:callback:getJson",
urlWFS = "http://domain.com:8080/geoserver/bla/ows?" + wfs + name + bboxFilter + attributeFilter +  format;

It works fine with the bounding box filter, but the attribute filter is simply ignored. What is the right syntax for that?


texts:
I'm trying to export a shapefile to a csv file using a cursor. I found this quesion asked and answered here: Exporting table to XYZ ASCII file via ArcPy?

However I only want to export several specific columns in a specific order. I have been able to get the .csv header to export correctly but can't get the cursor to provide the data for the rows in that same order.

import arcpy
# import csv module
import csv

# Set working file location
dataBase = r"C:\Users\vince.adams\Desktop\VLA_FinalProject\GEARY"

# Set output file location
outFile = r"C:\Users\vince.adams\Desktop\VLA_FinalProject\Tables"
tableEight = outFile + "\Table8.csv"

# S_GEN_STRUCT
genStruct = dataBase + "\S_GEN_STRUCT.shp"
fieldList = ["WTR_NM", "STRUCT_NM", "STRUCT_TYP", "LOC_DESC", "STRUC_DESC"]

#--first lets make a list of all of the fields in the table
fields = arcpy.ListFields(genStruct, fieldList)
field_names = [field.name for field in fields]

with open(tableEight,'wb') as outTable:
    writeOut = csv.writer(outTable)
    #--write all field names to the output file
    writeOut.writerow(fieldList)

    #--now we make the search cursor that will iterate through the rows of the table
    for row in arcpy.SearchCursor(genStruct, fieldList):
        field_vals = [row.getValue(field.name) for field in fields]
        writeOut.writerow(field_vals)

    del row

I'm using Arc 10.2.


texts: How to read OpenLayers event object features?

texts:
This is part of an interface that uses a box to select features on a WMS Layer. This then writes the name of the features to a div as a list.

var recycling = new OpenLayers.Layer.WMS(
    "Recycling Sites",
    wmsurl+'?',
    {layers: 'GeoserverWorkspace:recycling', 'format':'image/png', 'transparent':'true'},
    {'opacity': 1.0, 'isBaseLayer': false, 'visibility': true}
    );

 select = new OpenLayers.Layer.Vector("Selection", {styleMap:
            new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])
        });

map.addLayers([osmap, recycling, select])

var recycling_wfs =  new OpenLayers.Control.GetFeature({
            protocol: OpenLayers.Protocol.WFS({
            url: wfsurl,
            version: "1.1.0",
            featureType: "recycling",
            featurePrefix: "GeoserverWorkspace",
            srsName: "EPSG:27700",
            geometryName: "wkb_geometry",
            maxFeatures: 50,
            click: true
            }),
            box: true
});

map.addControl(recycling_wfs);

recycling_wfs.events.register("featuresselected",this, showInfo);

recycling_wfs.activate();

function showInfo(e) {

    // assign features object to a variable

    var features = e.object.features;

    //properties are accessed by using the attibutes property
    var selectedSites = "";

        for(var fid in features){
          var site = features[fid].attributes;
          var siteName = site.name;
          selectedSites += siteName + " ";
        }

    document.getElementById('responseText').innerHTML = "Recycling Site: " + selectedSites;
}

I am confused about the event object that is returned by the featuresselected event, to the callback function showInfo(). According to the OpenLayers documentation this should be an array. So if I select 2 features I should be able to access the first name by doing the following

var name = e.object.features[0].attributes.name;

However this returns an error of e.object.features[0] undefined. I am able to get the feature name if I use a For..In and iterate through the event.object.features object. Then use the resulting fid to access the feature, like so

var name = e.object.features[fid].attributes.name

Is this the correct way to access the event.object's features or have I missed somthing?


texts: QGIS 2.81. crashes when using OpenLayers

texts:
I just installed QGIS 2.8 Wien to get access to OpenLayers (I could not install it under 2.4, and had to uninstall completely and reinstall 2.8).
It was working apparently, but when I tried to get to the OpenLayers

OpenLayers Plugin not working in QGIS 2.8.1

The last comment in that question is that the plugin has something to do with a ./qgis2/python/plugins/openlayers_plugin/weblayers/htfolder, was of no help.

That question links to this one:

Where is the OpenLayers plugin in QGIS 2.6?

I checked out if the OpenLayers plugin was installed and activated in the "Installed Plugins" list. It is, but it doesn't show the OpenLayers  under the menu "Web" yet. I try to use "OpenLayers" under the "Plugins" menu, selecting "Add Google Satellite Layer" but QGIS crashes completely (closes).

I could not comment or answer in the mentioned questions, so I started a new question because the problem is still not solved.


texts:
In Mapbox, this example works on a PC but not on a touchscreen Android/iPhone/iPad

https://www.mapbox.com/mapbox.js/example/v1.0.0/opacity/

Is there a simple change to the code below to make it work on mobile devices?

==================

<!DOCTYPE html>

<meta charset=utf-8 />
<title>Opacity control</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

<link href='https://api.mapbox.com/mapbox.js/v2.4.0/mapbox.css' rel='stylesheet' />


texts: How do I make my raster layer transparent where there is no data?

texts:
I have just managed to load some TIF data into PostgreSQL 11 using raster2pgsql. This test contained just a few files.

When I load the resulting layer into QGIS, the data displays perfectly, expect all blank "no data" areas which are just black.

How can I stop this happening and make these areas transparent? Can I set this transparency back at database level, or can it only be manipulated within the end client (i.e.QGIS)? The ideal solution would be to sort this at database level, saving the hassle of changing it every time it's loaded into a new GIS session.


texts: The results of xy coordinate transform in R and arcmap are different

texts:
Some lon lat sites were converted to xy geographic coordinates in R (target geographic coordinate system EPSG:102012) :

     lon   lat
1 120.11 50.15
2 121.41 50.29
3 117.19 49.35
4 119.42 49.15
5 119.56 47.10

lonlat <- data.frame(lon=c(120.11,121.41,117.19,119.42,119.56),
                     lat=c(50.15,50.29,49.35,49.15,47.1))

inproj <- '+proj=longlat +datum=WGS84 +no_defs'
outproj <-  '+proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105
+x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0'

coordinates(lonlat) <- c(1,2)
proj4string(lonlat) <- inproj
xy <- as.data.frame(spTransform(lonlat,outproj))
xy
        lon     lat
1 1033208.0 5981950
2 1117747.6 6014440
3  848573.9 5862626
4 1006104.0 5867861
5 1056140.0 5654292

I added these lon lat points (WGS84 projection) to arcmap as shape points and changed to geographic projection. Then using Calculate Geometry in attribute table.The xy coordinate got as follows:

1033208.041 5981950.466
1117747.601 6014439.722
848573.8694 5862626.143
1006103.956 5867860.72
1062292.81  5684395.771

inproj Lambert_Conformal_Conic(LCC) ESRI WKT projection in arcmap:

PROJCS["LCC_WGS_1984",GEOGCS
    ["GCS_WGS_1984",DATUM
       ["D_WGS_1984",SPHEROID
       ["WGS_1984",6378137.0,298.257223563]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
       PROJECTION["Lambert_Conformal_Conic"],
       PARAMETER["false_easting",0.0],
       PARAMETER["false_northing",0.0],
       PARAMETER["central_meridian",105.0],
       PARAMETER["standard_parallel_1",30.0],
       PARAMETER["standard_parallel_2",62.0],
       PARAMETER["latitude_of_origin",0.0],
       UNIT["Meter",1.0]]

The results of first four sites are same between two methods, but there is bias in fifth point. That puzzle me.

Origin of this problem: I want to extract raster values to points.Fake codes as follows:

xy <- spTransform(lonlat,crs(raster))
values <- raster::extract(raster,xy)

I found some points values were NA. So I checked the point spatial distribution that overlap shape points and raster in arcmap. All points are in raster. but some the coordiates tranformed in R were out of raster. That is why the extracted values are NA.


texts:
hoping someone can help me with some styling on a map I've created. I'm rather inexperienced here, but feel rather close to success so hopefully this is a layup for someone who knows what they're doing.

Im basically trying to style my data exactly as I want it, which is uneven intervals. Is there a way to input what I'd like the intervals to be, as I'm able to when I create my legend, and have custom intervals color coded on the map? Currently, I have US counties and customer zipcodes populating the map for me, with a range of 1-1700 per county. I've got 7 buckets and equal intervals on my map, but not sure if quantiles, jenks, or categories would be better, not too much changes as I select them.


texts: How to target multiple specific layers when using -merge-layers

texts:
As part of trying to find an answer for How to -explode only one subfeature from only one feature? I created some fields in one layer only, then, after doing some operations, tried to merge that layer back with one other specific layer.

The docs for merge-layers say you can use target= and name= but I've not managed to get either to work:

Apparently target= can accept a comma-separated list of layers, but it seems not to work:

If I wrap each layer name in quotes:

-merge-layers target="'layer A', 'layer B'"

...I get an error suggesting it's looking for one layer called 'layer A', 'layer B':

  Missing target: 'layer A', 'layer B'

  Available layers: [list of layers here]

If I don't wrap each layer name in quotes:

-merge-layers target="layer A, layer B"

... then I get no errors, but nothing happens. It doesn't merge the layers.
If I add a name then it renames only the first layer in the comma separated list (regardless of what layer is selected). For example, this:

-merge-layers target="layer A, layer B" name="new layer"

...renames layer A to "new layer" but does nothing at all to layer B. So for whatever reason it seems to ignore the layers after the first in the comma seperated list.

I also tried using name like this:

-merge-layers name="'layer A', 'layer B'"

...which ran into this error, despite the fact I had used -filter-fields to remove these fields from the layers I was targeting:

  Field missing from one or more layers: someField

...but I later discovered that this was because I'd misunderstood the name parameter: this was simply the name the new merged field would take, so since I had no target defined, it was trying to merge all layers and was picking up on errant fields in layers other than the ones I was trying to merge.
I also tried using one layer name as a regular parameter and one as a target= but it didn't like that either:

-merge-layers 'layer A' target="layer B"

...gives this error:

  Unexpected option: layer A

I also tried two target options, but this doesn't work:

-merge-layers target="layer A" target="layer B"

All the usage examples I can find for -merge-layers use it to merge every layer.

What's the correct way to target and merge specific layers?


texts: My MapViewer stopped loading the layers unexpectedly and giving error 500

texts:
I have developed a  MapViewer with around 140 different layers with renewable energy information from the West African region which are all public and free for download. www.ecowrex.org/mapView The layers (vector files and raster files) are stored in GeoServer here:
www.ecowrex.org:8080/geoserver/web/?wicket:bookmarkablePage=:org.geoserver.web.demo.MapPreviewPage

Some days ago out of nowhere all of the layers have failed to load when selected. When I inspect the console for errors in the browser I see the error "500 (Internal Server Error)" each time it attempts to fetch the layer.

For some reason in the console when it tries to access the url of the layer, the port :8080 is not there after the host. But if I copy and paste the url to my browser and add the port myself I see the layer (ex: http://www.ecowrex.org:8080/geoserver/ecreee/wms?TILED=false&LAYERS=ecreee%3Aecowas_country_admin&STYLES=&FORMAT=image%2Fpng&SRS=EPSG%3A900913&TRANSPARENT=TRUE&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&BBOX=-7239270.7789323,-100001.97671781,6125590.7408137,3451568.1050302&WIDTH=1366&HEIGHT=363).

Can that be an indication of what is going wrong?


texts:
I have this WFS: http://sg.geodatenzentrum.de/wfs_vg1000?SERVICE=WFS&VERSION=2.0.0&REQUEST=GetFeature&TYPENAMES=vg1000:vg1000_sta&MAXFEATURES=5&BBOX=47.567208267268484,10.44518218531835,48.949175304644434,13.626305641159272,urn:ogc:def:crs:EPSG::4258&SRSNAME=urn:ogc:def:crs:EPSG:4258

As you can see, the version requested is 2.0.0 so the key for limiting the number of feature types as result should be COUNT and not MAXFEATURES which is for the earlier WFS versions. COUNT does not work but MAXFEATURES is working fine. Very strange...

when clicking on the same link fot the request, BBOX is actually not returning the elements in the BBOX but all the elements... the BBOX is not working properly.

Is there something that I am missing here with this service, or is the service in a bad shape?


texts:
My QGIS Installation is not loading after startup, cause of an error thrown. I have looked in similiar topics like QGIS fails to load and tried all of the options available. I changed to scipy 0.17 and 0.11, removed qgis entirely, removed qlibhull7, installed all again, removed .qgis2, tried to start with --noplugins. Nothing worked so far.

my default installation:

qgis 2.18
ubuntu 16.04
scipy 1.0.0

it's always the same error

Warning: loading of qgis translation failed [/usr/share/qgis/i18n//qgis_en_US]
Warning: loading of qt translation failed [/usr/share/qt4/translations/qt_en_US]
Warning: QCss::Parser - Failed to load file  "/style.qss"
QInotifyFileSystemWatcherEngine::addPaths: inotify_add_watch failed: No such file or directory
Warning: QFileSystemWatcher: failed to add paths: /home/cheeseman/.qgis2//project_templates
Warning: QLayout: Attempting to add QLayout "" to QgsPanelWidgetStack "mWidgetStack", which already has a layout
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
Warning: QVariantMap DBusMenuExporterDBus::getProperties(int, const QStringList&) const: Condition failed: action
QH6248 qh_lib_check: Incorrect qhull library called.  Caller uses reentrant Qhull while library is non-reentrant
QH6249 qh_lib_check: Incorrect qhull library called.  Size of qhT for caller is 8184, but for library is 2896.
QH6255 qh_lib_check: Cannot continue.  Library 'qhull 7.2.0 (2015.2 2016/01/18)' uses a dynamic qhT via qh_QHpointer (e.g., qhull_p.so)


texts: Manage representions by updating rules with original data

texts:
In case your not familiar with representations.
There is new since 9.2 a cartographic representation available.
It allows for rules based on data fields to determine symbology.

Admittedly it is very useful and powerful for single entity manipulation to show it how and where you want.
One problem with it is that there is only one way to get a representation from existing symbology.

Using "convert symbology to representation" it creates a new feature class (db only) and creates all the rules and symbology for the representation.  These symbols are not the normal esri symbols so there are some issues and restrictions for maintining that.
However at that point you have a (personal, file or sde) gdb with symbology in it.
Whenever you add that data to a new map the symbology is already there and if you have more than one rep you simply click between the reps you have defined.

As I said very powerful.
The REALLY big drawback I just found is that when you edit the original value that created the rule the rule is not updated.

An easy enough fix if you are symbolizing with only one field.
As luck would have we are not.
Now the question.

Has anyone used a script to update the rule_id field from more than one originating value?


texts: How to extract specific values with point coordinates from Kriging interpolations made in R?

texts:
By using R version 3.4.2 and the library "geoR", I made kriging interpolations for different variables (bellow I give an example of my process). I also made a matrix with the coordinates for 305 trees with distinct marks (species, DBH, Height) that are within the same space for the interpolations, as seen in the image attached (https://imgur.com/SLQBnZH). I've been looking for ways to extract the nearest value from each variable for each tree and save the corresponding values in a data.frame or matrix, but haven't been successful, and I can't find specific answers to this.

One thing I've been looking at is trying to convert the Kriging result into a Raster (.tif) and proceed from there. But Kriging interpolations are made out of vector data, so is it even posible?

I'm doing this so that I can latter use the data for spatial point patern analysis.

#Kriging####:
PG<-read.csv("PGF.csv", header=T, stringsAsFactors=FALSE)
library("geoR")
x<-(PG$x)
y<-(PG$y)

#Grid
loci<-expand.grid(x=seq(-5, 65, length=100), y=seq(-5, 85, length=100))
names(loci)<-c("x", "y")

mix<-cbind(rep(1,10000), loci$x, loci$y, loci$x*loci$y)

#Model
pH1.mod<-lm(pH1~y*x, data=PG, x=T)
pH1.kg<-cbind(pH1.mod$x[,3], pH1.mod$x[,2], pH1.mod$residuals)
#Transform to geographic data
pH1.geo<-as.geodata(pH1.kg)
#Variogram
pH1.vario<-variog(pH1.geo, max.dist=35)
pH1.vario.mod<-eyefit(pH1.vario)
#Cross validation
pH1.valcruz<-xvalid(pH1.geo, model=pH1.vario.mod)
#Kriging
pH1.krig<-krige.conv(pH1.geo, loc=loci, krige=krige.control(obj.model=pH1.vario.mod[[1]]))
#Predictive model
pH1a.yhat<-mix %*% pH1.mod$coefficients + pH1.krig$predict
#Exchange Kriging prediction values
pH1.krig$predict<-pH1.yhat
#Image
image(pH1.krig2)
contour(pH1.krig2, add=TRUE)

#Tree matrix####:

CoA<-read.csv("CoAr.csv", header=T)
#Data
xa<-(CoA$X)
ya<-(CoA$Y)
points(xa,ya, col=4)

TreeDF<-(cbind.data.frame(xa, ya, CoA$Species, CoA$DBH, CoA$Height, stringsAsFactors = TRUE))
m<-(cbind(xa, ya, 1:305))
as.matrix(m)

I tried to find the value of a point in space (trees [1:305]) through the minimum distance to a predicted value using the following code, (I suggest not running this since it takes too long):

for(i in 1:2){print(c(2:10000)[as.matrix(dist(rbind(m[i,], as.matrix(pH1.krig2$predict))))[i,2:10000]==min(as.matrix(dist(rbind(m[i,],as.matrix(pH1.krig2$predict))))[i,2:10000])])}


texts: SplineWithBarrier (Spatial Analyst) - different results between model builder and external python script

texts:
I'm puzzled by different results I get from application of SplineWithBarrier (Spatial Analyst) using Modelbuilder and by external Python scripting:

  import arcpy
  points_shp = "pointShape.shp"
  faults_shp = "faultsShape.shp"
  outGrid = "r_out"
  arcpy.gp.SplineWithBarriers_sa(points_shp, "id", faults_shp, "10", outGrid, "1")

This is the code I got from the model builder. While results there are as expected (smooth surface with discontinuities) I used it in an external script. The resulting surface derived from the external script, however, shows, that smoothing was skipped. Before I tried the version from model builder I followed the example from the ArcGIS help. Same result - the grid is created, but again smoothing step is skipped. My first code was:

  import arcpy
  points_shp = "pointShape.shp"
  faults_shp = "faultsShape.shp"
  outGrid = "r_out"
  arcpy.CheckOutExtension("Spatial")
  outSplineBarriers = SplineWithBarriers(points_shp, "id", faults_shp, 10, 1)
  outSplineBarriers.save(outGrid)

That's the syntax from the help:

SplineWithBarriers (Input_point_features, Z_value_field, {Input_barrier_features}, {Output_cell_size}, {Smoothing_Factor})

Why is smoothing not activated? Where is my mistake?


texts: We are using Mapbox Mobile SDK for Android application Development. We have used our own Maps Data SETs for POIs and Addresses.  These custom POIs are not available to us in Search and Creating Routes.  Is there a possibility of using Custom Maps with our own POIs and implement search through Cartodb?


texts: Create image from ESRI VectorTileServer using script

texts:
I have service from ESRI VectorTileServer

https://tiles.arcgis.com/tiles/P3ePLMYs2RVChkJx/arcgis/rest/services/Esri_Childrens_Map/VectorTileServer/tile/{z}/{x}/{y}.pbf

Map example:

https://tiles.arcgis.com/tiles/P3ePLMYs2RVChkJx/arcgis/rest/services/Esri_Childrens_Map/VectorTileServer?f=jsapi

There is possible create print image (jpg or png) for specified bbox coordinates and zoom level (scale)? What tools can I use for creating this images?

For creating images from ESRI VectorTileServer I finded this example:

https://developers.arcgis.com/javascript/latest/sample-code/sandbox/index.html?sample=widgets-print-vectortiles

but there is not available call print using for example python script.


texts: Move tool feature in digitizing toolbar doesn't work/ not available

texts: I can't move any polygons and the layer is in the edit mode, but the move feature button in the digitizing toolbar doesn't appear to be enabled to click it.


texts: How to make a raster layer transparent via MapServer WMS

texts:
I loaded a PostGIS raster data through OpenLayers and MapServer.

The code inside the OpenLayers is as following :

localWMS = new OpenLayers.Layer.WMS(
    "Local WMS layers",
    "http://127.0.0.1:8080/cgi-bin/mapserv.exe?MAP=C:/MapServer/raster.map",
    {
        layers: 'polygon',
        transparent: true,
        alpha:true
    },
    {
        projection:new OpenLayers.Projection("EPSG:4326"),
        units:'m',
        isBaseLayer: false
    }
);

map.addLayer(localWMS);

And here is my map file, I already set opacity to each class of the layer but nothing happens:

LAYER
    DATA "PG:host=*** dbname=*** schema=public table=london where='rid=3' mode='1'"
    NAME "postgis"
    PROCESSING "SCALE=AUTO"
    PROCESSING "NODATA=-9999"
    STATUS DEFAULT
    TYPE RASTER
    UNITS METERS
    CLASS
      NAME "Minor Impact"
      EXPRESSION ([pixel] >= 0 AND [pixel] < 1.29525335532 )
      STYLE
        ANGLE 360
        OPACITY 50
        COLOR 122 182 245
        SYMBOL 0
        WIDTH 1
      END
    END
    CLASS
      NAME "Moderate Impact"
      EXPRESSION ([pixel] >= 1.29525335532 AND [pixel] < 2.59050671063 )
      STYLE
        ANGLE 360
        OPACITY 50
        COLOR 76 230 0
        SYMBOL 0
        WIDTH 1
      END
    END
    CLASS
      NAME "Severe Impact"
      EXPRESSION ([pixel] >= 2.59050671063 AND [pixel] < 5.69911476339 )
      STYLE
        ANGLE 360
        OPACITY 50
        COLOR 252 228 91
        SYMBOL 0
        WIDTH 1
      END
    END
    CLASS
      NAME "Catastrophic"
      EXPRESSION ([pixel] >=5.69911476339 )
      STYLE
        OPACITY 30
        COLOR 232 16 20
      END
    END
END

The map looks like following :

Then, I tried to set one of the opacity to less than 1, the red color just disappeared.

But if I set the Opacity of the Layer to 50

LAYER
    DATA "PG:host=*******************'"
    NAME "postgis"
    PROCESSING "SCALE=AUTO"
    PROCESSING "NODATA=-9999"
    STATUS DEFAULT
    TYPE RASTER
    OPACITY 50
    UNITS METERS

Even more strange map appeared.

Does anyone know what's the problem about that? Or Does anyone know the proper way to set a transparent raster layer via MapServer?


texts:
I am trying to return the results of queryByShape call from a (large) raster file given a user specified polygon.

For example, when a user plots the following polygon on an openlayers layer as shown in this diagram, I send the coordinates of the polygon vertices to the python mapscript backend.

With python mapscript, I re-create the polygon as a mapscript shapeObj.

line = mapscript.lineObj()
line.add( mapscript.pointObj( 16.484375, 59.736328125 ) )
line.add( mapscript.pointObj( 17.4951171875, 55.1220703125 ) )
line.add( mapscript.pointObj( 24.74609375, 55.0341796875) )
line.add( mapscript.pointObj( 22.5927734375, 61.142578125) )
line.add( mapscript.pointObj( 16.484375, 59.736328125) )
poly = mapscript.shapeObj( mapscript.MS_SHAPE_POLYGON )
poly.add( line )

Then, I experiment with the queryByShape call like this:

mapfile = '/path/to/my/mapfile.map'
mymap = mapscript.mapObj( mapfile )

layer = mymap.getLayer(0)
layer.queryByShape( mymap, poly )

The returned result is an integer "1", which isn't what I was expecting.

Unfortunately, there isn't much by way of documentation for queryByShape so I am not sure what I am doing wrongly (http://mapserver.org/mapscript/querying.html#by-shape)

How does mapscript queryByShape work? I am actually trying to return a shape object where I can extract information about the selected shape (polygon defined by user) and determine

The total area of this shape,
The % in this shape which is of a specific RGB ("blue") value,
The % in this shape which is of a specific RGB ("green") value etc.

Here's the LAYER definition in my mapfile.

LAYER
    NAME 'GLOBCOVER_L4_200901_200912_V2.3.color'
    TYPE RASTER
    DUMP true
    TEMPLATE 'fooOnlyForWMSGetFeatureInfo'
    EXTENT -189.001389 -115.106069 188.998611 140.108847
    DATA '/media/GLOBCOVER_L4_200901_200912_V2.3.color.tif'

    METADATA
        'ows_title' 'GLOBCOVER_L4_200901_200912_V2.3.color'
    END

    STATUS OFF
    TRANSPARENCY 100

    PROJECTION
        'proj=longlat'
        'ellps=WGS84'
        'towgs84=0,0,0,0,0,0,0'
        'no_defs'
    END
END


texts: eCognition - prevent rounding of raster values?

texts:
I have a WorldView2 image, which I would like to segment and classify in eCognition. The image has been pre-processed to Top of the Atmosphere Reflectance, therefore hasvalues between 0-1.

The segmentation process in eCognition does not show any meaningful results and I think that is likely due to the fact that it does round decimal values. Is there a way to prevent that from happening? I would rather not have to convert to int.


texts: Reorder the features in a layer to display on infowindow

texts:
I am trying to display features from results on infowindow, I have written the following code for that, But I am not getting the output.

Please suggest me where to modify in my code to get first type "CIT" and next type "CNY" on infowindow.

promises.then(function (results) {
                    //results are in same order as were passed to the all constructor
                    var showFeatures = [];
                    if (results.length > 0) {
                        dojo.forEach(visibleLayers, function (layer, index) {
                            //set the info template for each returned feature
                            dojo.forEach(results[index], function (result, fIndex) {

                                if (results[index].length>0) {
                                        if(result.feature.layerId == 2){
                                            if (result.feature.attributes.munitype == "CNY"){
                                                result.feature = results[4].reverse();
                                            }
                                            else {
                                                result.feature = results[4];
                                            }

                                        }
                                        /*else{
                                        result = results[4];
                                        }*/

                                        else {
                                            result.feature.infoTemplate = layer.infoTemplate;
                                            showFeatures.push(result.feature);
                                        }

                                }

                            });
                        });

                        if (showFeatures.length > 0) {
                            app.mapUtil.map.infoWindow.setFeatures(showFeatures);
                            app.mapUtil.map.infoWindow.show(mapPoint);
                        }

                    }
                }, function (error) {
                    app.logger.logError("showInfoPopups Error searching for results", error);
                })


texts:
I am running this script . I get a 400 error message. The message is the following:

{"error":{"code":400,"message":"Unable to complete operation.","details":["Server environment property not found: requestProperties"]}}

I am not sure what it is. I am using the Extract Data Task (Server) tool for 10.2.


texts: How can I change more than one attibute at the same time?

texts:
I'm working in QGIS 2.0.1 Dufour. I started by creating a shapefile from an AutoCAD file. The problem is that in the AutoCAD, my layers had symbols like: "~; ç; ^". As the QGIS does not support thess symbols, my attribute table is really confusing, so I want to rewrite those attributes.

What I mean is, in one column, I want to change more than one line.
My problem is that I have more than 100 of each I want to rewrite. I thought it could be a way of changing more than one at the same time.

Is it possible?


texts: Choosing Geographic Transformation from D_WGS_1984 to D_Hughes_1980 in ArcMap?

texts:
I am using ArcMap 10.2 and I want to convert the shapefile projection(from GCS_WGS_1984 to NSIDC_Sea_Ice_Polar_Stereographic_North), however, their datum is different,the datum of GCS_WGS_1984 is D_WGS_1984 and NSIDC_Sea_Ice_Polar_Stereographic_North is D_Hughes_1980. Therefore, Geographic Transformation is must be filled（as is shown in Fig.1), however, I cannot find D_WGS_1984_To_D_Hughes_1980 in drop-down list of valid transformation methods, as Fig.2.

I want to know the content to fill in Geographic Transformation field.

                     Fig.1

                     Fig.2


texts: Changing projection from 0 360 to -180 180 for a simple raster

texts:
I'm trying to convert a raster with projections of (0, 360) to (-180, 180), and I've tried various SO solutions such as this, but it gave me a cropped output that doesn't have the western hemisphere.

The exact command that I've used is:

gdalwarp -t_srs EPSG:4326 -te -180 -90 180 90 prec_high_res.tif rec_high_res2.tif

Does anyone have any idea how to resolve this?


texts: How to create ROI in Semi-Automatic Classification plugin QGIS

texts:
while trying to classify landsat images, I get the following problem:

I click on the pixel I want to define a ROI, but the ROI polygon that the plugin creates appears in a different part of the image instead of appearing over the selected pixel...
The pink arrow shows the pixel I clicked to create the ROI and the blue arrow shows the ROI polygon that appears. Does anyone know how to fix this?

It was a projection problem: I projected the images to it's original projection and after reclassifying the image I changed the projection.


texts: Practical Building Extraction Algorithm Implementation

texts:
I've been reading many papers over 3D building extraction using LiDAR data and Aerial images. I see that there are many different algorithms for doing such things, such as this and this.

Those papers, as far as I can see, describe mostly abstract concepts and math. I am a software developer, so how could I practically implement said algorithms into a piece of software? Also, what algorithm is best for extracting building data with minimal vegetation interference?


texts: Using Layer Parameter variables in Field Calculator within QGIS Model?

texts:
If I define a Vector Layer as an input parameter to a QGIS model, then the parameter name is exposed as a variable in a Field Calculator within the model, however I am unable to use this variable as I am not sure what data type it is. For example with a Vector Layer parameter of InData, I have tried evaluating the following expressions into a text field in a temporary layer:

@InData => ''

layer_properties(@InData, 'name') => NULL

layer_properties(@InData, 'id') => NULL

parameter( 'InData') => NULL

This is the model:

Variables exposed in Field Calculator (the expression doesn't seem to be case sensitive):

My question is - how do you use the @indata variable? For example, can I use it to extract the CRS of the input layer?

A workaround is to use one of the @layer / @layer_id variables in a field calculator as the first tool connected to the input in the model, but then the result will need to be stored as an attribute if it is needed later in the model (see How to access input layer properties from QGIS' Modeler?).


texts:
I am trying to write a script that will go through a grid layer that has individual 535 grid squares.

I want to find a way that will tell me how much gas main is left to be checked in each grid. Any pipe newer than 2011 does not need to be checked.

I would like to write a script that could be run on a weekly basis to check the progress of the gas main layer because it is being edited by numerous users. To manually do this for grid square I usually do select by attributes with the query "Grid number = (#)", the layer has from 1 to 535, after I use select by location to select all gas main that is inside that individual grid square, then I remove the pipe that has prime and is newer than 2011, then export it to an Excel file.

I want the output of the script to create an Excel file with all remaining OIDs in the gas main layer to be validated for each individual grid square, so 535 different Excel sheets.

How do I get the value that is in the grid number field and store it as a variable to use in the select by attributes tool ("Grid_Number = (variable)" and change as it iterates through each row?

import arcpy
count = arcpy.da.UpdateCursor("keymap", "miles_left")

with count as cursor:
    if(row0]>=0 and row[0] <= 354):
        select by attributes(key map = %)
        select by location (gas main in kepy map)
        select by attributes (remove prime, remove new than 2011)
        export to excel
        update miles_left field


texts: Cross layer Filtering in Geoserver still getting error Could not parse CQL filter list. Function not found. Parsing :

texts:
Im trying to apply cross filter to retrieve the layer1 (where des= land) interects with layer2 as follows

http://localhost:2211/geoserver/wfs?service=WFS&version=1.0.0&outputFormat=application/json&request=getfeature&typename=ws1:layer1,ws1:layer2&CQL_FILTER=INTERSECTS(the_geom,querySingle('ws1:layer1','the_geom','des=land'))&srs=EPSG:4326

I looked at this question and this one. I'm still getting an error Could not parse CQL filter list. Function not found. Parsing : INTERSECTS(the_geom,querySingle('ws1:layer1','the_geom','descr_2=Urban')).

Where did I go wrong?


texts:
I have a really basic question which I can't wrap my head around.

I have a GeoTIFF of the world and a shapefile containing the oceans. The problem is that the GeoTIFF is centered at -180° longitude while the shapefile is centered at 0° longitude.

Looking at their SRS info using gdalsrsinfo it looks like they use exactly the same basic WGS84 projection.

How can I get them to overlap?


texts: Functionalities of expression builder in Autofield plugin available in QGIS 3.4?

texts:
Basically, what I try to do is automated calculation while editing (not with field calculator).

In QGIS 2.18 you could do automated calculations while editing with the Autofield Plugin. Since QGIS 3.X, you can do automated calculation of geometry (while editing) via the Properties > Standard values.

However, with the Autofield plugin you could also work with expressions based on existing field values. It seems that this is not possible in QGIS 3.4.: for some reason the fields and values don't show up in the expression window (Properties > Standard values), only NULL is available:

Does anybody know what I am missing, or if this functionality still exist in QGIS 3.4?


texts:
I'm not very experienced in this.

I am atmospherically correcting Landsat 8 images with ATCOR module in PCI Geomatica. I see that my final output reflectance values are not between 0.1 to 1.0, indicative of reflectance values. Also I read PCI manual and realized atmospherically corrected images are in percentage. I feel my reflectance values should range between 0.1 and 1.0. Any idea on this please (or can I just divide my image by 100 - i.e. converting percentage to decimal)?


texts: Where can I find the coordinate system used by GPS

texts: How do you determine which geographic coordinate system is used by a GPS device? For example, my Garmin GLO is outputting NMEA data with lat/lon. What system is represented by those lat/lon values? I think I heard the answer is WGS84 but how do I know / verify this? Is it just assumed all GPS NMEA data, unless otherwise projected, is WGS84?


texts: I created the slope layer in QGIS through "Raster > Analysis > Slope" and then used the function in GRASS r.slope.aspect, but how do I highlight a slope of 0° degrees to 60° degrees? I tried the input through the raster calculator 0 < slope < 60 but I obtain a completely black raster without any units.


texts: How do you unregister a table from a esri geodatabase?

texts:
How do you reverse the "Register with Geodatabase" operation?

I reckon it isn't uncommon to accidentally register (since it is the first option on Manage menu in ArcCatalog 10.2).

Is there a tool (or python function) that can reverse this?


texts: How to convert SPL2SMAP_S data (available in .h5 files) to *.tiff files?

texts:
I had downloaded some data from https://nsidc.org/data/SPL2SMAP_S
The description of the data is given in the page.
Using an earthdata login, one can download the data.
The data consists of *.xml and *.h5 files.

Before downloading the data i selected the region as is seen in the picture above.

I tried

gdal_translate -of GTiff HDF5:"SMAP_L2_SM_SP_1AIWDV_20200103T012207_20200103T010733_073E32N_R16515_001.h5"://Soil_Moisture_Retrieval_Data_3km/water_body_fraction_3km out.tiff

This output tiff when projected in QGIS was off

I am looking for water_body_fraction_3km actually.
I need help in creating a georeferenced tiff file (raster) containing this value.


texts: Error in Select Layer by Location python script?

texts:
I keep on getting an error on Select by location that the file cannot be found and the inputs are invalid, but I can't figure out what I'm doing incorrectly. I realize some of my variables unnecessarilly repeat themselves, which I plan to clean up, but it still doesn't explain why SelectByLocation isn't working.

# Import arcpy module
import arcpy

#Set Geoprocessing Environments
arcpy.env.overwriteOutput = True
arcpy.env.scratchWorkspace = "C:/PATH HERE/Testing"
arcpy.env.workspace = "C:/PATH HERE/ArcGIS/gdb.mdb"

# Parameters.
PRISM = arcpy.GetParameterAsText(0)
PRISM = "C:/LOCATION HERE"

#Unchanging Variables and Workspaces
Map = "C:/LOCATION HERE"
workspace = "C:/PATH HERE/gdb.mdb"
scratchWorkspace = "C:/PATH HERE/Testing"
PRISM_pts = scratchWorkspace + "/PRISM_pts.shp"
Pts_Lyr = scratchWorkspace + "/Pts_Lyr.lyr"
Map_Layer = scratchWorkspace + "/Map_Layer.lyr"
Prism_Lyr = scratchWorkspace + "/Prism_Lyr.lyr"

# Process: Raster to Point
arcpy.RasterToPoint_conversion(PRISM, PRISM_pts, "Value")

#Process: Make a Feature Layer for the Rasters
arcpy.MakeFeatureLayer_management(PRISM_pts, Pts_Lyr)

#Process: Make a Feature Layer for the Map
arcpy.MakeFeatureLayer_management(Map, Map_Layer)

# Process: Select Layer By Location
arcpy.SelectLayerByLocation_management(Map_Layer, "INTERSECT", Pts_Lyr, "", "")

#====================================
# Runtime error  Traceback (most recent call last):
# File "<string>", line 1, in <module>
# File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py",
# line 6618, in SelectLayerByLocation
# raise e ExecuteError: Failed to execute. Parameters are not valid.
# ERROR 000368: Invalid input data.
# ERROR 000732: Selecting Features:
# Dataset C:/PATH HERE/Testing/Pts_Lyr.lyr
# does not exist or is not supported
# Failed to execute (SelectLayerByLocation).
#
# =====================================

# Process: Spatial Join
arcpy.SpatialJoin_analysis(Map, Prism_Lyr, Spatial_Join_shp, "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "INTERSECT", "", "")

I don't quite understand why ArcGIS can't find the layer by location file, since it just created it. Am I missing a step?


texts: Creating waypoint coordinates on a grid within an area

texts: Do you know which GIS software or online tool supports returning GPS coordinates on a grid within an area? Kinda like waypoints for drone photographing 90 degrees looking down of an area for 3d modelling.