texts: Property with null-value disappearing from feature when creating feature-collection in Google Earth Engine?

texts:
When I create a feature with a property that has a null-value, that property disappears when adding the feature to a feature-collection. Why does this happen and how to prevent this?

An example (https://code.earthengine.google.com/5dd388c75939997f4d50d0e6628e1c22):

var dict1 = ee.Dictionary({'prop1': null, 'prop2': 400})

print(dict1);

var ft1 = ee.Feature(null, dict1)

print(ft1)

var ftc1 = ee.FeatureCollection(ft1)

print(ftc1)

"dict1" has two items, one of the values is "null". Creating a feature with this dictionary gives a feature "ft1", that has two properties. When adding "ft1" to a feature-collection, the property that has a null-value disappears from the feature.


texts: Create empty geometry in QGIS with Python

texts:
I would like to create a model point dataset in QGIS 3.0.3 with Python in the QGIS python console to do some testing. Specifically, I  want to create a point dataset that has a feature with an empty geometry, and, another feature with a null geometry.

What I have tried is:

layer = QgsVectorLayer('Point?crs=epsg:4326&field=id:string&field=type:string','ModelPoints',"memory")
pr = layer.dataProvider()

points = {}
points['Pt1'] = ['(1,1)', QgsPointXY(1,1)]
points['Pt2'] = ['(2,1)', QgsPointXY(1,1)]
points['Pt3'] = ['(4,1)', QgsPointXY(1,1)]
points['Pt4'] = ['empty', QgsPointXY()]
points['Pt5'] = ['null', None]

for id, point in points.items():
    pt = QgsFeature()
    pt.setAttributes([id,point[0]])
    if point[1] is not None:
        pt.setGeometry(QgsGeometry.fromPointXY(point[1]))
    print(pt.geometry().asWkt())
    pr.addFeature(pt)

layer.updateExtents()
QgsProject.instance().addMapLayer(layer)

In this code, the empty point geometry is intended to be created by:

pt.setGeometry(QgsGeometry.fromPointXY(QgsPointXY())

However, the result is:

pt.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(0,0))

So, my code is making a point with the vertex (0,0), but, whereas I need a point (or a geometry) with no vertices.

Note:  The approaches g = QgsGeometry.fromWkt('POINT') and g = QgsGeometry.fromWkt('POINT EMPTY') both result in a null geometry as shown by using print(g.isNull()).


texts:
I have points along along a line. Points that intersect the segment of the line to return attribute information. My line as three connected segments.

ESRI Docs says

  disjoint (second_geometry) Indicates if the base and comparison
  geometries share no points in common. Two geometries intersect if
  disjoint returns False.

with arcpy.da.UpdateCursor(in_point_along_line, ["SHAPE@", "Azimuth", "OID@"]) as cursor:
        for row in cursor:
            azimuth = get_azimuth_intersect(row[0], row[2])
            row[1] = azimuth
    def get_azimuth_intersect(in_geom, id):
        for geom in in_geom_azimuth_list:
            if in_geom.disjoint(geom[0]) is False:
                print id
                print in_geom.centroid
                print "Intersect"
                print geom[0].firstPoint
                print geom[0].lastPoint
                print geom[0].length

Result:

This one is reporting False, and is correct the point does intersect the line segment.
252501.7229 5653692.1437 NaN NaN The Point
252657.6907 5653761.9575 NaN NaN The line firstpoint
251107.9884 5653068.2849 NaN NaN The line lastpoint

The one is reporting it's also False? But it does not intersect...
252501.7229 5653692.1437 NaN NaN The Point
251107.9884 5653068.2849 NaN NaN The line firstpoint
251169.5741 5652467.8243 NaN NaN The line lastpoint

The one is reporting it's also False? But it does not intersect...
252501.7229 5653692.1437 NaN NaN The Point
251169.5741 5652467.8243 NaN NaN The line firstpoint
250415.1493 5651836.5709 NaN NaN The line lastpoint

What am I missing?


texts:
I am a novice when it comes to GIS systems and information. But, I am having a real problem trying to wrap my brain around a specific value that gdalinfo puts out.

Here is what gdalinfo shows from my GeoTiff file.

[bwalker@localhost tmp]$ gdalinfo Denver\ SEC\ 93.tif
Driver: GTiff/GeoTIFF
Files: Denver SEC 93.tif
       Denver SEC 93.tfw
Size is 16635, 12391
Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",38.66666666666666],
    PARAMETER["standard_parallel_2",33.33333333333334],
    PARAMETER["latitude_of_origin",37.5],
    PARAMETER["central_meridian",-107.5],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-380521.035404817666858,291529.854634510295000)
Pixel Size = (42.334666405531003,-42.334857259625188)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2015:06:05 10:25:35
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -380521.035,  291529.855) (111d57'15.54"W, 40d 2'50.37"N)
Lower Left  ( -380521.035, -233041.362) (111d41'25.51"W, 35d19'26.36"N)
Upper Right (  323716.140,  291529.855) (103d42'35.65"W, 40d 4' 8.14"N)
Lower Right (  323716.140, -233041.362) (103d56' 4.30"W, 35d20'39.77"N)
Center      (  -28402.448,   29244.246) (107d49'21.13"W, 37d45'47.67"N)
Band 1 Block=16635x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 255,255,255,255
    1: 255,255,0,255
    2: 255,0,255,255
    3: 255,0,0,255
....

I can not, to save my life, figure out what the "Output" line means?

Origin = (-380521.035404817666858,291529.854634510295000)

Does it have something to do with northing or easting? Can someone please point me to a reference that explains what it means. Or even explain it?


texts: Load gzipped HDF5 layer without decompressing using GDAL

texts:
I am interested in loading a VIIRS raster that's contained in a gzipped H5 file, but I get an error. I am using the following command to get the info:
gdalinfo HDF5:"/vsigzip/imagery/viirs/test3/VIIRS-SDR_VIIRS-Moderate-Resolution-Band-10-SDR_20200923_78694/SVM10_npp_d20200923_t0333337_e0339141_b46148_c20200923073914810642_noac_ops.h5.gz"://All_Data/VIIRS-M10-SDR_All/Radiance

But I get the following error:
ERROR 4: `HDF5:/vsigzip/imagery/viirs/test3/VIIRS-SDR_VIIRS-Moderate-Resolution-Band-10-SDR_20200923_78694/SVM10_npp_d20200923_t0333337_e0339141_b46148_c20200923073914810642_noac_ops.h5.gz://All_Data/VIIRS-M10-SDR_All/Radiance' does not exist in the file system, and is not recognized as a supported dataset name.

I know this path works because if I use it with the python bindings, I can load everything correctly. Also, the gdalinfo command works if I de-compress the h5 file first, but this is not practical for my application, and I would like to make it work with the compressed file.


texts: WMTS request with brutile and sharp map return wrong parameter

texts:
I'm using brutile and sharpmap to send wmts request and show maps.
here is my code:

    Public Shared Function Create() As ITileSource
    Dim uri As New Uri("http://91.98.30.14:8080/geoserver/gwc/service/wmts?SERVICE=WMTS&VERSION=1.1.0&REQUEST=getcapabilities&TRANSPARENT=TRUE&TILED=true")
    Dim req As Object = WebRequest.Create(uri)
    Dim resp As Object = req.GetResponseAsync()
    Dim tileSource As ITileSource
    Using stream As Object = resp.Result.GetResponseStream()
        Dim tileSources As List(Of ITileSource)
        tileSources = WmtsParser.Parse(stream)
        tileSource = tileSources.FirstOrDefault(Function(x) x.Name = "CH") '
    End Using
    Return tileSource
End Function

first of all, it is very slow to response and second, WmtsParser.parse returns "layer title" parameter for "layer name", and for "layer name" it returns "layer projection" to tileSources!
why do those issues happen?


texts: Replicating standard D3.js Conical US Albers projection in QGIS?

texts:
Oftentimes, when I see maps of the US rendered online, they use a custom version of the Albers projection, like so:

Alaska and Hawaii are off to the side, nicely arranged. I'd like to replicate this sort of thing in QGIS.

Would it be possible for me to do this?

Is this a D3-only projection?


texts: Outlining Kansas, adding in county boundaries and shapefile of oil and gas fields

texts:
I'm extremely new to ArcGIS Pro and still trying to understand the basics.
I need to map out all county boundaries in Kansas but I can't even figure out how to create a state boundary. After I create a state boundary and county boundaries, I need to add in a shapefile that I have on my computer that shows Oil and Gas fields around Kansas. After the fields are placed then I need to add in actual oil and gas wells using their coordinates (I assume that will be through the XY location). I am only focusing on the state of Kansas.
I've been trying to figure out how to do this on my own but have completely stuck.


texts:
I am trying to reclassify a slope (elevation) raster to a scale from 1 to 10 using QGIS GRASS r.recode with the following algorithm:

input_value_low:input_value_high:1:10

Unfortunately it will generate a layer with the scale from 1 to 9 and I have no clue why.
Any ideas?!


texts: How do you find the longitudinal and cross sectional profiles of polygons?

texts:
I am currently carrying out a research project involving submarine canyons. The canyons are mapped as polygons and I need to find the longitudinal and cross sectional profiles of the canyons. Each canyon also has a centre line, mapped as a polyline.
I know you can draw individual profile lines, but I have over 200 canyons and this would be extremely time consuming!

Does anyone know how it can be done by a more efficient method?


texts: Cannot add database connection using windows authetication SQL Server and ArcSDE 10.3

texts:
We can't seem to add a database connection from ArcSDE to sql server 2008 using ArcCatalog 10.3.1 from a specific computer with a specific username. When trying to add a database connection using this specific user's login (call it myuser), we encounter the error "Failed to connect to the specified server. Do you want to continue? Bad login". This occurs even when attempting to add a database connection using database authentication, but ONLY on this specified computer AND specified login

We have been trouble shooting this for the last day not and cannot seem to find a solution

To add more complexity to the problem, other registered users can add a database connection using their windows login on this computer and myuser can add a database connection using her windows login on other computers. Other users can also add a connection at the database level with the correct crudentials

Double checked the instance and database names for myuser
we have tried deleting and recreating the user from the instance level in sql
added sysadmin privileges to the user in sql server

We have tried just about everything but myuser is still unable to add a database connection on her computer.

Why is this?


texts: Choosing symbology classification scheme for Solar PV Electricity Generation?

texts:
I am presenting some data relating to rooftops suitable for solar PV and their potential PV electricity generation.

I am mapping the above data for the study area.

What would be the best symbology classification given the data? (Equal interval, Natural Break, Quantiles etc)


texts:
I have changed a shell launcher found at http://inasafe.linfiniti.com/html/id/developer-docs/platform_windows.html as to fit my wish: setting the system environment variables to be able to use PyQGIS modules from the default Python IDLE GUI. When I try to run the original version of the shell launcher it works, namely I succeed in importing (for instance) the "qgis.core" module, in DOS environment. When I launch my modified one (see below) the shell starts but the variables seem no to be set properly, as importing the module outputs the error:

ImportError: DLL load failed: Impossibile trovare il modulo specificato # (impossible to import the specified module)

My OS is Windows XP sp3, Python is 2.7.3 and QGIS is Lisboa v.1.8. Here's my personal SHELL LAUNCHER:

@echo off
SET OSGEO4W_ROOT=C:\QUANTU~1
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%OSGEO4W_ROOT%"\apps\Python27\Lib\idlelib\idle.pyw
@echo off
SET GDAL_DRIVER_PATH=%IDLE%\bin\gdalplugins\1.9
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass-6.4.2\lib
path %PATH%;"%OSGEO4W_ROOT%\apps\Python27\Scripts\"

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python;
set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\Python27\Lib\site-packages
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis

start "Quantum GIS Shell" /B "cmd.exe" %*

When I test whether if the variables are set properly, it seems they are, but still not succeeding in importing the modules.

import sys
for item in sys.path:
    print item

C:\Windows\System32
C:\Quantum_GIS_Lisboa\apps\Python27\Lib\idlelib
C:\Quantum_GIS_Lisboa\apps\qgis\python
C:\Windows\system32\python27.zip
C:\Quantum_GIS_Lisboa\apps\Python27\DLLs
C:\Quantum_GIS_Lisboa\apps\Python27\lib
C:\Quantum_GIS_Lisboa\apps\Python27\lib\plat-win
C:\Quantum_GIS_Lisboa\apps\Python27\lib\lib-tk
C:\Python27\ArcGIS10.1
C:\Quantum_GIS_Lisboa\apps\Python27
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages\PIL
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages\win32
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages\win32\lib
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages\Pythonwin
C:\Quantum_GIS_Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode


texts: pycsw and GeoNode — how to import metadata?

texts:
I am using pycsw as part of GeoNode 2.4.  I need to transfer metadata from an old app, which is powered by GeoNode 1.2 and GeoNetwork, but I can't find pycsw-admin.py file, which is supposed to import metadata.
Is this file named differently under GeoNode or is there some other way to import metadata?


texts: Mapbox personalized icons with cluster and filter checkbox

texts:
I have created a map using mapbox. I add icons from an Url and a type of "filter". But I want to add a marker cluster and in and also that the filter function with checkbox to select more than one option.

http://jsfiddle.net/MichelleCh/3du88qtr/7/

This is my example. I'll tried with the Leaflet marker cluster but... I'm totally new with code then, I don´t know in what place I should put the cluster code.

I think now, the most important form me is the checkbox for the filters. I tried to make put "true" or "false" to certain values in the properties of the icon and this at final. But aren´t checkboxes :(

myLayer.on('layeradd', function(e) {
    var marker = e.layer;
    var feature = marker.feature;

    var images = feature.properties.images
    var slideshowContent = '';

marker.setIcon(L.icon(feature.properties.icon));
    for(var i = 0; i < images.length; i++) {
        var img = images[i];
    }
    var popupContent =  '
' +
                            '
' + '
'+ feature.properties.title + '
'+'
' +
                            '
' + feature.properties.description +'
'

                        '
';

    marker.bindPopup(popupContent,{
        closeButton: false,
        minWidth: 400
    });

});

myLayer.setGeoJSON(geoJson)
.addTo(map);

$('.menu-ui a').on('click', function() {
      var filter = $(this).data('filter');
    $(this).addClass('active').siblings().removeClass('active');
    myLayer.setFilter(function(f) {

        return (filter === 'Development') ? true : f.properties[filter] === true;
    });
    return false;
});


texts:
I need to merge data from an OSM map with some other data (tiles) I get from a tileserver. I can import and display both layers well in QGIS and all zoom levels work fine. Now, I would like to export a region of this map as an .osm file so I can use it offline on my mobile phone (the app I am using on the phone can only load .osm files).

What I managed so far: I can export the merged maps from QGIS to an .mbtiles file. I can also extract that file to a folder structure or load the .mbtiles file as a single layer into QGIS again if that helps.

How can I get from the .mbtiles file to the .osm file, or more generally, how do I export the result of a merged OSM map with an tile-server map to single .osm file?

(additional info: I work with Ubuntu 14.04 and QGIS 2.18)


texts:
having build some (two) minor plugins, it's now time for me to link a few (perl-) applications to qgis. As there is no perl-interface in qgis I thought the best way to link both worlds together is a kind of messaging system. It should work like tis:

the perl-app sends a message to a certain folder
qgis starts a separete loop which looks every second, whether there is an incoming message in that folder
Qgis performs the necessary stuff an deletes the message again.
Or vice versa.

My question is how to start such a loop in qgis and how to do this at qgis-startup.

By the way - is there a possibility to "automate" QGis from external python applications? And if so is there any example code?


texts: Is it possible to edit a layer from Esri Personal Geodatabase in QGIS 3?

texts:
After a lot of research and effort I have successfully made a connection to my Esri Personal Geodatabase (.mdb file).

Now I want to edit some attribute data in the attribute table of a layer... BUT I can't because the toggle edit and the editing toolbar is greyed out.

So how do I edit a layer inside (.mdb) file that I have loaded into QGIS?!

I am using QGIS 3.

Don't ask me why I want to do it this way. (But If you wan to know why, it's because I want to edit the (.mdb)file in QGIS and then send it to someone who is using ArcMap who should be able to read/write to the file too. The other person has ArcMap and I have QGIS).


texts:
I have used the FeatureClassToNumPyArray tool in ArcPy, but the numpy array that is imported has only rows, with each of the fields in the original feature class as numpy.void type.
Here is a reproducible array.
b = np.array([(0,1),(1,1.5),(2,1.25),(3,1.28)],dtype=[('FID','<i4'),('Value','<f4')])

this returns a (4,) array, I want to reshape so that the array is (4,2) with the names 'FID' and 'Value', but can't quite figure out how to do this 'elegantly'.
I know that I can extract 'FID' and 'Value' as lists and then combine using np.array, but is that the best way to do it?
arr = np.array([b['FID'],b['Value']]).T


texts: Comparing layers in QGIS and finding overlay?s

texts:
I'm a newbie to QGIS. I have run into a thing I can not solve despite trying to search for a solution and working through the functions in the tool.

I have two layers I want to compare. One is filled with polygons, the other has polygons spread out. The first one is a map with real estate lots (pardon my English if this is not the right expression, what I mean is areas on the ground which is divided into different legal areas/entities). The second layer is houses.

What I want to have as a end result is a layer with the real estates where there are houses built upon them.

I've been looking for a function where I can compare the two layers and save all polygons from the first layer where there is a polygon in the second layer present. But so far I have only been able to find comparisons which cuts holes, find where there are intersections etc.


texts: Full list of objects not showing up after using Zonal Statistics as Table tool?

texts:
I am using the Zonal Statistics as Table tool in ArcMap 10.1 to compute the sum of raster cell values (brown scale in below image) in a given set of polygons (red in below image).

For an unknown reason, the output table is not including all the levels of the polygon file I am trying to get information about. I tried on a smaller set of polygons (as these polygons are buffers) and did not have the same problem. I have indeed asked the Zonal Statistics tool to ignore 'Nodata' (as shown in blue), so I am not sure why this is happening.
Suggestions?


texts: How to insert multiple images in the Atlas in QGIS?

texts:
I have a layer of CROSS SECTIONS (34) and the same number of IMAGES correspondent of each section. In the cross sections layer I have a column named PHOTOS where I defined this field as "ATTACHMENTS > IMAGES" in the attributes form.

The attribute table for the CROSS SECTIONS is this:

I made a ATLAS of it and I would like to see the images of each CROSS SECTIONS when I run the atlas. I tried to insert this HTML code:

  CROSS SECTIONS
  SECTION NUMBER =   [% "Cross-Sect"%]

to insert the images but it does not work as you can see in the image below:

Can you please tell me, how can I do this?
I´m using QGIS 3.4


texts: Leaflet vs Mapbox GL JS: JavaScript libraries comparison

texts:
I would like to do a comparison of Leaflet.js and Mapbox GL.js libraries. Therefore I would like to ask if you have some ideas of which parameters should I test.

In my initial research I decided to go with the following:

Learning aspects

Learning curve.
Frequency of updates
Community support
Documentation

Technical aspects

Performance : here I would like to compare for instance data consumption and loading times. (Does anybody know how to measure this in browsers?)
Interoperability : related with data format compatibility.
Compatibility : browser compatibility.
Raster and vector tiles.

Apart from my list I would like to know if you guys have any idea of what else can be compared in the context of mapping social related stuff (see the web apps at the end) for non-expert users with backgrounds different than computer science or geoinformatics.

The two web apps on which my question is based can be found here:

Leaflet
Mapbox GL JS


texts:
I'm using OpenLayers 2:10 Beginners Guide ebook. This code is APRA create a button click event in an opacity and zoom the map. But already rewritten several times and does not work. Nothing appears.

Where is the error?

<!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>
 <style type='text/css'>
    /*Custom Button*/
    .olControlCustomButtonItemInactive {
        background:#22dd22;
        border:5px solid #202020;
        cursor: pointer;
        height: 28px;
        width: 28px;
    }
 </style>
 <script type='text/javascript'>
 var map;
 //Função a ser chamada para criar o mapa na leitura da página
 function init() {
//Cria o objeto map
map = new OpenLayers.Map('map_element', {});
//Cria a primeira camada usando o serviço WMS, explicitamnete definida como base
    var wms_layer_map = new OpenLayers.Layer.WMS(
        'Base layer',
        'http://vmap0.tiles.osgeo.org/wms/vmap0',
        {layers: 'basic'},
        {isBaseLayer: true}
        );

//Adiciona as camadas ao mapa em uma martriz de objetos de camada
    map.addLayers([wms_layer_map);
 }

 //Função do botão para mudar opacidade e zoom ao evento do click
var custom_button_func = function(){
//Get a random coordinate from -90 to 90
var random_lon = Math.floor(Math.random() * 360) - 180;
var random_lat = Math.floor(Math.random() * 180) - 90;
if(map.layers[0].opacity === 1){
//If the layer opacity is 1 (fully opaque), then change it
and zoom
map.layers[0].setOpacity(.5);
map.setCenter(new OpenLayers.LonLat(random_lon,
random_lat), 3);
}
else{
//If the layer opacity is anything but 1, change it and
zoom
map.layers[0].setOpacity(1);
map.setCenter(new OpenLayers.LonLat(random_lon,
random_lat), 3);
}
};

//Cria o botão com a chamada para a função
var my_button = new OpenLayers.Control.Button({
    displayClass: 'olControlCustomButton',
    trigger: custom_button_func
});

//Painel onde será exibido o botão
var control_panel = new OpenLayers.Control.Panel({});

//Chama o addControl
control_panel.addControls([
    custom_button
    ]);

//Adiciona o painel ao mapa
map.addControl(control_panel);

//Move o mapa
control_panel.moveTo(new OpenLayers.Pixel(450,0));

 //Informações do centro do mapa
 if(!map.getCenter()){
    map.zoomToMaxExtent();
 }
 </script>
</head>
<body onload='init();'> <!-- Chama a função js init() -->
<!--Elemento HTML onde o mapa é exibido -->
 <div id='map_element' style='width: 500px; height: 500px;'>
 </div>
</body>
</html>


texts: List of ID's of all selected features as default-value

texts:
I would like to get a array with id's of all selected features of layer 'LayerB' (field-name: "id_of_layerB") as a default-value in a qgis-attribute-form of 'LayerA'.
I tried something like
aggregate(layer:='LayerB', aggregate:='array_agg', expression:= "id_of_layerB", is_selected('LayerB', ???))

but I have no idea how to realize this (especially what to put in in the "???"-part of the expression-part of the code.
I think, the "is_selected"-function checks agains just one feature, but is there a way tho get ALL (one or more) selected features of 'LayerB'?


texts: Making shapefiles converted from .KML appear in QGIS

texts: I have created shapefiles from Google Earth but they do not show in QGIS, whereas in ArcGIS Desktop they do. I have made sure they are at the top of the layers, not transparent and also I have changed their symbol colour. None of this has helped. I have also clicked 'Zoom to layer extent' and it zooms to where the layer should be. It is as if it is there, but not showing.


texts:
I have the following:

 var popProto = new OpenLayers.Protocol.WFS({
    url:  WFS_HOST,
    version: "1.1.0",
    featureType: layerName,
    featureNS: namespace,
    srsName: "EPSG:4326",
    propertyNames: ["grid_code"],
    defaultFilter: popFilter
});

which I wish to convert to GML so that I can chain it to a WPS request. I've been looking around for sometime, but I can't get a way out. Any suggestions on how I can convert this to GML format will be appreciated very much.

UPDATE:
Below is my attempt so far:

function createWPSExecuteRequest(namespace, layerName, popFilter) {
    var request = OpenLayers.Format.XML.prototype.write(new OpenLayers.Format.WPSExecute().writeNode('wps:Execute', {
        identifier: 'gs:Aggregate',
        dataInputs: [{
            identifier: 'features',
            reference: {
                mimeType: 'text/xml',
                href: WFS_HOST,
                method: 'POST',
                body: {
                    wfs: {
                        featureType: layerName,
                        version: '1.1.0',
                        featureNS: namespace,
                        srsName: "EPSG:4326",
                        propertyNames: ["grid_code"],
                        defaultFilter: popFilter
                    }
                }
            },
            identifier: 'aggregationAttribute',
            data: {literalData: 'grid_code'}
            .....more inputs here....
        }],
        responseForm: {
            rawDataOutput: {
                identifier: 'result'
            }
        }
    }));

    return request;
}

But the identifier: 'aggregationAttribute' and defaultFilter: popFilter is not being included in the generated GML. IS there anything I may be doing wrong?


texts:
So, I have two rasters showing land surface temperature where first is in range 18 - 27, and other 23 - 37. What I want is to make the same symbology for them so that these datasets can be comparable to each other. In this case they can't be compared since in first raster red color is showing values of cca more than 26, while on other raster values of 26-27 are represented by an orange color and so on for all the values. I would like them to be in stretched symbology since I have continuous data.

Can it be done?


texts: Cannot see GRASS Processing Error(s) in QGIS?

texts:
I'm trying to run r.neighbors in QGIS 3.4.1-Madeira and can't find the log that comes up with what looks like a red error message.  I've also tried a couple of other GRASS processing tools and they don't work so suspect it's a path/permissions problem of some kind.  This is on a Mac, High Sierra.

Trying to see the errors ... nothing in the logs via View - Panels - Log Messages.  Also under processing the box is checked to keep the dialog open after running an algorithm, but it still comes up then goes away before you can see what it says.


texts: Merge two maps only at the boundaries with python code

texts:
I have two maps one is normal cluster map and other is a grid map. I want that only at the boundaries the cluster map should follow grid map. Refer to the image, map in violet is my cluster map and map in green is my grid map. In this image grid map is on the outer part(1) but at some points grid map is on the inner part(2).

I am trying to do it in python using Geopandas but not getting any answer.The logic I am thinking is as follows The places where Grid map is on outer side I will take union and the places at which grid map is on inner side I will take intersection but i am not getting the condition logic .

Help me with a python code.


texts: Extract values to points works for only part of a raster

texts:
I am trying to extract values from multiple rasters using the "Extract Multi Values to Points" tool in ArcGIS 10.8.  The tool works until it gets to a certain row and then returns no data for the remaining rows.
Using the identify tool I have verified values for multiple point locations in the raster, double checked that all rasters and the shapefile are in the same projections, and the shapefile is not a multifeature point feature. I've tried the tool multiple times and fave even tried to extract the data from a single raster with the same results.


texts: Getting location of Milestones from Linear referenced data in ArcGIS Desktop?

texts:
I have shapefile of highways, where each segment has the information indicating the starting mile and ending mile.

Using these fields, I have created Routes, using the 'Create Route' gp tool from the ArcGIS Linear Referencing Toolbar.

From this route, I want to create a shapefile/featureclass indicating the location of each milestone (placed at the integer miles; i.e 5, 6, and so on).

How do I do this?

The only way I could find was to create a dummy table and then use the 'Make Route Event Layer'; But this requires to enter all the values in the table.

Is there a better way to do this?


texts:
I've been trying to run a viewshed analysis in ArcGIS 10.1 and everytime I've tried running it, it crashes the whole program, with no error. See the image below.

I've made sure it has OFFSETA field in right format, that it is one point only in the shapefile, and also that 3D analyst is ticked on.

Additionally I have also tried to run it as a Python script but it still produced the same result.

The elevation and XY distances are the same unit, in metres, though I haven't specified metres in the OFFSETA column. My input layers are a DSM raster that is quite big (85000x85000 roughly, 27.23 GB), and a point layer that is one point only. They are both in British National Grid projection. Im using ArcGIS 10.1 SP 1


texts:
I want to define projection for a Shape file. The projection system is:
"NAD_1983_StatePlane_Montana_FIPS_2500" which is referred in Qgis 1.9 as EPSG:102700.

The same projection system is not available in Qgis 1.8. Here it is designated as  "NAD 83/Mountana(ft)" and has EPSG:2256 and "NAD83/Mountana" and has EPSG:32100.

My problem is if I want to define projection in Qgis 1.8 what projection system should I use? EPSG:2256 or EPSG:32100?


texts:
I'm trying to install MAPNIK in ubuntu 16.04.
I downloaded the files from here: https://github.com/mapnik/mapnik, and follow the instructions from here:
https://github.com/mapnik/mapnik/blob/master/INSTALL.md
After running the ./configure option, just to check if I meet the requirements, I got the final message saying:
All Required dependencies found!

Overwritting and re-saving file 'config.py' ...
Will hold custom path variables from commandline and python config file(s)...

Configure completed: run `make` to build or `make install`

Now, when I try to build, I got this error:
scons: *** [srs/json/unicode_string_grammar_x3.o] Error 1
scons: building terminated because of errors.
Makefile:14: recipe for target 'install' failed
make: *** [install] Error 2

Any help with this error?


texts:
I made an insert into a table in this format.

insert into stg_place.stg_geo (data) values
(ST_GeomFromText('MULTIPOLYGON(((101.687764 3.134461,101.687582 3.131033,101.690898 3.132327,101.687764 3.134461)))',
                 4326));

This is a triangle and each of these points is in LONG LAT format.

Now, I was trying to confirm whether the polygon is correct. So , i did a test to find this long/lat in the polygon - 101.687764/3.134461 using the below query.

It didn't return any records. What is wrong in this ?

select * from stg_place.stg_geo where
ST_Within(ST_GeomFromText('SRID=4326;POINT('||101.687764||' '||3.134461||')'),data);


texts:
How to create vertical labels (each label in new line) like:

1998
150
pipe name

and not horizontal:

1998 150 pipe name


texts: How to fill the gap by using IDW(inverse distance weighting method) in R?

texts:
I'm beginner of R. I'm trying to fill the gap using idw, but I still had no idea about this.
The data in used is below:

The names sequence of columns are date/SO2/CO/O3/NO2/PM10/Lat/Lon. As you see, the gap values is existed in NO2 column. So I totally want to fill up this gap using idw through all location at the same date. It means like below picture num2.

Above picture has same date and different locations. So I want to use it to fill the gap by using idw.
I hope that you can understand what I'm trying to say (I'm Korean).
EDIT : Show my terrible code below
library(ggplot2)
library(rgdal)
library(dplyr)
library(tidyr)
library(gstat)
library(rgdal)
library(raster)
library(rgeos)
library(scales)
library(GISTools)
library(rgdal)
library(sp)
air_test <- read.csv("D:/aircondition/air_1Qu.csv", header = T)
str(air_test)
air_order <- air_test[order(air_test$측정일시),]
air_order <- air_order[,4:11]
a.fac <- as.factor(air_order$측정일시)
a.df <- data.frame(a.fac,air_order)
a.df <- a.df[,-2]
str(a.df)
names(a.df)[1] <-c("DATE")
a.df$ID <- as.integer(a.df$DATE)
str(a.df)
a.loc <- a.df[,7:9]
colSums(is.na(a.PM10))
str(a.df)

for (i in 1 : length(unique(a.df$ID))) {
  idw <- a.df[a.df$ID==i,]
    for (j in 2:6) {
    idw.air <- na.omit(idw)
    idw.air <- data.frame(idw.air[j], idw.air[,7:9])
    coordinates(idw.air) <- ~Lon + Lat
    proj4string(idw.air) <- CRS("+init=epsg:4326")
    idw.air_t <- spTransform(idw.air, CRS("+init=epsg:32651"))
    idw.final <- idw(formula = names(idw)[j] ~ 1, locations = idw.air_t,
                                      newdata = int.layer, idp = 2)
    idw.output.df <- as.data.frame(idw.final) %>% as.tibble()
    idw.output.df
  }
}

As you see, my codes are very mash. I really hope that the gap in my data is filled and I want to save the whole data which is fixed.
So, I totally want to be professional about R. Could you recommend something like books about R?
ERROR MESSAGES when I ran the code
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  :
  too many spatial dimensions: 226
In addition: Warning message:
In predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  :
  NAs introduced by coercion


texts:
I had to install "SemiAutomaticClassififcationPlugin" to complete one of my tasks. I tried to accomplish it via plugin manager but it ends up with error:

"Couldn't load plugin 'SemiAutomaticClassificationPlugin' due to an error when calling its classFactory() method".

Seems like it's not about the particular plugin, but about some more complex problem. Moreover, when I start my QGIS it always gives error "Python error: Couldn't load plugin 'processing'" but works fine, except couldn't run Python console. With the first error gives code that I attach below. Would be very thankful for any ideas about how to fix this.

Traceback (most recent call last):
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python\qgis\utils.py", line 335, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/misha_elyutin/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\__init__.py", line 65, in classFactory
    from .semiautomaticclassificationplugin import SemiAutomaticClassificationPlugin
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python\qgis\utils.py", line 685, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/misha_elyutin/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\semiautomaticclassificationplugin.py", line 43, in
    import numpy as np
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python\qgis\utils.py", line 685, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\misha_elyutin\AppData\Roaming\Python\Python37\site-packages\numpy\__init__.py", line 142, in
    from . import core
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python\qgis\utils.py", line 685, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:\Users\misha_elyutin\AppData\Roaming\Python\Python37\site-packages\numpy\core\__init__.py", line 23, in
    WinDLL(os.path.abspath(filename))
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\ctypes\__init__.py", line 356, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: [WinError 193] %1 не является приложением Win32

Python version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version: 3.4.15-Madeira Madeira, e83d02e274

Python Path:
C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python
C:/Users/misha_elyutin/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/misha_elyutin/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python/plugins
C:\PROGRA~1\QGIS3~1.4\apps\gdal2\pymod3\scripts
C:\PROGRA~1\QGIS3~1.4\apps\gdal2\pymod3\lib\site-packages
C:\PROGRA~1\QGIS3~1.4\apps\Python37
C:\PROGRA~1\QGIS3~1.4\apps\Python37\Scripts
C:\Program Files\QGIS 3.4\bin\python37.zip
C:\PROGRA~1\QGIS3~1.4\apps\Python37\DLLs
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib
C:\Program Files\QGIS 3.4\bin
C:\Users\misha_elyutin\AppData\Roaming\Python\Python37\site-packages
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\Pythonwin
C:/Users/misha_elyutin/AppData/Roaming/QGIS/QGIS3\profiles\default/python


texts: Clip raster by mask without change values

texts:
I am clipping raster by mask with python. My input raster owned Float32, however the output of clipping was Byte.

I don't know why? I suspect that transformation of tif to numpy.array is changing them.

I would like that output will have the same values from its origin image.

Here my code:

"""Clip a raster image using a shapefile"""

import gdal, gdalnumeric
import shapefile
import Image, ImageDraw

# Raster image to clip
raster = "/home/infra/PycharmProjects/map_ailton/dados/raster/resultado/custos_it_sema_2015xmunicipio_uc_multiplicacao.tif"

# Polygon shapefile used to clip
shp = "/home/infra/PycharmProjects/map_ailton/dados/vector/chapada/chapada_square"

# Name of clipped raster file(s)
output = "/home/infra/PycharmProjects/map_ailton/script/4558_06_Codes/output"

def imageToArray(i):
    """
    Converts a Python Imaging Library array to a gdalnumeric image.
    """
    a=gdalnumeric.numpy.fromstring(i.tostring(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

# Load the source data as a gdalnumeric array
srcArray = gdalnumeric.LoadFile(raster)

# Also load as a gdal image to get geotransform (world file) info
srcImage = gdal.Open(raster)
geoTrans = srcImage.GetGeoTransform()

# Use pyshp to open the shapefile
r = shapefile.Reader("%s.shp" % shp)

# Convert the layer extent to image pixel coordinates
minX, minY, maxX, maxY = r.bbox
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

# Multi-band image?
#Check this modification in script in: http://karthur.org/2015/clipping-rasters-in-python.html
try:
    clip = srcArray[:, ulY:lrY, ulX:lrX]

# Nope: Must be single-band
except IndexError:
    clip = srcArray[ulY:lrY, ulX:lrX]

# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

# Map points to pixels for drawing the county boundary
# on a blank 8-bit, black and white, mask image.
pixels = []
for p in r.shape(0).points:
  pixels.append(world2Pixel(geoTrans, p[0], p[1]))
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
# Create a blank image in PIL to draw the polygon.
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)
# Convert the PIL image to a NumPy array
mask = imageToArray(rasterPoly)

# Clip the image using the mask
clip = gdalnumeric.numpy.choose(mask, (clip, 0)).astype(gdalnumeric.numpy.uint8)

# Save clipping as tiff
gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)


texts: How to display information in a div outside the map?

texts:
I tried a bit around and created a simple map with two polygons(most of the code I got from the Interactive Choropleth Map Tutorial).
I added a "L.control()" to show information about the polygons in a div.
Contrary to the tutorial I would like to have this div outside the map.

To get a clue what I mean please look here: http://jsfiddle.net/Tbr9p/


texts: Can an ASCII Grid contain a CRS or projection?

texts:
I recently discovered the raster library in R. With that library one can define a raster's coordinate reference system (crs) or projection. For example, to set the projection as NAD83 UTM 10N:

crs(raster) <- CRS(+init=EPSG:26910')

I wrote that raster to file and imported into ArcMap, but Arc said that the file was missing a spatial reference system.

In general, is it possible for an ASCII file to contain/define its CRS or projection? If so, can that be done with R's raster library or with ArcMap?


texts: Taking an arbitrary object in Java and creating a shapefile out of it?

texts:
Let's say I have a data class in Java that essentially has

title
geo (WKT String)
Other irrelevant metadata

I want to create shapefile from lists of these types of objects using something like GeoTools. What are my options for being able to do this?

Essentially I'm creating an exporter for an internal data model we have.


texts:
I'm trying to create a hillshade surface from SRTM files downloaded form Earth Explorer, but the appearance is not the expected.

It has a grainy surface that I would like eliminate to have a smoothest surface. I'm running GRASS 6.4.3. to create this.

It happen mostly in coastal areas. I've tried to use r.denoise but the result is neither the expected, maybe isn't noise the cause of this.

Any idea?

If you need more info just ask.


texts: Transforming between ETRS89, WGS84, and Gauß-Krüger correctly using ArcGIS Desktop?

texts:
I have some projects where I work with this coordinate system:

  DHDN_3_Degree_Gauss_Zone_4
  WKID: 31468 Authority: EPSG
  Projection: Gauss_Kruger
  False_Easting: 4500000.0
  False_Northing: 0.0
  Central_Meridian: 12.0
  Scale_Factor: 1.0
  Latitude_Of_Origin: 0.0
  Linear Unit: Meter (1.0)

I usually get my data in ETRS89:

  ETRS_1989_ETRS-TM32
  WKID: 3044 Authority: EPSG
  Projection: Transverse_Mercator
  False_Easting: 500000.0
  False_Northing: 0.0
  Central_Meridian: 9.0
  Scale_Factor: 0.9996
  Latitude_Of_Origin: 0.0
  Linear Unit: Meter (1.0)

Or WGS84:

  GCS_WGS_1984
  WKID: 4326 Authority: EPSG
  Angular Unit: Degree (0.0174532925199433)
  Prime Meridian: Greenwich (0.0)
  Datum: D_WGS_1984
  Spheroid: WGS_1984
  Semimajor Axis: 6378137.0
  Semiminor Axis: 6356752.314245179
  Inverse Flattening: 298.257223563

Now I realised, that when I load my data into a ETSR89-project and to a Gauss-Krüger-project I have an offset of approximately 200 m between the two points. I used the transformation "DHDN_To_WGS_1984_5x" and "ETRS_1989_To_WGS_1984", respectively. (The "DHDN_To_WGS_1984_4_NTv2" does not give better results.) I found out about this using the calculate geometry function, which unfortunately calculates different WGS84 coordinates depending on the project coordinate system.

How can I make sure my data is on the right spot?


texts: Create PDF from map document (mxd) without Data Driven Pages

texts:
I've been trying to export map data (mxd) and combine into a PDF. I'm successful in creating the two PDFs, but can't combine them.

How do I correct this error?

  finalPDFDoc isn't defined

I'm very new to Python.

import arcpy
from arcpy import env
from arcpy.mapping import *

env.overwriteOutput = True
env.workspace = "C:\\ClassData\\Python_4080\\wk_8\\Lesson8_Data"

output_folder = "C:\\ClassData\\Python_4080\\wk_8\\Lesson8_Data\\Results"

#OPEN and work with map documents
mapDoc = MapDocument("C:\\ClassData\\Python_4080\\wk_8\\Lesson8_Data\\Lesson8.mxd")

#TODO define data frame

dataFrames = ListDataFrames(mapDoc)
dataFrame = dataFrames[0]       #gets first data frame to export

#output paths of exported maps. 1st is in data view and 2nd in layout view.
pdfPathdf = "C:/ClassData/Python_4080/wk_8/Lesson8_Data/Results/lesson8df"
pdfPathpl = "C:/ClassData/Python_4080/wk_8/Lesson8_Data/Results/lesson8pl"

#Export map in dataview, with data frame as third parameter
ExportToPDF(mapDoc, pdfPathdf, dataFrame, df_export_width = 1200, df_export_height = 800)

#Export map in layoutview, with constant as third parameter
ExportToPDF(mapDoc, pdfPathpl, "PAGE_LAYOUT", df_export_width = 1200, df_export_height = 800)
print "PDF has been generated at " + pdfPathdf

####### This doesn't run ###############

################################################################################
####################################################################################
#Add components to new PDFs. Merging the two views into one PDF

finalPDFDoc.appendPages(pdfPathdf)
finalPDFDoc.appendPages(pdfPathpl)
finalPDFDoc.saveACopy("C:\\ClassData\\Python_4080\\wk_8\\Lesson8_Data\\Results\\Final_PDF")

finalPDFDoc.close()
del finalPDFDoc


texts: NumPyArrayToTable - RuntimeError: unsupported time units. use M8[us]

texts:
I'm trying to convert a pandas dataframe (df) to a table in a file geodatabase using code below:

records = df.to_records(index=False)
arcpy.da.NumPyArrayToTable(records,r'C:\TEST.gdb\records')

Error message:

  RuntimeError: unsupported time units. use M8[us]

Using df.dtypes I can see that my date fields are type datetime64[ns]:

>>df.dtypes

field1              float64
field2              float64
field3              object
field4              float64
field5          datetime64[ns]
...

So I guess I need to change the type from datetime64[ns] to M8[us]. How do I do that?


texts: Basic WMTS example with GeoServer and OpenLayers

texts:
I want to set up a local WMTS with GeoServer. There are a couple of example data layers such as 'topp:states' which has the EPSG:4326 projection. I'm using the default configuration of GeoServer and testing it with OpenLayers like this:

import Map from 'ol/Map.js';
import View from 'ol/View.js';
import {defaults as defaultControls} from 'ol/control.js';
import {getWidth, getTopLeft} from 'ol/extent.js';
import TileLayer from 'ol/layer/Tile.js';
import {get as getProjection} from 'ol/proj.js';
import OSM from 'ol/source/OSM.js';
import WMTS from 'ol/source/WMTS.js';
import WMTSTileGrid from 'ol/tilegrid/WMTS.js';

const projection = getProjection('EPSG:4326');
const projectionExtent = projection.getExtent();
const size = getWidth(projectionExtent) / 256;
const resolutions = new Array(14);
const matrixIds = new Array(14);
for (let z = 0; z < 14; ++z) {
  // generate resolutions and matrixIds arrays for this WMTS
  resolutions[z] = size / Math.pow(2, z);
  matrixIds[z] = "EPSG:4326:" + z;
}

const map = new Map({
  layers: [
    new TileLayer({
      source: new OSM(),
      opacity: 0.7
    }),
    new TileLayer({
      opacity: 0.7,
      source: new WMTS({
        attributions: '...',
        url: 'http://localhost:8080/geoserver/gwc/service/wmts',
        layer: 'topp:states',
        matrixSet: 'EPSG:4326',
        format: 'image/png',
        projection: projection,
        tileGrid: new WMTSTileGrid({
          origin: getTopLeft(projectionExtent),
          resolutions: resolutions,
          matrixIds: matrixIds
        })
      })
    })
  ],
  target: 'map',
  controls: defaultControls({
    attributionOptions: {
      collapsible: false
    }
  }),
  view: new View({
    center: [-11158582, 4813697],
    zoom: 4
  })
});

But it shows me the map on the wrong location. What am I missing?


texts: Is the PostGIS documentation correct for ST_Project?

texts:
In example of the page ST_Project the select is

SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, 45));
                  st_astext
    ------------------------------------------
     POINT(0.63523102912532 0.63947233472882)
    (1 row)

But, in my enviroment this code give me an error like "Azimuth must be between -2PI and 2PI".

After a few tryings I solve with :

SELECT ST_AsText(ST_Project('POINT(0 0)'::geography, 100000, radians(45) ))

The sample in page tells the data type expected is in degrees, and the function expects radians.

May I using wrong PostGIS version?

Thanks.

BTW : My PostGIS version : 2.0


texts: I'm making a complex map with many labels (with QGIS 3.4.11). I manually arrange them for readability, attractiveness, etc.. However, I find that when I zoom in or out, or pan the view, their positions change. I don't want that. I want all of them to look the same (in the same positions), regardless of the current view. Under "Rendering" I do not have "Scale dependent visibility" checked; nor "Pixel sized-based visibility". I have checked "Show all labels..." Do I need a certain setting under "Data defined"?


texts: Arcgis Server - Highlight Identified Features

texts:
I am identifying features from an ArcGISDynamicMapServiceLayer

All the sublayers can be toggle on/off and when identifying, only the visible layers are queried/returned.  Below are the initIdentify and executeIdentifyTask functions.

   dojo.connect(map, "onLoad", function(map) {initIdentify(map);});

   function initIdentify(map) {
            //create identify tasks and setup parameters
            identifyTask = new esri.tasks.IdentifyTask("http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer");
            identifyParams = new esri.tasks.IdentifyParameters();
            identifyParams.tolerance = 3;
            identifyParams.returnGeometry = true;
            identifyParams.layerIds = wmsLayer.visibleLayers;
            identifyParams.layerOption = esri.tasks.IdentifyParameters.LAYER_OPTION_ALL;
            identifyParams.width = map.width;
            identifyParams.height = map.height;
            dojo.connect(map, "onClick", executeIdentifyTask);
        }

   function executeIdentifyTask(evt) {
            if (visibleLayers.length === 0){
                return;
            } else {
                identifyParams.geometry = evt.mapPoint;
                identifyParams.mapExtent = map.extent;
                identifyParams.layerIds = wmsLayer.visibleLayers;
                identifyTask.execute(identifyParams, function(idResults) { fillInfoWindow(idResults, evt); });
            }
  }

When the identify function is run, is there a way to highlight/select those features that are being identified?


texts: PostGIS queries result in "Symbol not found: __ZN5boost6system15system_categoryEv"

texts:
I seem to have lost the ability to use basic postgis queries. Normal postgres queries are unaffected.

Using Mac (High Sierra), using postgres 10.5 / postgis 2.4 (probably)

The same error happens accessing from QGIS, Postico, pgAdmin4, psql and using the Node.js pg library.

This seems to be a fairly new development - I successfully ran similar queries in Node about a week ago. I have no idea what's changed since then. I find it odd that the errors are related to SFCGAL as I don't recall ever installing this.

  PostgreSQL said: could not load library "/usr/local/lib/postgresql/postgis-2.4.so": dlopen(/usr/local/lib/postgresql/postgis-2.4.so, 10): Symbol not found: __ZN5boost6system15system_categoryEv
    Referenced from: /usr/local/opt/sfcgal/lib/libSFCGAL.1.dylib
    Expected in: /usr/local/opt/boost/lib/libboost_system-mt.dylib
   in /usr/local/opt/sfcgal/lib/libSFCGAL.1.dylib

I tend to install using Homebrew. Has anyone got a solution?

I suppose I can always set up a Docker container and use that, but I have a lot of databases that I'd rather not lose!


texts: X and Y Coordinates are not plotted correctly despite correct CRS

texts: I have converted a layer from MSSQL to X and Y points through QGIS geoalgorithms -> vector creation tools -> points layer from table. I set my target CRS as EPSG: 4326 (UK) but it does not plot it in the correct spot. Am I doing something wrong?


texts: Scaling Selection Symbols with Reference Scale using ArcGIS Engine?

texts:
With an ArcGIS Engine application, I am trying to figure out how to have the Selected Symbols scale with reference scale when selected using ArcObjects.  They are currently not scaling correctly in ArcGIS Engine.  Example, I select a stylized point that is using the "with this color" under layer properties (not the generic "with this symbol" option).

You can set this option in ArcMap via the "Select Options" dialog, and by checking "Scale selection symbols when a reference scale is set".  Saving the map does not persist this setting when an MXD is opened in the ArcGIS Engine application.  So in the example above, the stylized point will be colored, but it is only true to scale at the given reference scale when selected.

I have combed a lot of articles and found the various Selection interfaces (threshold tolerances, record counts, layer saving colors, etc.)

Where can I find this interface that controls the selection scaling (using ArcGIS Engine)?

I am about to go a graphics layer route, but I would prefer not to if the "easy" answer is out there.


texts: Ballpark transformation message-error for OSGB

texts:
This seems to be a repeat of previous CRS errors in QGIS 3.10+

New install of QGIS 3.10.5-A Coruña LTR on clean install of Linux Xubuntu 20.04 (install from Ubuntu-LTR PPA and apt). This OS install accesses the /home folder from a previous Xubuntu 18.04 install.

To create the error/problem:

Start a new project, QGIS defaults to EPSG:4326. No error messages.
Change project CRS to EPSG:27700 (British Ordnance Survey) and a warning message appears:

  Used a ballpark transform from EPSG:27700 to EPSG:4326

The same error occurs if I use EPSG:7405 (also British Ordnance Survey). Selecting random CRSs or custom CRSs does not produce an error message: the error appears to be unique to British Ordnance Survey projects.

Clicking on the error message Details gives:

  An alternative, ballpark-only transform was used when transforming coordinates between EPSG:27700 - OSGB 1936 / British National Grid and EPSG:4326 - WGS 84. The results may not match those obtained by using the preferred operation:

  Possibly an incorrect choice of operation was made for transformations between these reference systems. Check the Project Properties and ensure that the selected transform operations are applicable over the whole extent of the current project.

The project properties correctly state EPSG:27700. Moreover, this error occurs on a blank/empty project, therefore there have been no transformations attempted or performed.

Creating a template file set to EPSG:27700 and starting a new project from the template file produces the same error.

No other solution(s) suggested.

Question:
What is the practical effect of this error? I need to work to as high a level of accuracy as possible and certainly to better (less than) than 5m error. Do I need to worry about the error or not? Can I resolve the error?

QGIS 'About':

QGIS version
3.10.5-A Coruña

QGIS code revision
984615fe1e

Compiled against Qt
5.12.8

Running against Qt
5.12.8

Compiled against GDAL/OGR
3.0.4

Running against GDAL/OGR
3.0.4

Compiled against GEOS
3.8.0-CAPI-1.13.1

Running against GEOS
3.8.0-CAPI-1.13.1

Compiled against SQLite
3.31.1

Running against SQLite
3.31.1

PostgreSQL Client Version
12.2 (Ubuntu 12.2-4)

SpatiaLite Version
4.3.0a

QWT Version
6.1.4

QScintilla2 Version
2.11.2

Compiled against PROJ
6.3.1

Running against PROJ
Rel. 6.3.1, February 10th, 2020

OS Version
Ubuntu 20.04 LTS

Active python plugins
mmqgis;
LAStools;
GeoCoding;
quick_map_services;
processing;
db_manager


texts: Share FME workflows with other organizations

texts: Can we publish FME Workbenches to GitHub, to enable sharing these with other organizations.A custom transformer might be a solution?


texts: Setting alias for featurelayer name in ArcGIS Pro?

texts:
I have a layer in ArcGIS Pro that I want to separate technical name, and display name for the layer.
Reason for this is that the web client I develop has some functionality that fails if the layer name contains certain characters. So I need to have a technical name, and a different display name.

I see in the javascript code when debugging that the featurelayer object contains name and description, where the description contains the value that is set as alias on the featureclass table in the geodatabase.

I want to set this differently in the service when publishing this from ArcGIS Pro. How can I do this?


texts: Handling data gap in new SoilGrids datasets

texts:
I am working with the latest SoilGrids clay and OCS datasets, but notice a gap in Russia (180deg) when I export the datasets.
Do you have a recommendation how to fill this gap?
Did something happen when exporting?

To export I open the dataset in QGIS and export the dataset in Homosoline projection
When I import into GEE it shows the below gap, with values of 0.


texts: I installed pgAdmin 4.4. In the query tool I wanted to load a query that is saved on our server but this seems to be impossible because only the local drives c: d: etc. are provided. I tried to insert the IP address or the drive name e.g  x: but this doesn't work. What can I do?


texts:
I'm trying to write a python script with ArcPy that draws arcs and I've achieved something with creating a buffer and then clipping it with a bounding box so I only keep the length I need (I thought it was easier that way than calculating from my center point and radius, vertex by vertex from my start angle to my end angle). The problem is that the buffer too creates polygon buffers and since I only have ArcVIew (soon ArcEditor) and I can't use the polygon to line tool (licenced to ArcInfo).

Is there a way to convert a polygon to a polyline somehow with ArcPy or create directly a polyline buffer?


texts: Getting current windows username in pyQGis

texts:
I need to obtain the current windows user in a python script inside QGIS.

Is this possible somehow?


texts: Getting geometry of drawing in ArcGIS Server JavaScript API

texts:
To me, ArcGIS Help and Online Resources is always full of unneccessay information and missing some important points. For example there are lots of samples for drawing but I couldn't find any thing for steps post drawing like saving on Database or Map!

By the way, I am trying to Save some user inputs (Points or Multipoints) on the Map to MySQL database. I have a Table in database as :

create table Points (
name VARCHAR(20) PRIMARY KEY,
location Point NOT NULL,
description VARCHAR(200),
SPATIAL INDEX(location)
);

Now, I would like to save Geometry of the drawing point(s) in to this table. Can you please point me to an sample, tutorial or give me an idea how to get this information from the drawing classes?


texts: How to pass only selected features in python toolbox

texts:
I have a feature class in a file geodatabase. I want to select features from the feature class and export their attribute to an excel sheet. I need a python tool in ArcGIS which does this. My problem is, how can i pass only selected features to arcpy cursor.

In the post (Getting list of selected features in ArcGIS for Desktop using Python code?
) i read that cursor object only returns selected rows for a layer. But this does not seem to work for a feature class in file geodatabase. When i pass a feature class to search cursor it returns all the features in the feature class even when there are selected features in ArcMAP.

If you need more info feel free to ask.

code to pass the feature class to python toolbox

def getParameterInfo(self):
    param1 = arcpy.Parameter(
        displayName = "Input Feature",
        name = "in_layer",
        datatype = "DEFeatureClass",
        parameterType = "Required",
        direction= "Input")
    param1.filter.list = ["Polyline"]

   params = [param1]
   return params

def execute(self, parameters, messages):

    fc = parameters[0].valueAsText

    fields = [f.name for f in arcpy.ListFields(fc)]
    l = [[row.getValue(field)  for row in arcpy.SearchCursor(fc)] for field in fields ]


texts:
I run gdal2tiles on debian machine with 4 cores:

sudo /usr/bin/gdal2tiles.py -r cubic -a 0,0,0 -z 10-15 /home/adm/topo/ekb.vrt  /var/www/tiles

In htop I see:

Why does gdal use only one core in work?
Gdal version is 2.1.2.


texts: I am using QGIS 2.14 on Mac OS X 10.11.3; every time I try to open a raster file, in any format, QGIS crash, any suggestions?


texts: How to Count Features in a Feature Type in FME?

texts:
I am reading around 120 spatial tables in FME.  I would like to export a spreadsheet or CSV using FME that contains the feature name (table name) and the number of records in that table.

So far, I've tried linking a counter to each feature type and the Statistics Calculator to get the max value for the counter.  I then try to write the CSV with only fme_feature_type and max attributes.  I'm getting the max value output okay, but the fme_feature_type is not writing to the CSV.

Any ideas?  Also, do I need to add a counter to all the tables separately?  Or can I connect them all to a counter and then split it up somehow by fme_feature_type?


texts: gdal_calc raster calculcator : "any" syntax

texts:
I'm trying to combine 2 raster layers using gdal_calc in order to use them in QGIS 3.8.
One of the condition that has to be met could be translated as "if the value of raster A is 1 and the value of raster B is either 10 or 11 or... or 800 (I've got 20 different values or so), then keep the value of raster B".

I managed to achieve what I wanted using several logical_and :

logical_and(A==1,B==10)*B + logical_and(A==1,B==11)*B + ... + logical_and(A==1,B==800)*B

But it's tedious to write and hard to read for other users. Is there a way to get the same result using a more compact syntax ?
I tried using

logical_and(A==1,any(B==[1,10,11,12,33,35,41,42,43,44,51,52,53,222,223,224,602,700,702,710,800]))*B and logical_and(A==1,B==any([1,10,11,12,33,35,41,42,43,44,51,52,53,222,223,224,602,700,702,710,800]))*B

But neither of these gave me the same result as my previous lengthy expression...


texts: One vector data layer with different styles in QGIS

texts:
I have one layer of vector data that I would like to show in two different styles. Is there a way of accomplishing this in QGIS?

I've solved it now by duplicating the the vector layer and add separate styles for the two layers. That works, but there is one problem with this; When I edit an attribute in one of those two layers, the edited data does not appear in the other layer and the categorized symbols do not longer represent the updated data. I have to hit "Refresh" (F5) to make the edited data update in the other layer. This takes time, since there are many other layers in the project and is not a good solution.

Is there a better way of doing this? I'm pretty sure that this can be done in MapInfo and ArcGis where it's possible to have multiple themes (styles in QGIS) for one layer.

The layers are in WFS format.


texts: How can I add a customised Carto Map Marker via Carto Map Moblie SDK.UWP?

texts:
I'm implementing a Universal Windows Platform (UWP) app, and I m using the Carto Map Mobile SDK (UWP). However, I don't know how to add a .png image as a Map Marker programmatically. Here is my code:

        /// Preparation - create layer and datasource

        // projection will be needed later
        Projection projection = map.Options.BaseProjection;

        // Initialize an local data source - a bucket for your objects created in code
        LocalVectorDataSource datasource = new LocalVectorDataSource(projection);

        // Initialize a vector layer with the previous data source
        VectorLayer layer = new VectorLayer(datasource);

        // Add layer to map
        map.Layers.Add(layer);

        /// Now we real adding objects

        // Create marker style
        MarkerStyleBuilder builder = new MarkerStyleBuilder();
        builder.Size = 20;
        BinaryData iconBytes = AssetUtils.LoadAsset("Z:/FolderName/ProjectName/Assets/markers_mdpi/mapmarker.png");

        byte[] bytearray = iconBytes.GetData();
        int size = Marshal.SizeOf(bytearray[0]) * bytearray.Length;

        IntPtr pnt = Marshal.AllocHGlobal(size);

        builder.Bitmap = new Bitmap(pnt, true);

        MarkerStyle style = null;

        style = builder.BuildStyle();

        // Create a marker with the style we defined previously and add it to the source
        Marker marker = new Marker(position, style);
        datasource.Add(marker);

The Carto Map official technical document didn't help at all, and here is the screenshot . However, when I installed the official SDK.UWP via the Nuget, there aren't any relevant functions that mentioned in the document in the library.

Can anyone help me solve this problem?

Otherwise, it is meaningless for me to create this UWP app further.


texts: Can I made a legend permanent so that it doesn't disappear when the layer isn't visible?

texts: I have a map with four layers, and I want to have one legend for all the layers. If I add a legend to each layer, then when all four layers are visible, there will be four of the same legend. And if I add a legend to only one layer, say layer 1, then when layer 1 is hidden, so is the legend. Is it possible to create a legend associate with layer 1 and have it visible in other layers even if layer 1 is hidden?


texts:
I recently tried to install geopandas and shapely into my ArcGIS Pro Python 3.7 environment, but instead of utilizing the package manager within Pro, I foolishly tried to conda install into the base environment using miniforge3 (for Windows 10). I then tried to fix a dependency issue by uninstalling and reinstalling geos, but now my ArcGIS base environment is hopelessly broken, as it cannot seem to access numpy, mkl, or gdal.
Somehow, the process of installing/reinstalling these packages caused Python to upgrade from 3.7.9 to 3.7.10, and now my attempts to conda install --revision 0 are being met with a litany of PackagesNotFoundErrors. Pro tells me that my base environment is broken when I go to the package manager, but it won't allow me to make any modifications to the base environment.
Is there anything I can do short of a full reinstallation of Pro?


texts: How can I copy selected polygon features to a polyline feature class with ArcObjects?

texts:
I would like to copy selected polygon features to a polyline feature class. (Similar to how polygon features can be copied and pasted inside an edit session in ArcMap)

Until recently I've been using the pPolygon2:ExteriorRingCount, pPolygon2:InteriorRingCount method - iterating through the polygon rings and converting them into polylines and so on...  As described in this post I have a suspicion that the pPolygon2:ExteriorRingCount method might not be working in ArcObjects 10.1 (VS 2010 Express, VB.NET)

Although this approach has been working flawlessly until version 10.0 perhaps I am over-complicating things.  What is the best approach for copying polygon features to a polyline feature class with ArcObjects?


texts: Remove border of tiles while Smooth the image by convolving with the boxcar kernel in Google Earth Engine

texts:
I am trying to remove the black border of each tiles in Google Earth Engine while doing "Smooth the image by convolving with the boxcar kernel" on sentinel 2.

Here's my code https://code.earthengine.google.com/155b6b5e2a3429d198f33429b2ff6349

I am not sure what i am missing there, can someone please help me

// Load and display an image.
var image = ee.ImageCollection('COPERNICUS/S2').filterDate('2018-09-01','2018-09-20');
Map.setCenter(-121.9785, 37.8694, 11);
//Map.addLayer(image, {bands: ['B4', 'B3', 'B2'], max: 5000}, 'input image');

// Define a boxcar or low-pass kernel.
var boxcar = ee.Kernel.square({
  radius: 7, units: 'pixels', normalize: true
});

// Smooth the image by convolving with the boxcar kernel.
var smooth = function(img){
  return img.convolve(boxcar);
}
Map.addLayer(image.map(smooth), {bands: ['B4', 'B3', 'B2'], min:0, max: 5000, format: 'png'}, 'smoothed');


texts: Is it possible to run a model in ArcMap 10.3 and export directly to a MySQL table? Preferably avoiding ArcSDE as well.


texts: Spatial distance join between polygon and points using R?

texts:
I'm figuring out how to do a Intersection (Spatial Join) between point and polygons from shapefiles. My idea is to get the closest points and those points that match completely inside the polygons. In ArcGIS there's a function for match option named CLOSEST and they have defined by: "The feature in the join features that is closest to a target feature is matched. It is possible that two or more join features are the same distance away from the target feature. When this situation occurs, one of the join features is randomly selected as the matching feature."

I have a function to intersect points into polygons, it was kindly contributed by Lyndon Estes at the r-sig-geo list and the code works very well when all the polygons have filled all areas. The second case is known as a Spatial join distance and in ArcGIS is know as INTERSECT when match_option is CLOSEST, as ArcGIS does. So, you can modify the minimal distance between the point and the polygon when the area is not filled by all polygons, but I don't have access to ArcGIS.

Here's the data (direct to my public ftp) and the function of the first INTERSECT:

library(rgeos)
library(sp)
library(maptools)
library(rgdal)
library(sp)
xy.map <- readShapeSpatial("http://www.udec.cl/~jbustosm/points.shp")
manzana.map <- readShapeSpatial("http://www.udec.cl/~jbustosm/manzanas_from.shp" )

IntersectPtWithPoly <- function(x, y) {
# Extracts values from a SpatialPolygonDataFrame with SpatialPointsDataFrame, and appends table (similar to
# ArcGIS intersect)
# Args:
#   x: SpatialPoints*Frame
#   y: SpatialPolygonsDataFrame
# Returns:
# SpatialPointsDataFrame with appended table of polygon attributes

  # Set up overlay with new column of join IDs in x
  z <- overlay(y, x)

  # Bind captured data to points dataframe
  x2 <- cbind(x, z)

  # Make it back into a SpatialPointsDataFrame
  # Account for different coordinate variable names
  if(("coords.x1" %in% colnames(x2)) & ("coords.x2" %in% colnames(x2))) {
    coordinates(x2) <- ~coords.x1 + coords.x2
  } else if(("x" %in% colnames(x2)) & ("x" %in% colnames(x2))) {
    coordinates(x2) <- ~x + y
  }

  # Reassign its projection if it has one
  if(is.na(CRSargs(x@proj4string)) == "FALSE") {
    x2@proj4string <- x@proj4string
  }
  return(x2)
}

test<-IntersectPtWithPoly (xy.map,manzana.map)

Sharing some ideas with Lyndon Estes, he told me this:

I think the easiest thing to do would be to put a buffer around each of the points (you could specify 50 m if it is in projected coordinates), converting them to polygons, and then your task becomes an intersection of two different polygon objects.

I haven't done this type of operation in R, but I suspect you could find your answer with the following functions:

library(sp) ?over

library(rgeos) ?gBuffer ?gIntersects

Maybe someone else who has a better idea on polygon to polygon intersects/overlays could suggest the method.

I know that this functions could help to achieve it.

library(sp)
?over

library(rgeos)
?gBuffer
?gIntersects

On the other hand, I think perhaps making the points bigger could help??... but I'm not sure!


texts: Finding intersection point between two arcs?

texts: I am trying to find the intersection between two arcs. I have done this for line segments using shapely package. Does anyone have any advice?


texts: Export to PDF from ArcMap has Point Marker Symbols Appearing as Black Dots?

texts:
I am having trouble exporting MXDs to PDF.  If I export and MXD to PDF some point symbols appear as black dots (not the colored dots they should be) and others do not.  This happens in different data sets in the same map.

I have embed fonts and convert marker symbols checked in the Export to PDF tool. I do not use non-USA thousand separators.

I have read this ESRI literature about diagnosing export problems:

http://support.esri.com/es/knowledgebase/techarticles/detail/17783

Based on that document link above I can export the MXD to EMF, open the EMF in MS Paint and the EMF looks correct.
IF I add the EMF to an MS Word document all the point symbols appear as crescents!?

I had the IT department rebuild my Citrix profile, reinstall the PDF drivers, and I did a reinstallation of ArcGIS 10.3 for Desktop.  The first map exported after all this work looked fine, but all subsequent exports revert back to the black dot problem.


texts:
I would like to add a citation on the map (like for example the "powered from esri")  that would mention the researchers that created the displayed data, so if someone scrolls or zooms in out, it will be visible at all times.

Is it possible?


texts: Importing OSM data to pgRouting enabled database including turn restrictions

texts:
I recently used the tool pypgroutingloader to import OSM data into a pgRouting database. For a small amount of data, this tool works perfectly. Most routing restrictions are present and it modifies the costs of various different road types right during the import. After adjusting the underlying speed profile, I was quite impressed by the results of the pgr_trsp() function!
Then I tried to load some more OSM data into a database (a whole country). The python-script easily eats up all my RAM (16 GB) and crashes within minutes. This happens even for medium large parts. Furthermore it is not possible to execute the script incrementally using smaller parts. The script is not able to merge nodes with the same OSM ID. It internally assigns new unique IDs (starting from gid = 1 every time). Since I'm completely new to python programming, I'm not able to rewrite the script in a way to make this possible.
Is there any tool with similar capabilities? (turn restrictions + "speed profile")
If not, how could one achieve similar results using other tools, like osm2pgrouting?


texts: Converting polygons to/from raster coordinates with PostGIS

texts:
I'm looking to convert polygons from one set of imagery to a another (imageryA -> world -> imageryB) and I'm struggling to come up with an in-database solution. I have a solution in python, but we would really benefit from something that could be hooked into a trigger.

The ST_WorldToRasterCoord/ST_RasterToWorldCoord functions look like the perfect solution, but I'm having trouble getting them to work with polygons.

Right now I'm trying to use ST_DumpPoints to unpack the polygon but I'm having difficulty sticking the points back together as my array_agg causes an "aggregate function calls cannot contain set-returning function calls" error. Here's what I have so far trying to compute the world-coordinate polygon; I'm sure it's not pretty:

select poly_set_a.id,
       ST_MakePolygon(
           ST_MakeLine(
               array_agg(
                   ST_MakePoint(
                       ST_RasterToWorldCoordX(imagery_a.raster,ST_X((ST_DumpPoints(poly_set_a.image_geom)).geom)::int),
                       ST_RasterToWorldCoordY(imagery_a.raster,ST_Y((ST_DumpPoints(poly_set_a.image_geom)).geom)::int)
                       )
                   )
               )
           ) as world_poly
    from poly_set_a inner join imagery_a on (poly_set_a.image_id = imagery_a.id)
    where poly_set_a.id=1234;

Rasters seemed like a great solution but I can't find anything about using them with polygons. Maybe breaking it down into a function would be better.

Is it going to be easier to cobble the transform together from ST_Transform and ST_Scale?


texts: Installing SciPy for Semi-Automatic Classification Plugin of QGIS?

texts:
When I open QGIS, this message shows:

  Semi-Automatic Classification Plugin: Please, restart QGIS for
  executing the Semi-Automatic Classification Plugin. Possible missing
  dependecies: SciPy

When I try to reinstall or uninstall the plug in, or even only click on it, it always crashes.
I have already reinstalled QGIS but didn't help.
I have OS X 10.8.5 and QGIS 2.14.8

I have to do a land cover change analysis for a project and wanted to use NDVI, I had all my hopes on this plug in.


texts: Will ArcGIS Pro have (better) ModelBuilder?

texts:
An earlier Q&A (What are Esri's new "ArcGIS for Professionals" and "ArcGIS Pro"?) provided some reassurance that ArcPy code would have re-use value with ArcGIS Pro:

  ArcGIS Professional Extensibility. Here are our current thoughts about extensibility for ArcGIS Pro. You can write and run Python
  scripts in ArcGIS that call geoprocessing tools and use an exhaustive
  suite of scripting functions available in the ArcGIS Python API,
  ArcPy, to automate your GIS tasks.

However, it makes no mention of whether ArcGIS Pro will enable model building to automate workflows using a GUI that is either (equivalent to) the existing ModelBuilder or a major improvement on it to address previous enhancement requests from others and I.

Has anyone heard the road ahead for ModelBuilder (or at least model building) in ArcGIS Pro?


texts: Problems coloring features of file geodatabase in QGIS

texts:
I managed to open a file geodatabase (.gdb) in QGIS. The extensions of my data are correctly shown and the attribute table shows all values correctly. However, when I double-click on the layer to color a specific column based on its values it does not classify the values as it used to with normal shapefiles. I attached a picture to show what I mean.
Opening the attribute table and clicking on a random row does highlight the specific feature on the map so I assume the data is read in correctly.
Any idea on what to improve?
I use Qgis 2.8 and gdal 1.11.2 on Ubuntu 14.04


texts: Is there a way to get correctly aligned print export with OpenLayers plugin?

texts:
I am a beginner using QGIS (1.8.0 Lisboa) on a Mac (OS X 10.6.8 Snow Leopard). I work on a file with two layers in the following sequence (top to bottom):

a vector polygon layer with a few outlines of various buildings in France and Italy (CRS: EPSG:4326 - WGS 84 - "on the fly" re projection is activated )
a Google satellite layer uploaded via open layer plugin.

Sometimes when changing the scale (scrolling in and out) the Google layer updates a little wrong in size or position compared to the polygon layer. But as soon as I "pan" the map everything reconstructs right.
Now the problem I have concerns the print composer.

The result of any kind of export - pdf, or image or svg - of my map (vector layer + Google layer in the background) differs from what is shown in the preview of the print composer:
the vector layer and the Google layer positions are mismatching and in a different scale. I can send some pictures showing screenshots of the print composer preview and the final results if required.

I tried almost everything:

new files with different CRS settings
"cache" update preview
"render" update preview
different paper sizes,
different scales,

The problem appears only using an open layer map (Google, OpenStreetMap...). If I try the same using two vector layer (shapefile) - one with the buildings, another one with the city in background - everything works out great.

...exasperating! Since it is a real important research project I would be very grateful for some help and tips.


texts: Table of Latitude and Longitude with Associated Elevation

texts: I need to create a table that allows me to retrieve the elevation value for a given GPS coordinate pair. Currently, I am just doing this for New Zealand. I'm working in QGIS using a GeoTIFF file, and have already extracted the contours from the raster. This has given me a table of elevation values for the contour id's, but I need a table of lat/long pairs and their elevation values. Is there a way to extract this data from the GeoTIFF as a table?


texts: Converting coordinates: Feet to Degrees-Minutes-Seconds

texts: I have a point shapefile and I want to get the coordinates for each point in Degrees-Minutes-Seconds format. When I use Add XY Coordinates (Data Management) it gives me the coordinates in Feet format. Any idea/ thought?


texts: Custom GeoServer rendering transformation

texts:
Following on from my previous question here, I'm trying to create rendering transformation using the GeoScript Groovy API. I've got most of the code written but I can't figure out how to wrap it up for GeoServer to discover.

I've found the process package in the GeoScript API but, to me, it looks like it's for an existing process. I've also looked at this Python example, but the way the inputs and outputs are defined is different to Groovy. Could someone please provide a short example?


texts:
I've created a line by mouse click event when I click on the second point I need two lines at an angle of 28 degrees from the first line how should I do it?.

IPoint point1;
        IPoint point2;

     protected override void OnMouseDown(ESRI.ArcGIS.Desktop.AddIns.Tool.MouseEventArgs arg)
    {

        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
        IGraphicsContainer gContiner = mxdoc.FocusMap as IGraphicsContainer;
        IDisplayTransformation dispTransformation =
        mxdoc.ActiveView.ScreenDisplay.DisplayTransformation;

        if (point1 == null)
        {

            point1 = dispTransformation.ToMapPoint(arg.X, arg.Y);
            IMarkerElement markerElement = new MarkerElementClass();

            IElement element = (IElement)markerElement;
            element.Geometry = point1;
            gContiner.AddElement(element, 0);
            mxdoc.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics,
            null, null);
        }
        else if (point2 == null)
        {

            point2 = dispTransformation.ToMapPoint(arg.X, arg.Y);
            IMarkerElement markerElement = new MarkerElementClass();

            IElement element = (IElement)markerElement;
            element.Geometry = point2;
            gContiner.AddElement(element, 0);
            mxdoc.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics,
            null, null);
        }

        else
        {
            //
            IPolyline polyline = new PolylineClass();
            IPointCollection pointColl = polyline as IPointCollection;

            pointColl.AddPoint(point1);
            pointColl.AddPoint(point2);

            // Create a line element, this is the graphic that will get added to container

            IElement element = new LineElement();
            element.Geometry = polyline;

            ILineElement lineElement;
            lineElement = element as ILineElement;

            gContiner.AddElement(element, 0);
            mxdoc.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);

        }
    }


texts: Checking whether lat/long is in sea using Google Maps API?

texts: I am using the Google Maps API, in which as of now we are plotting the track of a vessel track, but the system is now allowing us to update even a Lat/Long which is on land. Is there is any way to stop this and check that the system should allow only updates to a lat/long that is in the sea?


texts:
I am just starting to look into the QGIS project.

I wish to plan a roadbed on my property. I would like to start with a map with topographic information (height, direction, distance) and draw the roadbed on the map.

I would then like to be able to read the gradient and distance of the planned roadbed as well as curvature and curve radii. Maybe also calculate the amount of mass to be removed or filled to support a roadbed.

Would QGIS be a suitable tool to do this?

If so, I would like to get involved in the project and learn to use the program.

At this stage I am not even sure where I could find a response to my question, if any.


texts: clip and/or intersect tools giving me the wrong results in qgis 2.18.16

texts:
I have my own digitilized polygon (called clipRahmen, north Asia), I wanted to clip all vector layers (points, polygons, lines) within this polygon which is an area north Asia, I used "clip" and "intersect", even the tool from the processing toolbox (under Geoalgorithms, Vector, Polygons) "Polygon intersection", but they are giving me the same result, which was the result NOT within this polygon that I have created. Data that I'm using is downloaded from Natural Earth, corner coordinate is user-defined: 100 000, layer coordinates WGS84.

In Screenshot u can see that the yellow points (populated places, 10 Mio. -Natural Earth data) are clipped/intersected, violet points are original-not clipped, but are not within the Polygon, blue rectangle is my polygon (north Asia). Is there any way to get correct clip result without updating the qgis?

I used the same tool (clip) for same data (points, lines, and polygons) in the same version of QGIS (2.18.16) for another area- was working (clipping and clipped layers were WGS84, corner coordinate was 100 000: user-defined), now is the same, I've just clipped to different area (north Asia), but is not working.

result:


texts: Creating scene layer package in ArcGIS Pro?

texts: I have a multi-patch layer in ArcGIS Pro 2.0 and I want to share the layer to ArcGIS Online. I cannot figure out how to save the multi-patch layer as a scene layer package. Under ArcGIS Pro Share tab, I can only see project, map, layer and geoprocessing package.


texts:
I exported a model as script using the properly option in QGIS 3.8.
The Script works perfectly, but the output name is changed.

Heres a print of the model:

I want that the script´s output name, from the Processing Toolbox, be the one that I designed in the  def initAlgorithm section, and not the given default name.

In the model, the output name comes as "output_name"
and in the script, the output name comes as "Reprojected".

Here's the converted Script:

    from qgis.core import QgsProcessing
    from qgis.core import QgsProcessingAlgorithm
    from qgis.core import QgsProcessingMultiStepFeedback
    from qgis.core import QgsProcessingParameterVectorLayer
    from qgis.core import QgsProcessingParameterFeatureSink
    from qgis.core import QgsCoordinateReferenceSystem
    import processing

    class Rename(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('shapeinput', 'ShapeInput', types=[QgsProcessing.TypeVector], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Output_name', 'output_name', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
        results = {}
        outputs = {}

        # Reproject Layer
        alg_params = {
            'INPUT': parameters['shapeinput'],
            'TARGET_CRS': QgsCoordinateReferenceSystem('EPSG:4326'),
            'OUTPUT': parameters['Output_name']
        }
        outputs['ReprojectLayer'] = processing.run('native:reprojectlayer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Output_name'] = outputs['ReprojectLayer']['OUTPUT']
        return results

    def name(self):
        return 'rename'

    def displayName(self):
        return 'rename'

    def group(self):
        return 'rename'

    def groupId(self):
        return ''

    def createInstance(self):
        return Rename()


texts: Find records with highest value based on unique ID

texts:
I'm trying to find a way to return all records that have the highest value in a field for each set of unique IDs. The attached table shows what I'm trying to achieve. The highlighted records would be the ones I want to have selected so I can export them to a new table.

For each set of unique IDs (Int_ID) I want to select those that have the highest value in (MPO_Lanes). I need to keep all the attributes, so Case items in Summary Statistics does not appear to work.

If your curious why I want to do this...for each Street Intersection (Int_ID), I want to identify which Street Legs (ExtID) have the most lanes (MPO_Lanes) and therefore identifying which Street Legs would be considered the "Major" street at that intersection.

I'm using ArcGIS 10.4.1, advanced.


texts: Finding details of both overlapping shapes using ArcGIS Pro without Advanced level license

texts:
After Understanding Union analysis output using single feature class with overlapping polygons in ArcGIS Pro I now have a better understanding of that tool.
However, I have a single polygon Feature class with about 33,000 shapes. I'm looking for overlaps.  I want to be able to say "Shape x overlaps with Shape y and Shape z"
Union (Analysis) doesn't do this. Intersect (Analysis) which may requires an Advanced level license but I am using ArcGIS Pro 2.6 with a Standard level license.
How else can I do it?
I've looked at this solution by NeilAyres: https://community.esri.com/t5/data-management-questions/find-overlapping-features/td-p/6948
In his code he says
poly_dict = {}
# union has "ORIG_FID" as part of its output
with arcpy.da.SearchCursor(fc, ["ORIG_FID", "SHAPE@", "SHAPE@XY"]) as Cur:

Mine doesn't. There is no "Orig_FID" in my union analysis output. So the rest of the solution is meaningless. Maybe I'm using different software?
On that same page, KenBuja says to try the "Find Overlapping Features tool" which seems to be a plugin. But not one that I can get working with ArcGIS Pro.


texts: Add spatial index to versioned feature class

texts:
I try to add spatial index to versioned feature class (Oracle SDE) using toolbox "Add Spatial Index" in ArcGIS 10.6 . When I try it as Data owner I have error message -"Tool must be run as owner of SDE.Default version". But when I try to run it as SDE user I have error message "Tool must be run as Data owner".

So I am little bit confused. Do you have any idea how to solve this problem?

I can do this in ArcCatalog as Data Owner, but I need to run via toolbox as batch process.


texts: ArcScene Extrusion Is not proportioned correctly?

texts:
I am new to ArcScene so I have no idea what I am doing.

In ArcGIS Desktop I have a buildings layer that I would like to make into 3D.

I created a field in ArcGIS Desktop with the building heights and used that field as the Extrusion value.

When I press apply the building are like really high and obviously not correct. I tried doing /10, /100, and even /1000, but then building either were still too high or were 2D again.

I do not know what to do at this point. I am using GCS_WGS_1984 as my coordinate system.


texts: Showing only outer boundaries of polygons / features in QGIS

texts:
I started with the Natural earth admin 0 - Countries vector layer. Based on this, I created  a new layer (using the 'Collect Geometries' tool) that only contains 26 features, representing major world regions (colors in the figure show regions; in this example North Africa has been selected and is shown in yellow).
Now I would like to get rid of the country-boundaries within the regions. Is there a way to only show the outer boundaries of each feature? Or another way to accomplish what I want?
Map of 26 world regions; each region is a feature in a new layer based on classifying countries in the original layer as belonging to one of the regions


texts: On ArcGIS online, can you create a web app that will retain auto updating functionality from a web map that auto updates using csv files?

texts: Does making your web map an application on ArGIS online cause loss of auto updating for the layers in the web app?


texts:
I have a dataset of village boundaries in ESRI format.
with the following projection information,

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG",4326]]

Then , I tried to overlay(with the help of ggplot) it on an image from google static maps using get_map in R(which is basically Google static maps).
For this village 'Chhajju Nagala',I clearly see an offset(village is actually to the lower left).

This problem is with many villages.
I want a general solution for all them.
Should I change the projection or towgs84 parameters?
What to do?


texts: QGIS not loading jpg with coordinate information

texts:
I have downloaded some jpg files that came with jpgw cordinate information but when I bring the data into qgis using the Data Source Manager it doesn't ask me for the projection and it doesn't place the file in the correct spot.

I have tried to save the file as wld, jpw and jgw with no luck.

Examples can be downloaded from https://qimagery.information.qld.gov.au/ for free - The ones I tried were 1974 for Poruma island.

I am running
QGIS version
3.10.1-A Coruña
QGIS code revision
ef24c526da

I tried the 'fixes' in What causes wrong placement of JPGs in QGIS 2.8.1? as well.


texts:
I have downloaded a hgt file from EarthExplorer but when I try to load it on QGIS I get an error "GDAL provider cannot open GDAL dataset".

I was under the impression that hgt files were supported. It this a known issue? Is there any workaround?


texts: how to reduce Geoserver WMTS rendering time?

texts:
I am publishing 100 Gb of aerial image til on GeoServer using pyramid technique. I have made GeoServer in a production environment according to GeoServer manual "Running in a production environment".
I have prepared data using gdal_retile:

gdal_retile.py -s_srs EPSG:27700 -v -r bilinear -levels 4 -ps 500 500 -co 'TILED=YES' -co 'COMPRESS=JPEG'  -targetDir OUTPUT *.jpg

Input jpg image details:

  Driver: JPEG/JPEG JFIF Files: SS0096.jpg
         SS0096.jgw Size is 4000, 4000 Coordinate System is `' Origin = (200000.000000000000000,197000.000000000000000) Pixel Size =
  (0.250000000000000,-0.250000000000000) Metadata:
  EXIF_BitsPerSample=8 8 8   EXIF_ColorSpace=65535   EXIF_Compression=1
  EXIF_DateTime=2016:09:08 07:20:28   EXIF_Orientation=1
  EXIF_PhotometricInterpretation=2   EXIF_PixelXDimension=4000
  EXIF_PixelYDimension=4000   EXIF_PlanarConfiguration=1
  EXIF_ResolutionUnit=2   EXIF_SamplesPerPixel=3   EXIF_Software=Adobe
  Photoshop CS2 Windows   EXIF_XResolution=(72)   EXIF_YResolution=(72)
  Image Structure Metadata:   COMPRESSION=JPEG   INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr Corner Coordinates: Upper Left  (
  200000.000,  197000.000) Lower Left  (  200000.000,  196000.000) Upper Right (  201000.000,  197000.000) Lower Right (  201000.000,
  196000.000) Center      (  200500.000,  196500.000) Band 1 Block=4000x1 Type=Byte, ColorInterp=Red   Overviews: 2000x2000,
  1000x1000, 500x500, 160x160   Image Structure Metadata:
      COMPRESSION=JPEG Band 2 Block=4000x1 Type=Byte, ColorInterp=Green   Overviews: 2000x2000, 1000x1000, 500x500, 160x160   Image Structure
  Metadata:
      COMPRESSION=JPEG Band 3 Block=4000x1 Type=Byte, ColorInterp=Blue   Overviews: 2000x2000, 1000x1000, 500x500, 160x160   Image Structure
  Metadata:
      COMPRESSION=JPEG

Windows Server configuration is Intel Xeon E5, 2.20GHz (2 processor) and 32gb RAM. But when I am accessing this layer using WMTS service, it is taking a lot of time (around 20-30 minutes) to load completely.

Can anyone suggest me how to reduce rendering time?


texts: Converting ms-access database into easy to use spatial database

texts:
I've recently come into ownership of an MS-Access database (2007 accdb), which I want to be able to use in a GIS. I know a little bit about GIS (I've used QGIS and ArcGIS quite a bit) but very little about databases. There is currently no spatial data in it but I do have various shapefiles which I want to link it to and ideally merge all these files/databases into one spatial database.

If I had access to ArcMAP I'd probably just set up an ESRI Geodatabase, but that's not available. From what I've read PostGIS might be something I could use but as I say I know very little about databases beyond MS-Access.

Is there a user friendly system someone could point me in the direction of? It would need to be compatible with MapInfo 7 (don't ask!) or QGIS. Any help/advice would be much appreciated.