texts: Using PyQGIS in standalone scripts without crashing

texts:
I am new to the PyQGIS API.

For developing I use my Spyder IDE with QGIS and GRASS modules. In order to make my script executable I just run my code with the QGIS python.exe in my bat file. This all works so far. However, none of the PyQGIS cookbook code snippets worked and for some reason the only thing that did the job was this example I could adapt (posted several times by user @Germán Carrillo). Also it took me quite a while to figure out that the order how you import your modules matters (PyQt4.QtGui before qgis.core breaks the code).

So for example I want to create a simple hillshade from a DEM. This works with my code from my Spyder IDE if I run it once. When I run it another time, it causes python to stop working. When I run it as a bat file, I still get the hillshade, but again I get the message 'python.exe has stopped working'. This is not clean and helpfull at all if I want to develop something more complex and let it run in an automated way. This answer didn't work for me unfortunately. Also, either way the processing module copies its folder structure (only the folders though, not its content) into my working directory together with a qgis-auth.db file. This also seems not clean at all.

Here my code with a download link provided to a small DEM sample file:

# get small DEM from https://github.com/GeospatialPython/Learn/blob/master/dem.zip

# Prepare the environment
import sys
from qgis.core import *
from PyQt4.QtGui import *
app = QApplication([], True)
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()
print "QGIS successfully Initialised"

# Prepare processing framework
# https://gis.stackexchange.com/questions/129513/how-can-i-access-processing-with-python
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins') # Folder where Processing is located
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

# Write your code here to load some layers, use processing algorithms, etc.
rasterLyr = QgsRasterLayer("dem.asc", "hillshade")
rasterLyr.isValid()

general.runalg("gdalogr:hillshade", rasterLyr, 1, False, False, 1.0, 1.0, 315.0, 45.0,"qgis_hs.tif")

# Exit applications
QgsApplication.exitQgis()
QApplication.quit()

Again, this code work in my Spyder IDE when run once, but causes python to stop working when run another time. It also works when run from a bat file, but also crashes my python.exe. Any ideas? That way I can't replace ArcPy with PyQGIS unfortunately.

The problem definitely seems to be related to how QApplication and QgsApplication is initialised (or not correctly ended?). I just tried an NDVI calculation without the messy processing library and I still get the same crashing errors.

Yes, most definitely a Qt problem that is caused by not properly closing the app. This was reported quite frequently. Unfortunately, none of the proposed solutions work for me (or I just don't know how to implement correctly).


texts: Combining proximity criteria within Location Allocation?

texts:
I am using ArcGIS 10.5.
I will explain the situation first:

There are 134 139 forest centroids that serve as 'Demand point' to source Bioenergy conversion plants in the state of Queensland Australia. There are 24 Existing Bioenergy plants 'Facilities' that 'Require' demand. Finally there are 847 potential new locations for a bioenergy 'Facility' that serve as a 'Candidate' plant.

For my research I want to know how many Forest is being allocated to the 24 existing facilities within 80km radius of each. This can be easily done by the Location Allocation under Network Analist.
Additionally I want to find which of the 847 can also serve as a facility without taking resource away from the 24 existing facilities or each other.

Therefor I want to select only those facilities out of 847 that are closest to:

sawmills and towns (points);

distribution network and rail (line).

When facilities are within each others 80 km radius I want to select the most feasible one and allocate the resource (forest) to the selected.

Ideally I want to rerun the location allocation network analysis with the 24 existing and a small selection of new facilities (out of the 847 proposed).

Any ideas on how to make a selection withing or outside the location allocation tool.


texts: Why does the peak of a mountain appear bright on Radar images?

texts:
Why does the top of a mountain appear bright on a SAR (radar) image.
I know that the peak of a mountain is subject to a layover effect, thus appearing closer on the radar image than the bottom part of the mountain. And, that the back part of a mountain is dark (shadow) because it is not reached by radar's radiations.

The only reason I see that makes the top of a mountain look brighter than the bottom part is because the slope of the mountain is narrower on slant range compared to its real dimensions on ground range. This implies that backscatter return from many objects is returned inside the same cell (more intensity). But I'm not sure if that is really the exact reason.

See Sentinel-1 SAR image of mountain (to understand what I'm talking about):


texts: Adding and saving layers from QGIS console

texts:
I have a data-set of .tif satellite images and .geoJSON labels corresponding to them.
I am trying to iterate over the satellite images,plot the geoJSON labels and save them as .jpeg images.

However,only the first image is being saved correctly.The others are either completely blank or have very little of the original .tif image on them.
I am very new to QGIS.

Here is my code

for el in [1,3,4,7,8,9]:
    iface.addRasterLayer("Path/RGB-PanSharpen_AOI_2_Vegas_img"+str(el)+".tif", "Image") #Adding Image
    iface.addVectorLayer('Path/buildings_AOI_2_Vegas_img'+str(el)+'.geojson','Label','ogr') #Adding labels
    eMenu =iface.viewMenu()
    eMenu.actions() [11].trigger() #Zooming in to layer
    iface.mapCanvas().saveAsImage("Desktop\Im"+str(el)+".jpeg",None,"jpeg") #Saving
    QgsMapLayerRegistry.instance().removeAllMapLayers() #Clearing canvas

I am using QGIS 2.14.20-Essen


texts: Finding regions in raster with rasterToPolygons (R)

texts:
I have a raster file where all pixels either are 0 (background) or 1 (see below, available at https://drive.google.com/file/d/1s2g8srXSTL6t-znF9nWbwK8f0C0K5FNk/view?usp=sharing).

I want R to find all regions of connected pixels and give them a unique ID. Then, I want to generate a shapefile via rasterToPolygons with one polygon for each region.
From what I can see, clump does the region grouping well. My code so far is:

test <- raster("C:/Users/.../sinks.tif")
regions <- clump(test)
regions_shp <- rasterToPolygons(regions, na.rm = T)

But that rasterToPolygons gives me a single polygon for every pixel and does not take into account the clump ID. What am I doing wrong?


texts: How do I calculate the pixel areas for the ee.Reducer.frequencyHistogram() calculation?

texts: I have been following this tutorial: Reducing Land Cover Data Sets in Earth Engine , but cannot calculate the pixel values into areas before exporting the tables. I am very new to Google Earth Engine and cannot find the answer elsewhere.


texts: Change the color attribute of a mulitline layer

texts:
I'm very new in world of the shapefiles, i have a code that add point in a multiline layer, but when i open with ArcGis each layer have a different color, and i want that all the layer have the same color automatically. I have this code but it isn't work.

How i can fix it?

def create_line(coords,scale):
    line=ogr.Geometry(ogr.wkbLineString)
    for coord in coords:
        coord[0]=scale[4]+(coord[0]*scale[0])
        coord[1]=scale[5]+(scale[3]*coord[1])
        line.AddPoint(coord[0],coord[1])
    return line

.....

            new_shape=str('Vias '+`key`+`z`)

            new_shape= new_shape + '.shp'

            #Crear una nueva fuente de datos y capa para los elementos detectados
            if os.path.exists(new_shape):
                driver.DeleteDataSource(new_shape)
            ds= driver.CreateDataSource(new_shape)
            print "El shapefile creado es " + new_shape

            if ds is None:
                print 'No puede crear el shapefile'
                sys.exit(1)

            #Crea una capa de elemntos
            layer=ds.CreateLayer(new_shape, geom_type=ogr.wkbMultiLineString)

            #Añade un campo id al shapefile de salida para los elementos encontrados
            fieldDefn=ogr.FieldDefn('id', ogr.OFTString)
            layer.CreateField(fieldDefn)
            z=z+1
            multiline=ogr.Geometry(ogr.wkbMultiLineString)

        line=create_line(data[key][key1],scale)

        multiline.AddGeometry(line)
        h=0
        if line:
            featureDefn=layer.GetLayerDefn()

            #Crea un nuevo rasgo
            feature= ogr.Feature(featureDefn)
            feature.SetGeometry(multiline)
            feature.SetField('id',h+1)
            feature.SetColor(QtGui.QColor('blue'))

            #Añade el rasgo de la capa de salida
            layer.CreateFeature(feature)

            h=h+1


texts: Leaflet - Can I use a different background tile provider depending on zoom level?

texts:
This may look strange, but I cannot made up my mind between providers of satellite tiles for Leaflet. ESRI sattelite looks beautiful, but cannot reach the zoom level I need. HERE is not nice on low zoom levels, but goes as deep as I want.

Can I use ESRI for low zooms and HERE for high ones?

How do I do this?

Do you think it is wise?


texts: Updating shapefile metadata to show up in ArcCatalog?

texts:
I am attempting to update the metadata for a variety of shapefiles to change the information under 'Summary' (or 'Purpose'), 'Description' (or 'Abstract'), etc. in ArcCatalog's Description tab.
I am using ArcMap 10.4 and am working in C#.
I can already update the Metadata itself and I can see that my changes have been applied to the shapefile's .shp.xml file, but this doesn't consistently change the info in ArcCatalog's Description tab. Looking at several shapefiles I have been given, it seems that there is some other setting that dictates where the ArcCatalog pulls its information from.
For example here are two different paths for two different shapefiles from where ArcCatalog gets the Description:

idinfo/descript/abstract
dataIdInfo/idAbs

The values for those elements are unique to the file, so there is nowhere else the data could be getting pulled from.
What do I need to set in my shapefiles to make ArcCatalog look at a consistent location for it's sections?


texts: How to split a county layer, but keep all surrounding counties per layer in QGIS

texts:
I have a shapefile of all US counties.
I need to split the shapefile by county, and save in a county file in the file type GeoJSON. This makes about 3k files.
That's easy.
The hard part is that for each county shapefile, I also need the surrounding counties in the shapefile.
For example, the McHenry County shapefile will have the polygon representing McHenry County, but also all the touching counties (i.e. Lake, Boone, deKalb....etc.).  The Lake County shapefile would have the touching counties (i.e. Mchenry, etc..)
The resulting shapefiles will then be exported out has GeoJSONs.
How do I capture the "surrounding counties" along with the target county, and save off as GeoJSON?


texts: HTTP 500 error while open http://localhost:8080/geoserver/web in windows 10

texts:
Im just now installed the GeoServer in my system . The installation process fine.After the installation i have started the GeoServer then i tried to use geo server web admin page but it is showing some error.in google chrome This is a proxy server. Does not respond to non-proxy requests.

And the internet explorer shows that

**HTTP 500 error

That’s odd... Microsoft Edge can’t find this page

This page can’t be displayed, because this site’s server might be under maintenance or there could be a programming error.
**

I have uninstalled the GeoServer and the i change the port 8080  to 8085 then also not working well. this time the error is changed ie, dial tcp [::1]:8085: connectex: No connection could be made because the target machine actively refused it.


texts: Calling layers from Google Maps API without success?

texts:
I'm following the examples of a free on OpenLayers, in the code below I am calling various layers of google maps api. Have checked with the code book and is correct but nothing appears on the browser.

Tested with one at a time and appeared, does anyone know what is wrong in the code?

Also I've taken the spec version of the api call to google maps and did not work. Nothing appears.

    <!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>

 <!-- Chamada a Api Google Maps especificando a versão 3.2 -->
 <script src="http://maps.google.com/maps/api/js?sensor=false&v=3.2"></script>

 <script type='text/javascript'>

 //Variável map
 var map;

 //Função a ser chamada para criar o mapa na leitura da página
 function init() {

 // A classe Map espera dois argumentos "map_element" que o ID do elemento HTLML geralmente uma DIV e opções de mapa {key:value}
 map = new OpenLayers.Map('map_element', {});

 //Função de inicialização  - google.maps.MapTypeId.HYBRID (híbrida)
 var google_hybrid = new OpenLayers.Layer.Google(
    "Google Hybrid",
    {type: google.maps.MapTypeId.HYBRID}
    );

 //Função de inicialização  - google.maps.MapTypeId.TERRAIN (topológica)
 var google_physical = new OpenLayers.Layer.Google(
    "Google Physical",
    {type: google.maps.MapTypeId.TERRAIN}
    );

 //Função de inicialização  - google.maps.MapTypeId.TERRAIN (topológica)
 var google_satellite = new OpenLayers.Layer.Google(
    "Google Satellite",
    {type: google.maps.MapTypeId.SATELLITE}
    );

 //Função de inicialização  - google.maps.MapTypeId.ROADMAP (street layer default)
 var google_streets = new OpenLayers.Layer.Google(
    "Google Streets",
    {
        );

//Adiciona as Google Layers ao mapa
map.addLayers([google_hybrid,google_physical,google_satellite,
    google_streets]);

 //Camada de controle que vai mostrar as camadas no mapa
 map.addControl(new OpenLayers.Control.LayerSwitcher());

 //Verifica se o mapa tem um ponto central e o extende a sua extensão máxima
 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: 1000px; height: 900px;'>
 </div>

</body>
</html>


texts:
I'm trying to rename a layer while converting it from a shapefile to sqlite. From looking at the ogr2ogr documents, -nln name should do the trick. I've seen this similar post and this very similar post  but neither will work for me.

ogr2ogr -f "sqlite" output.sqlite input.shp -nln "new_layer_name"


texts: ArcEngine: Drawing graphics returned from BackgroundWorker

texts:
My application is drawing graphics objects on the map(as simple as adding IElement to GraphicsContainer).
When I create graphics elements and add them to GraphicsContainer in BackgroundWorker DoWork method, ActiveView Refresh method does not show any added graphics.
Are there any triks to take into account when drawing on the map in different thread? (if at all possible)
Any help will be greatly appreciated.


texts: arcpy.mapping.Export creates some labels but not others

texts:
I've created a simple python toolbox (pyt) to export map layouts to multiple formats.  My problem is that some labels which are defined in the layout are exported, but others are not.  For example point labels export but contour labels dont.  When I export from within ArcMap all labels are exported.

Question: Is there any argument missing from my "arcpy.mapping.Export" statements that would affect labeling?

I had problems posting the whole toolbox but it can be downloaded from here: Export MXD Toolbox

Here's the code that does the exporting.

# mxd name
mxd = arcpy.mapping.MapDocument(parameters[1].valueAsText)
# workspace
env.workspace = parameters[0].valueAsText
# workspace text to use in paths
workSpace = parameters[0].valueAsText
AddMessage("The workspace is " + env.workspace)

# resolution
dpi = parameters[2].valueAsText
# filename
name = parameters[3].valueAsText
# option for formats to export
createPDF = parameters[4].valueAsText
createTIFF = parameters[5].valueAsText
createEPS = parameters[6].valueAsText
createAi = parameters[7].valueAsText

# Export PDF
if createPDF == "true":
  AddMessage("Exporting PDF")
  arcpy.mapping.ExportToPDF(map_document = mxd, out_pdf = workSpace + "\\" + name + ".pdf", resolution=dpi,
           convert_markers=True, layers_attributes="LAYERS_AND_ATTRIBUTES")
else:
  AddMessage("You chose not to create a PDF")

# Export TIFF
if createTIFF == "true":
  AddMessage("Exporting TIFF")
  arcpy.mapping.ExportToTIFF(map_document = mxd, out_tiff = workSpace + "\\" + name + ".tif", resolution=dpi, tiff_compression="LZW")
else:
  AddMessage("You chose not to create a TIFF")

# Export EPS
if createEPS == "true":
  AddMessage("Exporting EPS")
  arcpy.mapping.ExportToEPS(map_document = mxd, out_eps = workSpace + "\\" + name + ".eps", resolution=dpi, image_quality="BEST",
           picture_symbol="VECTORIZE_BITMAP", convert_markers=True)
else:
  AddMessage("You chose not to create an EPS")

# Export to Ai
if createAi == "true":
  AddMessage("Exporting Ai")
  arcpy.mapping.ExportToAI(map_document = mxd, out_ai = workSpace + "\\" + name + ".ai", resolution=dpi, picture_symbol="VECTORIZE_BITMAP",
           convert_markers=True)
else:
  AddMessage("You chose not to create an Ai")


texts:
I would want to have a simple WPS client, that would work with GeoServer. QGis (1.8.0, Lisboa) WPS plugin (v1.1.0) doesn't work*. Any client that is SIMPLE to use, user friendly, will do (the one built in GeoServer as demo is perfect, but I'd like a third party client).
It could either be a desktop client, or a web client. I tried using OpenLayers but I got a .js file that had to be incorporated to a web page and I am no programmer to know how to do it. A simple webpage with cute little buttons is what I want.

If anybody has knowledge of an existing open-source program that does support OGC standards (WPS mainly), or simple tutorials on how to set up a web client, I would be grateful to read your answers.

Thanks

Running on Win XP.

*I posted another question regarding this error.


texts: Applying SLD style to raster with two bands?

texts:
I have a raster with two bands. Both bands only have ones and zeros. I want to apply the following style:

Band1    Band2    Color
  0        0      #000000
  0        1      #ff0000
  1        0      #ffff00
  1        1      #ffffff

So I want to apply a style considering values on both bands, I searched on GeoServer documentation but I didn't find any similar example


texts:
I am currently trying to load a MPK package using ArcGIS runtime 100.2.1 for .NET (WPF). The package comes from external sources, and I do not have ArcGIS Pro/ArcMap at the moment as I am trying to do a basic evaluation for now.

Apparently, I thought it would work as I saw the documentations of ArcGIS runtime 10.2 stating the support of "MPK" files, then gradually I realized that "MPK" seems to have disappeared in the documentations of ArcGIS runtime 100 onwards, and the documentations started to mention about "MMPK".

What are the differences between "MPK" and "MMPK"? Can they be used interchangeably, or be converted without the use of ArcGIS Pro/ArcMap?


texts: Iterate through rasters in subfolders using Model Builder

texts:
I have a directory folder of precipitation data; within this folder are subfolders for years 1901-2014. In each of these subfolders are 12 .asc grid files; each file represents a month of the year (i.e. Jan through Dec).

I am trying to create a model that will iterate through each of these subfolders and perform a calculation on the 12 grid files (for example: calculate total annual precipitation). I may want to perform other calculations, but for now would like to just have the model working.

I have tried the Iterate Workspace function with the Raster Calculator, and the Iterate Rasters function with the Raster Calculator, but could not get either to work how I wanted.

I also tried Iterate Rasters with Collect Values, and Cell Statistics - but that is also not working.

Below are the models I have tried:


texts: Open-source software solutions for an Esri user

texts:
I am currently a student and have been using Esri software pretty much exclusively for my GIS needs. I occasionally use R for raster creation, especially as outputs of statistical models - which I usually bring back into Arc to work with.

I will be graduating in a year or so and may not continue to have access to ArcMap.

I am wondering what programs are recommended in the place of Arc. I currently use R for most of my raster work, but find that it can not always handle large rasters.

Of the other open source GIS software programs, is there one that has the most similarity to Arc (in capabilities and UI), or one (or some) that is easiest to learn coming from that background?

Is it common to use multiple programs for different functions, or do most people have only one or two programs they use?


texts: Can I Change .QGS File Association from One Version to Another?

texts:
I have used OSGEO4w as my installer for QGIS. I've noticed that it loads both v1.8 and v1.9. I don't have any issues with that except that I prefer to use the more stable v1.8.

I've noticed that when I open a pre-existing project file (.qgs) it defaults and opens in 1.9 (i.e. I open the project by double clicking on the project file. This is different than opening the QGIS desktop interface first and then loading the project file. The latter method works because I open v1.8 first). I checked my file association (Windows 7) for .qgs and it has 'qgis-dev' which is fine, but I'm wondering how I can get it to open in 1.8 instead of 1.9?

Thanks,
Mike


texts: Change the values of NA cells in a raster by using a geographic subset of the raster as a condition

texts:
I would like to replace the NA cells that are within a specified extent (here rectangular area defined by drawExtent in the code bellow) of a raster object. My idea was to use the functions cellsFromExtent and extract to extract the NA cells within an extent of the raster and to assign the value 250 to these cells. Finally, the NA cells of the raster that are within the extent should have the value 250 and the other NA cells that are outside the extent would keep the value NA. So, How can I change the values of NA cells in a raster object by using a geographic subset of the raster as a condition ? Here is the beginning of a code:

r <- raster(ncol=10, nrow=10)
values(r) <- sample(1:8,ncell(r),replace=T)
r[c(5)] <- NA
r[c(20)] <- NA
r[c(43)] <- NA
plot(r)

e <- drawExtent()
z <- extract(r, cellsFromExtent(r, e))

From z, how can I assign the value 250 to the NA cells in the raster (I tested r[is.na(z)] <- 250 but this doesn't work) ?


texts: OSM Map quality deteriorates when changing CRS - QGIS

texts:
I have a problem with basemaps for my drawings in QGIS.

By basemaps I mean OpenStreetView, Google Maps etc. I'm aware that most of them (if not all) are natively set in EPSG:3857, and they look just fine in that CRS. However, when I'm required to use a different coordinate system in my project - for example EPSG:32630 - they understandably get a little distorted.

Are there any basemaps that would accommodate to the chosen CRS without the loss of quality?

Below I enclosed screenshots of OSM in 3857 and 32630:

3857:

32630:


texts:
I'm trying to extract elevation data from a DEM file but I can't make it work.

I want to extract contour lines from a given area, and label those lines with the proper elevation value. The problem I'm facing right now is that I managed to extract the contours with the contour algorithm, but the elevation values are wrong.

The steps I follow are:

Download the DEM file (.asc) and add it to my QGIS project.
Run the contour algorithm on the DEM layer. Since I'm working with a 2x2 DEM file, I set the interval between lines to 2 (not sure if that's correct) and I also set the atribution name to "ELEV" (default)
Label the output vector layer with the contours to show the "ELEV" values on each contour line.

Then, I used the "identify features" tool on the DEM layer, and clicked a random point that I knew the real elevation to check if the DEM had bad data, but it was correct. So no issues with my DEM.

On the other hand, the "ELEV" label value was not correct. I was getting an error of like ~80 units (the label value was 90 and the real value was ~170)

So there's something wrong with my procedure, but I don't know how to fix it.

Additional info:

SRC  EPSG:25831 - ETRS89 / UTM zone 31N - Proyectado
Tipo de datos (data type)    Float64 - Número de coma flotante de 64
bits
Banda 1   STATISTICS_MAXIMUM=322.522 STATISTICS_MEAN=171.58152588021
STATISTICS_MINIMUM=87.393 STATISTICS_STDDEV=43.912078992413
STATISTICS_VALID_PERCENT=100

I feel like the "Banda 1" data is the one that I need, and the one that shows when I click on a random point using the "identify features" tool.


texts: What's the syntax for an AttributeFilter in FME?

texts:
I want to use AttributeFilter and I need to filter data by ATTRIBUTE1 with possible values of ATTRIBUTE2. What is the syntax for it?

I tried:

attribute2
&attribute2
$attribute2
'&attribute2'

etc. but wasn't successful.


texts: Accessing Noder functionality through C API of GEOS?

texts:
All this stuff looks rather useful:
http://geos.osgeo.org/doxygen/classgeos_1_1noding_1_1Noder.html

However, from the GEOS README file:

The C++ interface should not be used directly; the geos project
views it as a bug for another program to use the C++ interface or even
to directly link against the C++ library.

Is there a C interface to the Noding functionality that I'm missing?  I can't find it in geos_c.h.


texts: readFeatures is returning the wrong (or no) projection

texts:
I'm trying to change my map from using JSON to JSON-P for the source of data from my GeoServer. Unfortunately, what should be a vector layer covering the entire North Atlantic displays as a point in the Gulf of Guinea--so a classic projection error. But I can't see what I'm missing.

TL;DR: If I actually request the data from Geoserver using SRSNAME=EPSG:900913 I don't have a problem, but I can't figure out why I didn't need to do that originally!

The original vector layer is created by:

var vectorLayer = new ol.layer.Vector({
    source: new ol.source.Vector({
      url: '//'+$(location).attr('host')+'/geoserver/SAHFOS/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=SAHFOS:CPRStdAreasIntRes&maxFeatures=50&outputFormat=application/json',
      format: new ol.format.GeoJSON()
    })
});

and works fine. The 41 features fit exactly to the coastlines of a world map.
Trying to do it with JSON-P:

    var stdAreas = new ol.layer.Vector({
        // Source retrieving WFS data in GeoJSON format using JSONP technique
        source: new ol.source.Vector({
            loader: function(extent) {
                $.ajax('http://domain/geoserver/SAHFOS/ows', {
                    type: 'GET',
                    data: {
                        service: 'WFS',
                        version: '1.0.0',
                        request: 'GetFeature',
                        typeName: 'SAHFOS:CPRStdAreasIntRes',
                        maxFeatures: '50',
                        outputFormat: 'text/javascript'
                        },
                    dataType: 'jsonp',
                    jsonpCallback:'callback:loadFeatures',
                    jsonp: 'format_options'
                })
            }
        })
    });

    var loadFeatures = function(response) {
        var source = stdAreas.getSource();
        source.addFeatures(new ol.format.GeoJSON().readFeatures(response, {dataProjection: 'EPSG:4326'}));
        var fmt=new ol.format.GeoJSON().readProjection(response)
    };

fmt is EPSG:4326,  so that dataProjection shouldn't be necessary, but nevertheless, the layer still ends up being a point in the Gulf of Guinea.


texts:
I created a web map with the qgis2web plugin. The map simply contains 2 polygon and 1 point layer. Since I just wanted a static map with these 3 layer I don´t use any basemap.

The final map should be placed in a frame on website and thats where things started to get dicey. The map doesn´t scale correct when the window size of the browser is changed or on mobile devices.

I assume since the map is in a “sub”-frame of the website the

width=device-width

from meta doesn´t kick in since the frame always stays in the same proportion to the Website.

I´ve looked around for a solution and found this:

EventListener

I would work fine, but I would prefer a more dynamical solution rather than defining the possible widths for all zoom levels..

Note that the map will always be on a fixed zoom level (all zoom controls are disabled). The zoom level should be the maxmium possible with in the frame, so that all Layer-content is visible.


texts: Twitter2QGIS plugin error "EOF occurred in violation of protocol"

texts:
I installed twitter2qgis plugin and also tweepy through pip manually and then, I put access token etc value from my developer twitter account and when I apply I get this error:

And this is written in it:

SSLError: [Errno 8] _ssl.c:504: EOF occurred in violation of protocol Traceback (most recent call last): File "C:/Users/HAMMAD AHMAD/.qgis2/python/plugins\geotweet\geotweet.py", line 321, in run stream.filter(track=keywords) File "C:\Users\HAMMAD AHMAD\AppData\Roaming\Python\Python27\site-packages\tweepy\streaming.py", line 445, in filter self._start(async) File "C:\Users\HAMMAD AHMAD\AppData\Roaming\Python\Python27\site-packages\tweepy\streaming.py", line 361, in _start self._run() File "C:\Users\HAMMAD AHMAD\AppData\Roaming\Python\Python27\site-packages\tweepy\streaming.py", line 294, in _run raise exception SSLError: [Errno 8] _ssl.c:504: EOF occurred in violation of protocol
Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)] QGIS version: 2.18.1 Las Palmas, 3e73840
Python Path: C:/PROGRA~1/QGISLA~1/apps/qgis/./python/plugins\processing C:/PROGRA~1/QGISLA~1/apps/qgis/./python C:/Users/HAMMAD AHMAD/.qgis2/python C:/Users/HAMMAD AHMAD/.qgis2/python/plugins C:/PROGRA~1/QGISLA~1/apps/qgis/./python/plugins C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg C:\Python27\ArcGIS10.3\python.exe C:\PROGRA~1\QGISLA~1\bin\python27.zip C:\PROGRA~1\QGISLA~1\apps\Python27\DLLs C:\PROGRA~1\QGISLA~1\apps\Python27\lib C:\PROGRA~1\QGISLA~1\apps\Python27\lib\plat-win C:\PROGRA~1\QGISLA~1\apps\Python27\lib\lib-tk C:\PROGRA~1\QGISLA~1\bin C:\Users\HAMMAD AHMAD\AppData\Roaming\Python\Python27\site-packages C:\PROGRA~1\QGISLA~1\apps\Python27 C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\PIL C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\win32 C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\win32\lib C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\Pythonwin C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg C:\PROGRA~1\QGISLA~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg C:/Users/HAMMAD AHMAD/.qgis2//python


texts: Getting only the exterior boundary of a 'donut-shaped' GeoSeries object

texts:
GeoSeries.boundary() finds the entire object's perimeter, so I supposed that wouldn't work.
I made my GeoSeries donut-like shape using unary_union of smaller GeoSeries objects.


texts: automatically zoom map to maximum extent of wms layer

texts:
I try to automatically zoom map to maximum extent of layer after adding
layer to map.
I tried to code like this.

var layer = OpenLayers.Layer.WMS(.......);

map.addLayer(layer);

map.zoomToExtent(layer.getDataExtent());

But it doesnt' work. I also changed like this

map.zoomToExtent(layer.getMaxExtent());

It still doesn't work. What can I do to solve this issue.


texts: how to build topology with 3d linestring in pgrouting

texts:
I have some 3D linestrings in a postgis database, and I want to use pgrouting for calculating the shortest path with Dijkstra.

When I execute select pgr_createTopology('route2', 0.000001, 'geo3d', 'id'); ,
there is the error:

  Geometry has Z dimension but column does not

How can I create topology with 3d linestring by pgrouting?


texts: Running Point To Raster has output coordinate system different from input?

texts:
I need to transform a DEM from a Point-Feature-Class (Shape) to raster (GRID). For that I'm using the tool 'point to raster' and it works perfectly.

The only problem is that it seems to change the coordinate system, so that my resulting GRID has another system than my Shape.

The Shape has DHDN_3_Degree_Gauss_Zone_3 (which I want to keep), but the resulting GRID has Deutsches_Hauptdreiecksnetz_Transverse_Mercator.

Can anyone tell me why that is and how I can prevent it?

Or is there any way to transform the Transverse_Mercator-system back to Gauss?

I fill in the tool dialog as follows:
input features: my point FC which has the DHDN_3_Degree_Gauss_Zone_3 CS
value field: my z-values (height of the DEM)
output raster dataset: destination folder
cell assignment: mean
priority: none
cellsize: 1
and set the output coordinate system to 'same as input' in the environments settings.


texts: Line feature reshape tool reshapes unselected lines

texts:
I have QGIS 2.18.18 installed in my PC. I need to reshape large number of roads/tracks to match with Imagery. QGIS Reshape tool reshapes roads that are not even selected. Any suggestion?

Below are the screenshots for your reference

After making selection I modified the road using reshape tool

As a result, tool reshaped unselected adjacent road


texts:
I am using QGIS, version 3.10.14 and the semi-classification plugin and want to create a training input, to classify different land-use-types. I followed the manual (https://readthedocs.org/projects/semiautomaticclassificationmanual-de/downloads/pdf/latest/, page 163).
When I define a training input, it does not appear in the line within the SCP Dock.
When I try to save my classes, I get the following error:

Information [3]: Select a SCP training input; input is not loaded

This question has already been asked, but there was no answer. How can I define the training input?


texts: Exit code 139 while using ReadAsArray in OS X

texts:
I'm trying to load TIFF image of DEM as a numpy array using method ReadAsArray() on Mac OS:

dataset = gdal.Open('DEM.tif')
band = dataset.GetRasterBand(1)
array = band.ReadAsArray()

This code crashes Python with statement

Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

According to this question https://stackoverflow.com/questions/23704040/exit-code-139-when-perfroming-image-subtraction, I get error 'Invalid memory access'. But in this manual I found SIGSEGV as signal for error 'Invalid memory reference'. What does it mean and how to fix it? For Windows this code works fine.

Method ReadRaster() also works fine, but I need to get a bumpy array of my data.


texts: Rendering a custom bounding box to an image file using python and QGIS?

texts:
Edit: the problem was merely a typo in the coordinates (see comment on definition of p). Otherwise the code below should work for other situations, too, if adapted properly. The weird pixelated/grey/colorful images arise, when the rendered coordinates do not lie on the layer extend.

I have tried to adopt some code from the PyQGISDeveloperCookbook (p.37) to render a QGIS raster layer to a png image file (see bottom for layer's metadata). One thing that I have tried to change is to modify the rendered extent to a custom bounding box bbox centered in a given point p:

    from PyQt4.QtGui import QImage, QPainter
    from PyQt4.QtCore import QSize
    import os

    # raster layer in CRS UTM 33N (epsg 25833)
    rlayer = iface.activeLayer()
    # radius for image in meters
    radius = 50
    # center point
    p = [33401477,5826362] # <- should be 334014.77 for x-coordinate
    # bounding box
    bbox = QgsRectangle(p[0]-radius, p[1]-radius, p[0]+radius, p[1]+radius)
    # image
    img = QImage(QSize(200,200), QImage.Format_ARGB32_Premultiplied)
    ## create painter
    painter = QPainter()
    painter.begin(img)
    render = QgsMapRenderer()
    ## set layer set
    lst = [rlayer.id()]
    render.setLayerSet(lst)
    ## set extent
    render.setExtent(bbox)
    ## set output size
    render.setOutputSize(img.size(), img.logicalDpiX())
    ## do the rendering
    render.render(painter)
    painter.end()
    ## save image
    img.save("D:\\Daten\\QGIStmp\\bbox1.png", "png")

The resulting image does not really look the way I expected (though interesting :-)):

I tried hard (generating more interesting images), but I couldn't figure out what the problem could be.

Querying Layer Metadata:

>>> rlayer.metadata()

Response:

    u'
Treiber
\n
GDAL provider
\nECW
ERDAS Compressed Wavelets (SDK 5.0)
Datensatzbeschreibung
\n
D:/Daten/QGIStmp/Nordost/dop20_400_5826.ecw
\n

\nCOLORSPACE=RGB
\n
\nCOMPRESSION_RATE_TARGET=20
\n
\nVERSION=2
\n

Kanal 1
\n
Kanal 2
\n
Kanal 3
\n
Dimensionen
\n
X: 10000 Y: 10000 Kan\xe4le: 3
\n
X : 5000,Y 5000

X : 2500,Y 2500

X : 1250,Y 1250

X : 625,Y 625

X : 312,Y 312

X : 156,Y 156

Ursprung
\n
400000,5.828e+06
\n
Pixelgr\xf6\xdfe
\n
0.2,-0.2
\n
Leerwert
\n
*Leerwert nicht gesetzt*
\n
\n
Datentyp
\n
Byte - Acht Bit vorzeichenlose Ganzzahl
\n
Pyramiden\xfcbersichten
\n

R\xe4umliches Bezugssystem des Layers
\n
+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
\n
Layerausdehnung (in ursp\xfcnglicher Projektion des Layers)
\n
400000.0000000000000000,5826000.0000000000000000 : 402000.0000000000000000,5828000.0000000000000000
\n
\nKanal
\n
Kanal 1
\n
Kanal Nr
\n
\n1
\n
Keine Statistik
\n
\nNoch keine Statistik gesammelt
\n
\nKanal
\n
Kanal 2
\n
Kanal Nr
\n
\n2
\n
Keine Statistik
\n
\nNoch keine Statistik gesammelt
\n
\nKanal
\n
Kanal 3
\n
Kanal Nr
\n
\n3
\n
Keine Statistik
\n
\nNoch keine Statistik gesammelt
\n'


texts: Is there a way to wholesale rewrite all the metadata from gdal.GetProjection() back into a raster?

texts:
I am using gdal (Python) to get metadata out of a geo tiff using GetProjection.  I do some manipulation on the raster in openCv, and i would like to then re-merge the altered raster data with the meta data string that is returned by GetProjection.  I have been able to individually set parameters with specific methods from osr, but it requires me to parse the long string I get from GetProjection and then find the corresponding setter from osr.  Is there some way to do this all in one fell swoop?

here is my reader function --

def read_ortho4align(filename):
    # open ortho in gdal to get affine transform from pixels to geo space
    imgeo = gdal.Open(filename)
    xfrm = imgeo.GetGeoTransform()
    geoinfo = imgeo.GetProjection()
    print("xfrm for image = ", xfrm)
    imgeo = None #free-up memory
    # open ortho (again) with opencv for image content manip
    img = cv2.imread(filename, cv2.IMREAD_LOAD_GDAL) #I think this is where big file reading is failing
    imgGray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    edgemask = make_edge_mask(imgGray)
    return img, imgGray, edgemask, xfrm, geoinfo

And here is where I re-constitute my geotiff (snippet)

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, bands, gdal.GDT_Byte, options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF',])
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(32615)  #TODO - fix hardcoded
outRaster.SetProjection(outRasterSRS.ExportToWkt())


texts:
Can I import a KML file with styles to QGIS 3? (As of now version 3.4.1)

I have a .kml file with cycling routes that have a different color for the time between sunset and sunrise.

When I import the .kml to QGIS all styles are overridden.

Is there a way to display or even edit the preset KML features in QGIS?

A trivial sample file: https://gist.github.com/mkpaa/a6a6f433366019eabdc2811ecafba2af


texts: OpenLayers has wrong overlay positioning while zooming out

texts:
I'm quite new to OpenLayers. Right now I'm using it with OSM for a project and I'm facing what seems an issue to me: when I set an Overlay while zoomed in to max resolution, then, when I zoom out, the marker looks like it have moved to a different location. If I zoom in again everything is correct.

I this a bug in OpenLayers?

I'm using projection EPSG:3857 all the time and the marker positioning is set to "bottom-center".

Here is an image of the issue:

From left to right:
- 1. what I set
- 2. what I see when zooming out
- 3. What should be

Any Idea?


texts: Why is the highlighting on these features obscured?

texts:
I have a function that highlights features when the user enters search text, and that text is found in the attributes of the given features.

The code basically works; the problem is in the display. I want the highlighting to appear clear and unbroken on top of the map, but it ends ends up partly obscured by the other layers of the map. I tried using the .bringToFront method, to no avail.

How can I get the highlighting to appear strong and clear on top of the map?

Here's the code:

document.getElementById('searchText').value = "";
var searchText = document.getElementById('searchText');
var geoJsonLayer = L.geoJson('').addTo(map);

function goBtn() {
        map.removeLayer(foundFeatures);
        Zones.find().layers('0').text(searchText.value)
                .run(function(error, featureCollection, response){
                    foundFeatures = L.geoJson(featureCollection, {
                        style: function(){
                          return {
                            color: 'black',
                          weight: 5
                          };
                        }
                    }).addTo(map);
          foundFeatures.bringToFront();
          map.fitBounds(foundFeatures.getBounds());
          });
    }

Here's the result:

You can see that the pieces of the feature boundaries that do NOT border another feature actually do appear unobscured, but the boundaries that are shared with another feature are greyed out.


texts:
Shapely with IronPython 2.7
Will Shapely run with IronPython 2.7?


texts: Is it convention to index a raster band from 1?

texts:
Both GDAL and QGIS's raster calculator index raster bands starting at 1.  This is different than the usual 0 indexed convention that one finds in programming.  And despite Googling, I'm unable to find if band indexing from 1 is a standard convention, or a byproduct of GDAL's convention.  Does anyone know if it is standard convention to refer to a raster's band as 1 based indexing?

For context, I'm developing an API for raster processing and I need to allow the API to express which bands to operate over.  My intuition would be to index them from 0, but I wonder since these industry standard tools do otherwise.


texts: Applying style from QGIS to OpenLayers VectorLayer?

texts:
I have a vectorlayer exported from QGIS and added to an openlayers map. I now want to style this layer according to the style I use in QGIS. It consists of numerous line features and they should be styled according to an attribute. Since I want to apply a lot of different styles it would be a ton of work to write those styles manually.

Is there a better workflow than manually adding the styles to the layer?

I thought about applying a .sld to the layer, but openlayers does not seem to support this. I know it is possible to add the .sld to a Layer in Geoserver, but Geoserver is not an option for me here.

I use Openlayers version 5.3.


texts: Custom data source for GeoServer raster image

texts:
I'm working on a project where a custom binary data file contains latitude, longitude, and temperature data.  I'm interested in creating a DataStoreFactory plugin to enable GeoServer to load this information, and then use its own layers capabilities to create a coverage raster image.

I've spent some time looking through tutorials and the GeoTools site, and fighting with some dead links.  My progress so far has been to create an eclipse project with the following class:

package org.geotools.giops.plugin;

import java.io.File;
import java.io.IOException;
import java.io.Serializable;
import java.util.HashMap;
import java.util.Map;

import org.geotools.data.AbstractDataStoreFactory;
import org.geotools.data.DataStore;
import org.geotools.data.DataStoreFinder;

public class MyCompanyDataStoreFactorySpi extends AbstractDataStoreFactory {

    private static DataStore dataStore_ = null;

    private synchronized static DataStore getDataStore() throws IOException
    {
        if ( dataStore_ == null )
        {
            File file = new File("example.shp");
            Map map = new HashMap();
            map.put( "url", file.toURL() );

            dataStore_ = DataStoreFinder.getDataStore( map );
        }

        return dataStore_;

    }

    public DataStore createDataStore(Map<String, Serializable> arg0)
            throws IOException {
        System.out.println("createDataStore: args:");
        for(String key : arg0.keySet())
            System.out.print("key=" + arg0.get(key).toString() + ", ");
        System.out.println("");

        return getDataStore();
    }

    public DataStore createNewDataStore(Map<String, Serializable> arg0)
            throws IOException {
        System.out.println("createNewDataStore - weird");
        return getDataStore();
    }

    public String getDescription() {
        return "MyCompany Formatted Data Store";
    }

    public Param[] getParametersInfo() {
        // let's not return null at this point in the game
        Param p[] = { new Param("doit"), new Param("hmm") };
        return p;
    }

}

I've also added a META-INF/services/MyCompanyDataStoreFactorySpi file:

org.geotools.giops.plugin.GIOPsDataStoreFactorySpi

My expectation is that I can jar this up together with the META-INF folder, and then drop it into the GeoServer's webapps/WEB-INF/lib folder and have it pick it up as a new type of data store.

Currently, I don't see any indication in the logs that the jar is being noticed, or looked at.  Perhaps I need some of the UI modifying code?  I don't really know.  While my dataSource is a shapefile at this point, the end goal is to load lat/long/temp data and render a raster image.

I've noticed that DataSource GeoTools classes only reference feature data.  I'm worried I'm going in the wrong direction here.

Do you have any tips on how I should be proceeding with getting my custom datastore visible in the GeoServer's admin pages?


texts: Geoserver heatmap in meters rather than pixels

texts:
Is it possible to get the radius' density kernel of the heatmap in meters rather than in pixels?
http://suite.opengeo.org/opengeo-docs/cartography/rt/heatmap.html

If not, is there other alternatives to draw heatmaps efficiently from my postgresql db and send it to leaflet? (I tried heatmap.js but it is quite slow)


texts:
In the following CSS styling:

* {
   mark: url(smileyface.png);
   mark-mime: "image/png";
   mark-rotation: 45;
}

How can I make the URL value dynamic (to be read from a layer's field image_url)?

I've tried the following:

* {
   mark: url([image_url]);
   mark-mime: "image/png";
   mark-rotation: [rotation_angle];
}

But didn't work.
The server console says:

WARN [renderer.style] - Unable to read image at file:/Applications/GeoServer.app/Contents/Java/data_dir/workspaces/tst/styles/%5Bdescription%5D : Can't get input stream from URL!

The auto-generated SLD is:

<sld:PointSymbolizer>
            <sld:Graphic>
              <sld:ExternalGraphic>
                <sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="[image_url]"/>
                <sld:Format>image/png</sld:Format>
              </sld:ExternalGraphic>
              <sld:Rotation>
                <ogc:PropertyName>rotation_angle</ogc:PropertyName>
              </sld:Rotation>
            </sld:Graphic>
</sld:PointSymbolizer>

Although the same works for mark-rotation.


texts: QGIS Extract Nodes with M-Values for Linear Referencing

texts:
I have a MultiLineStringZM layer in a sqlite database, and I'm trying to visualize the measures or m-values at the vertexes. I've tried looking up information on how to do this in QGIS, and about all I've been able to gather is that this is not possible directly from the linestring layer and that the points need to be extracted to a separate layer.

I've used Vector --> Geometry Tools --> Extract nodes to create a multipoint layer representing the vertexes of my multlinestring layer, but the process loses the m-values of the vertexes. I need the m-values preserved by either saving the m-value as an attribute of the point, or something else?

Internally we have a command line tool that converts linestrings to a point shapefile with the m-values stored as an attribute on each point, and I've used that to verify that there are m-values assigned to the vertexes, and I could use that if I have to, but if possible it would be nice if this could be done directly inside of QGIS.

EDIT - Repeating what I've said above, but re-emphasizing the fact that we do have a command line tool that can achieve the results that I'm looking for that uses the GDAL libraries, so a solution showing just a partial answer in PyQGIS is not the answer I'm looking for. I'm looking for a built in tool, a plugin ready made for QGIS, or a complete script that can extract (not create/generate) and visualize m-values from a MultiLineStringZM or LineStringZM geometry.


texts: Footprints for Sentinel-1 GRD imagery as GeoJSON or shapefiles

texts:
I made this script to create a visualization of the images available for an area of interest:

var north_demosite = ee.Geometry.Rectangle([13.18,58.53,16.99,61.33])
Map.addLayer(north_demosite, {color: 'FF0000'}, 'North Demosite');

// Filter the collection for the VV product from the descending track
var collectionVV = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterBounds(north_demosite)
    .select(['VV'])
    .filterDate('2017-03-15', '2017-09-15');

    //.filterDate('2018-03-15','2018-09-15');
print("printing VV collection")
print(collectionVV.size())
print(collectionVV);

// Filter the collection for the VH product from the descending track
var collectionVH = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.eq('instrumentMode', 'IW'))
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
    .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))
    .filterBounds(north_demosite)
    .select(['VH'])
    .filterDate('2017-03-15', '2017-09-15');

    //.filterDate('2018-03-15','2014-09-15');

print("printing VH collection")
print(collectionVH.size())
print(collectionVH);

function addImage(image) { // display each image in collection

  var id = image.id
  var image = ee.Image(image.id);
  var array = id.split("/");
  var date = array[2];
  var date = date.slice(0, 8);
  print(image)

  Map.addLayer(image)
}

collectionVH.evaluate(function(collectionVH) {  // use map on client-side
  collectionVH.features.map(addImage, {min: -25, max: 5})

})

collectionVV.evaluate(function(collectionVV) {  // use map on client-side
  collectionVV.features.map(addImage, {min: -25, max: 5})

})

It works fine! But I would need the information of this query in a GeoJSON or shapefile format. Is this possible? I haven't seen anything in the documentation of Google Earth Engine


texts: Python: Point in polygon, boundary and vertex check (ray casting)

texts:
So I have written some code to check whether a point [(x,y)] lies within a polygon [(x,y), (x,y), (x,y), (x,y)] but the code fails if the point lies on a boundary or a vertex, I need this to be classed as lying within the polygon. This is the code:

def areasign(poly):
    sumarea = 0
    for i in range(0, len(poly)-1):
        xi = poly[i][0]
        yi = poly[i][1]
        xj = poly[i+1][0]
        yj = poly[i+1][1]
        sumarea = sumarea + ((xi*yj)-(yi*xj))

    if sumarea == 0: return 0
    if sumarea < 0:  return -1
    if sumarea > 0:  return 1

def lineintersection(xy, ij):
    xyi = []
    xyj = []
    ijx = []
    ijy = []
    xyi.extend(xy)
    xyj.extend(xy)
    xyi.append(ij[0])
    xyj.append(ij[1])
    xyi.append(xy[0])
    xyj.append(xy[0])
    ijx.extend(ij)
    ijy.extend(ij)
    ijx.append(xy[0])
    ijy.append(xy[1])
    ijx.append(ij[0])
    ijy.append(ij[0])
    a = areasign(xyi)
    b = areasign(xyj)
    c = areasign(ijx)
    d = areasign(ijy)
    if (a!= b) and (c!=d):
        return True
    else:
        return False

def openpointfile(fname):
    p = []
    f = open(fname)
    for line in f:
        line = line.replace('POINT(', '')
        line = line.replace(')', '')
        vals = line.split(' ')
        res = map(float, vals)
        ctuple = tuple(res)
        p.append(ctuple)
        print ctuple
    f.close()
    return p

points =         openpointfile('C:/Users/Dan/Documents/Informatics/Data/Point_On_Line.txt')

def openpolygonfile(fname):
    p = []
    f = open(fname)
    for line in f:
        line = line.replace('POLYGON((', '')
        line = line.replace('))', '')
        s = line.split(',')
        poly = []
        for coord in s:
            vals = coord.split(' ')
            res = map(float, vals)
            ctuple = tuple(res)
            poly.append(ctuple)
        p.append(poly)
    f.close()
    return p

polys = openpolygonfile('C:/Users/Dan/Documents/Informatics/Data/Polygon_On_Line.txt')

def pointinpoly(points, polys):
    infinity = tuple([100000000.0, 100000010.0])
    sum = 0
    testline = (points, infinity)
    for i in range(0, len(poly)-1):
        start = poly[i]
        end = poly[i+1]
        line = (start, end)
        l = lineintersection(line, testline)
        if l == True:
            sum = sum + 1
    if sum % 2 == 0:
        return  False
    else:
        return True

pointcount = 0
for point in points:
    pointcount = pointcount + 1
    polycount = 0
    for poly in polys:
        polycount = polycount + 1
        l = pointinpoly(point, poly)
        res = 'Point ' + str(pointcount) + ' lies within Polygon ' + str(polycount) + ' = ' + str(l)
        print l
        print (res)

Like i said before this doesnt work if the point lies on a boundary or a vertex, i have this function that i found on here that seems to work but i can't find a way to get it my files into the correct format for the function and to get it to iterate oer a list of polygons.

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

Any help on how to build that function into my existing code or any help on a new function that could be used to check if the point lies on the boundary of a polygon or on a vertex


texts: Finding Relative Orbit number from Sentinel 1 image

texts:
python
after = ee.Image(IC_post.sort('system:time_start').limit(1))
properties = after.getInfo()
properties['relativeOrbitNumber_start']

I am trying to find the relative orbit of a Sentinel 1 image from its metadata to filter an image collection to only that geometry.

I think that .getInfo() should return a dictionary of metadata and that the key ['relativeOrbitNumber_start'] should return the correct relative orbit but it just throws an error.

What am I messing up?

]2


texts: ArcMap label font size vs callout font size

texts:
Within ArcMap 10.1:

When converting labels to annotations I have specified a font size of 14. I have also created a number of callout text annotations within the map and specified exactly the same font, size and properties as the converted label annotations.

In the data view everything appears to be displaying at the same size. However when I go into layout view I find that despite both entities having a font size of 14 they are displaying at noticeably different sizes!

How can I ensure that all of my annotations and callouts display at the same size, without having to experiment with different font sizes?

Takk,

Andy


texts: Min. weighted average distance from multiple points

texts:
I am a complete beginner in mapping-related topics. Still, I´ve been wondering if anyone know this problem. It might be well known, I just don´t know its name.

Initial situation:

nx : x geographical locations

Wx: weight per location

Dx: distance between the location nx and the optimal location O

Problem:

Where is the geographical position O, where the Σ (Dx*Wx) is at a minimum?

Numerical example:

n1: Rheinfelden; Wx: 20%

n2: Zurich; Wx: 40%

n3: Berne; Wx: 30%

n4: Lucerne; Wx: 10%

Result: somewhere around Aarau in Switzerland

Does anyone know the common name for this spatial problem?

I have found a tool that can do this on a straight line ("as the bird flies", http://www.geomidpoint.com/) but I would like to be able to calculate this on a street routing level. Is anyone aware of such a tool?


texts: problem merging polygons based on a column value - using Spatialite

texts:
As I am unfamiliar with any tools I decided to try my luck with PostGis.
I imported my data.
and executed the query :

create merged as
select delomrade_code, st_union(geometry)
from grunnkrets
group by delomrade_code

it seems to be doing a much better job on my data.
first it executed in minutes not in hours.
and the result set has only 6 MULTIPOLYGONS out of 1546 as oppose to over 200 I got with the Spatialite.

I'm OK with this result and will try to dig deeper and understand why I couldn't get these 6 areas as I wanted.

I'm new to GIS.
I have seen relevant and possibly duplicate questions but can't solve this problem.

I'm using Spatialite.
I have a table with the basic statistical area units of norway (Grunnkrets) - 13990 polygons.
I added a column to the table named "delomrade_code" with the code of groups of these smaller polygons - 2224 unique values.

I'm trying to convert the these 13990 polygons to 2224 bigger polygons containing these smaller polygons.

first I tried running the query :

select delomrade_code, GUnion(Geometry) as geometry
from Grunnkretser
group by delomrade_code

The result geometry column has mostly polygons (success) but a few hundred lines contain MULTIPOLYGONS.
I need the result to be polygons only.
what do I need to do in order to try to avoid the multipolygons in the result ?

I suspect there is some minor difference in LatLng of point in adjoining polygons and therefore GUnion fails to merge these polygons.
I then tried running :

select delomrade_code, GUnion(Buffer(Geometry, 0.02)) as resulting_polygon
from grunnkretser
group by delomrade_code

results were better, more of the 2224 rows in the result were polygons but still not all.

The documentation of the Buffer function says :

  return a geometric object defined by buffering a distance d around
  geom, where dist is in the distance units for the Spatial Reference of
  geom

How do I find my :

  ...distance units for the Spatial Reference of
  geom


texts:
Using ArcMap 10.0

Trying to export a photo that has been taken with a device which is stored inside of a Feature class table. I have tried to export the data using "Right Click>Data>Export Data" This only exports the other data inside the table, such as the date taken and another attribute.

I am trying to export this layer into QGIS so that the raster layer is still available. So you can see the point, click on it and the photo will show, all in QGIS. Unsure on what needs to be done or if this is even possible.

The raster column is removed when exporting it to .dbf or .shp file extension. Something I am missing here?


texts: GeoPandas Intersect function gives different result to Shapely

texts:
Using this example dataset of country geometries, I seem to get a different result on whether or not the USA and Mexico intersect when I do the calculation
with GeoPandas as opposed to directly with Fiona and Shapely (upon which Geopandas depends). I am reasonably experienced with Fiona and Shapely, but quite new to GeoPandas - does anyone know why this might be?

Consider the following code snippet:

from geopandas import read_file

# load the shapefile of countries
world = read_file('TM_WORLD_BORDERS-0.3.shp')

# extract usa and mexico
usa = world[ (world.ISO3 == 'USA') ]
mexico = world[ (world.ISO3 == 'MEX') ]

# test intersection
print(usa.intersects(mexico))

This gives:

119    False
208    False
dtype: bool

However, if I do the same thing directly with Fiona and Shapely:

from fiona import open
from shapely.geometry import shape

# load the shapefile of countries
with open('TM_WORLD_BORDERS-0.3.shp') as countries:

    # extract usa and mexico
    for feat in countries:
        if feat['properties']['ISO3'] == 'USA':
            usa = shape(feat['geometry'])
        elif feat['properties']['ISO3'] == 'MEX':
            mexico = shape(feat['geometry'])

# test intersection
print(usa.intersects(mexico))

Then I get:

True

Clearly, I would expect both code snippets to evaluate to True. I had wondered if Geopandas handled MultiPolygons differently (as USA is a MultiPolygon, whereas Mexico is a Polygon), but I can't seem to get a positive intersection by isolating the individual polygon either.


texts: Possible data corruption when copying geodatabases to and from multiple hard drives

texts:
I hope to manage an ArcGIS Pro project and geodatabase between my personal laptop, external hard drive, and a shared office computer.

For the current project, I am accessing files from one geodatabase as well as shapefiles in a separate folder. The below graphic displays my file structure on my external hard drive.

If I copy and paste the entire folder (including the shapefile folder and geodatabase) from the external hard drive to the desktop of either computer, will the geodatabase be corrupted?


texts: L.esri.featureLayer shows up and sometimes not

texts:
I'm new to the leaflet/esri subject so please forgive me if it's just a simple question.

I wrote my first webmap-app. It uses different sources. The main aim is to show two basemaps hosted by MapBox with an AGS layer coming from a different source.

Mapbox does not make any trouble. The problem is the AGS layer.
Its address is:

http://geoportal1.stadt-koeln.de/ArcGIS/rest/services/Stadtplanthemen/MapServer/9

My script shows four different behaviours. I tested it with two different AG Servers:

Test with http://geoportal1.stadt-koeln.de/ArcGIS/rest/services/Stadtplanthemen/MapServer/9
a: It shows up in IE10 correctly
b: Not showing up in FF 33
Test with the same script but a different AGS http://services.arcgis.com/rOo16HdIMeOBI4Mb/arcgis/rest/services/Heritage_Trees_Portland/FeatureServer/0
a: Shows up in IE 10
b: Shows up in FF 33

Well, what is the error in this script? The geoportal1's AGS seems to work properly. FF understands the script. But what's not working? If I use L.esri.dynamicMapLayer instead of L.esri.featureLayer the layer shows up. But for further uses I need the L.esri.featureLayer.

I also published it here:

https://github.com/Esri/esri-leaflet/issues/412

Hopefully someone can help!

<!DOCTYPE html>

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

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


texts:
I am using OpenLayers 5 for visualizing some polygons I get from a GeoJSON file.
I would like to style them with markers at the center of each, but no idea how to start. Styling them as polygons is ok.
Now I would like to iterate through all the features, and get their center, so that would be a layer to style.


texts:
I want to re-project and resample a netcf file but I am getting a usage problem:

gdalwarp  -t_srs EPSG:3410 -ts 1383 586 NETCFD:"C:\data\product.nc":gsfc output.tiff


texts: Python split table by "State" field and export to individual tables?

texts:
So I have never posted here before.  I have used the site a lot for help so many thanks to all the great contributors.  So, I have used ArcGIS online to perform batch geocoding (18,000) records and it works great.  Extremely fast and accurate.  It does however chew up a LOT of credits.  This is something I cannot sustain in the long term.  To avoid that credit drain I would like to split my table out by the "State" Column and geocode individual states and then merge back together.  The free online geocoder can handle smaller datasets and by breaking down by state I think I could acheive this.  So far I have:

import arcpy
fc ="C:\GISProjects\ALLMembers\LEE_AllMembers.gdb\LeeMembers_Sort"
field = "State"
myList=set([row.getValue(field) for row in arcpy.SearchCursor(fc)])
for row in myList:
    print myList

This works and returns the set of values.  I guess I"m wondering now how to take those values and arcpy.TabletoTable_Conversion?  From there I could geocode each separate table and merge the results.

any help is much appreciated!


texts: Getting width of road segment given its end points using Leaflet?

texts:
I want to highlight a road segment on a Leaflet Map. I have the coordinates of the end point of the road segment, so could plot it as a linestring. But the problem is how can I get the width of the road segment, so as to highlight it completely?

Is there some Leaflet function / library which can make this kind of queries?

One way could be to use this approximate method  from Obtaining shape or width of road in GeoJSON format?, but is there a better method? Also this method won't work with zoom-in / zoom-out, since the width should change with zoom-in / zoom-out, but hard-coding width of different types of road will not give this effect.


texts: Cutting rectangle set (not squares) into triangles using ArcGIS for Desktop?

texts:
I have some large features that I need to split into smaller pieces for terrain alignment through vertices in another program. I figured triangles would be the best solution.

So basically I want to turn this:

into this:

Sadly the latter image is a mockup. I had some success with the fishnet tool and intersecting to create rectangles, but I can't get triangles out of it. Rotating the net doesn't work because the rectangles aren't square. Any ideas that don't involve manual cutting?


texts:
I have managed to work the QGIS Web Client on my server with the basic buttons and functions.

Now I would like to add "extra" geoprocessing tasks (for example create buffer), based on WPS.

I understand how to add new button to the QGIS Web Client toolbar by the Customizations.js file, but I can not find any working example about how to define the button's function as a geoprocessing task, based on a public WPS (for example Zoo project's WPS capabilities, which I can use in QGIS Desktop).

Is it possible to do? If yes, could you please help me with the basics or show me a working example and code behind it?


texts: Select features of layer which is not on top in OpenLayers

texts:
I have posted this question before on the general StackOverflow forum:
https://stackoverflow.com/questions/16278569/select-features-of-layer-which-is-not-on-top

In short my problem is that I have a baselayer, an overlay layer and a WFS layer in OL.
The base and overlay layers are WMS using GWC in GeoServer.
I want to select a feature of the WFS layer but KEEP the overlay layer on top.

By binding a selectFeature control to the WFS layer moves the layer on top.
Using

map.setLayerZIndex(overlayLayer, map.Z_INDEX_BASE["Feature"]+10);

I move the overlay layer back.
But now I can no longer select a feature of the WFS layer.

Adding the overlayLayer to the SelectFeature control doesn't work because it is a WMS layer.

How can I select a feature on the WFS layer which is not on top?
I've tried the GetFeature control, but that requires a roundtrip to the server and I don't want.

I've made an sample of my problem at http://www.smartdakscan.nl/webapps/assen/selectFeature.html

Any help is much appreciated.


texts: How can I convert Bing's "quadtree" tile addresses to ZXY tile addresses in Python?

texts:
Bing uses a quadtree-based addressing scheme to come up with a single number (or all-digit string) as the address of a tile. Each digit/character (starting from left / most significant) identifies the branch to take along the tree, starting from the root node. This gives one digit per (slippy-map) zoomlevel:

(source: s-msft.com)
Google Maps, OpenStreetMap and others use the ZXY addressing scheme, where a tile is identified by 3 numbers:

the (slippy-map) zoomlevel (z, not to be confused with a z-coordinate like elevation)
the 0-based count of the tile from left (western) edge to the right (east) (x)
the 0-based count of the tile from upper (northern) edge down (south) (y)

Determining the z from a quadtree address is trivial, but how do I calculate x and y of the ZXY scheme from a quadtree address? I'm using Python.


texts: GeoNode credentials not valid for OGC services

texts:
I have GeoNode 2.6.3 installed in ubuntu 16.04 server (standard instalation). The web portal works fine with authentications, view and download permissons, etc. But when i try to access via WFS/WMS services on QGis, Geonode credetials are not valid. I can access with no credentials and see the public layers, and i can access using the geoserver credentials, but not with geonode ones (tried with diferent users).
geoserver.log file shows:

2017-11-13 12:14:47,164 WARN [geofence.cache] - org.geoserver.geofence.cache.CachedRuleReader$NoAuthException: Can't auth user [gpadmin]

(gpadmin is the geonode admin username)

Geofence configuration is as default, and not service security rules added. I can access geoserver adminUI with geonode credentials.
I´ve read this link too but can´t figure out the problem.

Any ideas or tips to find the problem?


texts:
I'm trying to change from Geotools to Geotoolkit but even the last code change is from Jan 7, 2019 its last release is from Oct 27, 2017.

The maven repositories are even older.

How can I use an up to date Geotoolkit version from maven?


texts: Two Different Python Version in QGIS and in Macbook

texts: How can I fix this and make sure that I'm installing Scipy, Matplotlib, GDAL, and Numpy in the correct python version if my QGIS 3.8 uses Python 3.6 and my Macbook Terminal is using Python 2.7? I'm trying to fix this for the SCP plugin.


texts:
I would like to complete the following map tool application that allows only to add some point features, by other functionalities like prompting the user to fill attribute form while editing for each point and save or rollback changes, together with maintaining the same layer styling.

from qgis.core import *
from qgis.gui import *
from PyQt5.QtGui import *
from PyQt5.QtCore import *

class AddPointTool(QgsMapTool):
  def __init__(self, mapCanvas, layer):
    QgsMapTool.__init__(self, mapCanvas)
    self.layer=layer
    self.setCursor(Qt.CrossCursor)

  def canvasReleaseEvent(self, event):
    point=self.toLayerCoordinates(self.layer, event.pos())
    fields=self.layer.dataProvider().fields()
    feature=QgsFeature()
    feature.setFields(fields)
    feature.setGeometry(QgsGeometry.fromPointXY(point))

    self.layer.addFeature(feature)
    self.layer.updateExtents()
    self.layer.triggerRepaint()

canvas=iface.mapCanvas()
layer=iface.activeLayer()

layer.startEditing()
layer.triggerRepaint()
tool=AddPointTool(canvas,layer)
canvas.setMapTool(tool)


texts: How to save (Commit) a temporary table in Mapinfo?

texts:
I tried the following:

Commit Table Query1 As "C:\Temp\XYZ.TAB" TYPE Query Interactive

It does not work.


texts: Search and Replace text in All Fields in QGIS 3

texts:
I want to replace a string in ALL FIELDS / Columns of table without having to do it using field calculator in every field manually using  replace(  "RoadName1" ,'Wharf','WharfXXX').

I have

layer = iface.activeLayer()
searchText = "Wharf"
replaceText = "WharfXXX"
for field in layer.fields():
    fieldName=field.name()
    print (fieldName)
    for feature in layer.getFeatures():
        inText=str(feature[fieldName])
        #print (inText)
        if searchText in inText:
        #try QgsStringReplacement
            QgsStringReplacement(inText,replaceText, False, False)
        #try Field Calculator

            print ("Replaced %s with %s" %(inText,replaceText))
        else:
            #print ("not found")
            pass

It seems to work but I am not seeing any changes in the dataset.

I can't find any example implementations of https://qgis.org/api/classQgsStringReplacement.html

Does this function do what I think it's supposed to -i.e. replace text when it matches?

I guess I could configure field calc to use these inputs but the function above seems cleaner in the long run.

processing.run("qgis:fieldcalculator", {'INPUT':'test.tab|layername=test','FIELD_NAME':'RoadName','FIELD_TYPE':0,'FIELD_LENGTH':10,'FIELD_PRECISION':3,'NEW_FIELD':False,'FORMULA':' replace( \"RoadName\" ,\'Warner\',\'WarnerXXX\')','OUTPUT':'memory:'})


texts:
What am I missing, trying to get Geojson object to get passed with a style, but its not working and there is no console error as well. What do I do?

var olview = new ol.View({
  center: [0, 869099],
  zoom: 7,
  minZoom: 2,
  maxZoom: 20
})

myMap = new ol.layer.Tile({
  source: new ol.source.TileWMS({
        url: 'http://localhost:8080/geoserver/wms',
        params: {'LAYERS': 'basemap', 'TILED': true },
        serverType: 'geoserver'
      })
})
/*`var myMap = new ol.layer.Tile({
      source: new ol.source.OSM()
    });`
*/
// assuming you want to try it without having geoserver, use this

var geoJsonObject = {
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [
          -1.6149902343749998,
          6.118707747190845
        ]
      }
    }
  ]
};

var geojson = new ol.source.Vector({
  features: (new ol.format.GeoJSON()).readFeatures(geoJsonObject)
});

var iconStyle = new ol.style.Style({
        image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
          anchor: [0.5, 46],
          anchorXUnits: 'fraction',
          anchorYUnits: 'pixels',
          src: 'https://openlayers.org/en/v4.2.0/examples/data/icon.png'
        }))
      });

var geojsonVector = new ol.layer.Vector({
  source: geojsonVector,
  style: iconStyle
});

var map = new ol.Map({
  target: 'map',
  layers: [myMap, geojsonVector],
  view: olview
});


texts:
I have installed GRASS GIS 7.6 and when launching it, the next error shows:

ERROR: The APPDATA variable points to directory which does not exist, ask your operating system support
Exiting...
Press any key to continue . . .

How can I fix this error?


texts: Big ImageMosaic is not showing in Geoserver

texts: I have created an ImageMosaic with many big pictures, but in GeoExplorer or GeoWebCache this layer is not displayed. If i create an ImageMosaic with only 3 images (from the big ImageMosaic), then the layer displays. Why?


texts: Inconsistent Value Using Calculate Value in ModelBuilder?

texts:
I have a set of models that use Calculate Value to store file locations as strings for use in model variables and outputs, and am having an issue with the model either hanging on to the first value it calculates for each, or not using the value properly later.

For example, my first model takes an input raster (Input DEM), parses its path (stores as DEMdirstr), feeds the path to a calculate value (expression = r"%DEMdirstr%"; not using the "r" would give us issues with files starting with "t"), which outputs the value as a folder data type to be used in Create Directory (which won't take DEMdirstr on its own), and the value from Create Output Directory is HYA_Workspace, which we attempt to use as in inline variable in output names throughout the model.

If I run the model  by double-clicking it and supply parameters in the dialog box it will create the directory, but the first geoprocessing tool (Pit Remove) will fail because the inline variable for %HYA_Workspace% is being returned as a "1".

Executing (Pit Remove): PitRemove I:\InterpolationResearch\t03\TTR\t03dem.tif 12 1\P2v1_q02\q02_fel.tif

If I run the model by supplying the locations within the model it will complete, but it was ALWAYS use that first value calculated for HYA_Workspace. It will not recalculate that value.

We feel like we've been careful to set preconditions as appropriate. I don't understand why the directory creation would work, but then the directory created couldn't be used by the next step in the model. I'm hoping there's a simple solution we're missing.


texts: Creating feature from empty space in ArcGIS for Desktop?

texts:
If you see below, I would like to create features from the empty space. What is the easiest way to complete this? The freehand tool was sloppy and difficult.

I am using ArcMap 10.2. I have a Basic License, but my company has Advanced available if needed.

The colored areas are part of one feature class, and I would like to include the empty areas as part of that feature class. It can be either one feature or many.


texts: Accessing CSV data directly from website

texts:
I want to do what this person has asked here.  I am using QGIS 3.16
I have copied what was suggested and what seemed to work for the OP, however I am not getting any information displayed.  The data is a CSV file with Ontario COVID data as well as lat and long info.  Excel reads no problem and putting that URL into a web browser will download the CSV file.  Ideally I would like to read it directly as it changes daily.
Also, is there a way to debug this in QGIS?  I would expect to get some kind of error but it just loads the feature name based on filename with no data.  It does append the layer name from the VRT so it must be reading something.
Here is my VRT file:
<OGRVRTDataSource>
<OGRVRTLayer name="COVID Cases">
    <SrcDataSource relativeToVRT="0">"/vsicurl/https://data.ontario.ca/dataset/f4112442-bdc8-45d2-be3c-12efae72fb27/resource/455fd63b-603d-4608-8216-7d8647f43350/download/conposcovidloc.csv"</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="Reporting_PHU_Longitude" y="Reporting_PHU_Latitude"/>
</OGRVRTLayer></OGRVRTDataSource>


texts: Get intersecting and overlapping point between two lines in ArcGIS?

texts:
I have a line feature class that contains two lines with an overlap and intersecting points. The overlap is circled in the picture.

Intersection tool in ArcGIS returns the intersecting point, where two lines cross each other but it is not returning any point at the start of overlap between those two lines.

How can I get all intersecting points and the point where overlapping of two lines start?


texts:
I need to install selenium to automate an ArcGIS Online workflow. I have Python 3x installed with pro in this path:
C:\Program Files\ArcGIS\Pro\bin\Python

These are the subdirectories:
08/03/2020  09:38 PM    <DIR>          .
08/03/2020  09:38 PM    <DIR>          ..
06/09/2020  02:28 PM                35 .condarc
11/12/2019  11:33 AM                 0 .nonadmin
08/03/2020  09:38 PM    <DIR>          conda-meta
11/12/2019  11:34 AM         6,973,824 cwp.exe
08/03/2020  09:38 PM    <DIR>          DLLs
08/03/2020  09:37 PM    <DIR>          envs
08/03/2020  09:38 PM    <DIR>          Lib
08/03/2020  09:38 PM    <DIR>          pkg-metadata
08/03/2020  09:38 PM    <DIR>          pkgs
08/03/2020  09:38 PM    <DIR>          res
08/03/2020  09:38 PM    <DIR>          Scripts

I've tried setting environmental variables, running setx PATH “%PATH%;C:\Program Files\ArcGIS\Pro\bin\Python\Scripts”
And I continue to get an error (pip is not recognized as an internal or external command" when trying to install selenium using pip.
What do I need to do to make pip installs work? Also, is there a place I need to have the selenium and related java sdk files installed to make the pip install work so I can run selenium commands in Python?


texts: How to relate Non-Spatial table to Spatial table in two different Postgis DBs?

texts:
I need to relate a Non-spatial table to Spatial table, but both are in different database.

Spatial_table(dbname:dist)
 gid
 district_name
 district_code
 geom

Non_Spatial(dbname:census)
 ID
 district_code
 population
 male_popu
 female_popu

Can anyone please suggest me, How to relate the above to tables to get the query result for the population of specific district?

Also can anyone tell me about the difference between Joining and Relating of two tables.

Thanks in advance..


texts: Getting GeoServer full datastore layer list by REST

texts:
I have installed GeoServer 2.6.1 with an Oracle datastore which has some of the tables of the database published - the others will be published in the near future.

The documentation says that It is possible to access to the published layers of the datastore by REST GeoServer service like following:

http://localhost:8080/geoserver/rest/workspaces/[WORKSPACE]/datastores/[DATASTORE].html

As far as I know, this REST method only lists all layers of the datastore which are published. I would like to know if could be possible to get all layers (published and non published yet) from a specific Geoserver's datastore.


texts: Calculation of elevation data from raster file(.Tiff) in QGIS

texts:
I am trying to find elevation data of a particular coordinate(lat/long) from a raster file(EPSG:4326) using QGIS C++ APIs. I am able to render a raster file on QgsMapCanvas, which means GDAL is working in my application.  Which APIs should be used to extract elevation value using lat/long value.
Any lead will be valuable even if it is in PyQGIS.


texts: Rasterio read outshape resample does not give correct array

texts:
I have a population grid and want to decrease the resolution by a factor 2.  So I would like to get a raster that has 2 times less rows and columns, a resolution that is 2 times as large as the original one, while retaining the same number of population.
This should be a straightforward operation, but the population numbers change incorrectly if I follow the documentation. When the resampled array is written to file, the numbers are incorrect.
import rasterio
import numpy as np
import matplotlib.pyplot as plt
from rasterio.enums import Resampling

factor = 2

with rasterio.open("mwi_ppp_2020.tif") as dataset:

    arr = dataset.read(1,masked=True)
    #Show total population of Malawi
    print(arr.sum())

    data = dataset.read(1, masked=True,resampling=Resampling.bilinear,
        out_shape=(
            int(dataset.height * factor),
            int(dataset.width * factor)))

     #Show total population of array
    print(data.sum())

     # scale image transform
    transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[1]),
            (dataset.height / data.shape[0])
                                )

Shouldn't the resample method yield an array with resampled cell values to get the same number of total population? Any tips and hints would be appreciated!


texts: Aggregates Functions - QGIS - Calculating Area Explained

texts:
I've been reading through question and answers on the Aggregates Function, but have been unable to obtain a clear understanding and solution to my problem. I have a polygon grid ("Grid 4km") and I have another vector polygon layer ("S1 Grid Intersection Strata"). I want to calculate the area of the polygon layers that are within a grid cell. I intersected the vector polygon layer with the grid so that you can select on polygons within the grid.

The steps I have taken so far:

I calculated the $area in a field called Area in "S1 Grid Intersection Strata".
I worked in the field calculator for "Grid 4km" to create a new field called "S1 Area" using the following:

aggregate('S1 Grid Intersection Strata', 'sum', "area", filter:=intersects($geometry, geometry(@parent)))

This seems to work, but when I select on individual examples and calculate by hand the sum of the area is not what is showing up in the new "S1 Area" field. The image below shows the polygons selected within the grid with the area values. The sum of the areas in the "S1 Grid Intersection Strata" do not add up to the value in the "Grid 4km" "S1 Area" field.

I am hoping that someone might be able to provide a good explanation of what the @geometry component in the formula above is doing and how I can fix this to sum up the area of the intersecting polygons in the grid.


texts:
How do I get CMV to start up more efficiently.
Exploring ArcMap GENERAL start services immediately check box (to un check it).
ArcMap Caching is set to Dynamic from data, levels 9 - 17, 5MB estimated RAM.
Looking for ways to delay start of CMV layers and widgets so CMV appears sooner.


texts: What is the altitude difference between WGS84 and NAD27 Zone II?

texts:
I have a point cloud in a piece of software and have exported the point in both WGS84 and NAD27 State Plane California Zone II (in US survey feet). The WGS data looks like this, with elevation in meters:

-121.330272 38.547287 5.661467
-121.330272 38.547287 5.646508
-121.330272 38.547287 5.661565
-121.330272 38.547287 5.657426

where as the NAD27 Zone II data looks like this:

2191840.924977 321431.770306 44.012074
2191840.804111 321431.731749 43.962994
2191840.798730 321431.632417 44.012394
2191840.725906 321431.877987 43.998812
2191840.639473 321431.718561 44.001784

I don't understand why the elevation is different and what units it is meant to be in. It's not a simple change from meters to feet. Something else is going on?

EDIT:

I found out that the projection file for the NAD27 Zone II projection used is this in case that helps:

PROJCS["NAD27 / California zone II",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982138982]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",39.83333333333334],PARAMETER["standard_parallel_2",38.33333333333334],PARAMETER["latitude_of_origin",37.66666666666666],PARAMETER["central_meridian",-122],PARAMETER["false_easting",2000000],PARAMETER["false_northing",0],UNIT["Foot_US",0.30480060960121924]]

EDIT:

Example code showing difference in converting first point of data.

import pyproj

p1 = pyproj.Proj("+init=EPSG:4326")
p2 = pyproj.Proj("+init=EPSG:26742")

a = (-121.330272, 38.547287, 5.661467)
b = (2191840.924977, 321431.770306, 44.012074)

print "Source point", a
print "Target point", b
print 'Converted: ', pyproj.transform(p1, p2, *a) # different :()

The two coordinate systems are:

Geographic Coordinate System: WGS 84 (EPSG::4326)
Geodetic Datum: World Geodetic System 1984 (EPSG::6326)
Ellipsoid: WGS 84 (EPSG::7030)
Prime Meridian: Greenwich (EPSG::8901)
Angular Units: degree (EPSG::9102)

and

Projected Coordinate System: NAD27 / California zone II (EPSG::26742)
Projection Method: Lambert Conic Conformal (2SP)
    Latitude of false origin: 37.6667
    Longitude of false origin: -122
    Latitude of first standard parallel: 39.8333
    Latitude of second standard parallel: 38.3333
    Easting of false origin: 2e+06
    Northing of false origin: 0
Geographic Coordinate System: NAD27 (EPSG::4267)
Geodetic Datum: North American Datum 1927 (EPSG::6267)
Ellipsoid: Clarke 1866 (EPSG::7008)
Prime Meridian: Greenwich (EPSG::8901)
Linear Units: US survey foot (EPSG::9003)


texts:
I am working in a project in which I need to use PyQGIS and H5PY. The problem comes  when I load qgis library and h5py, for example:

import qgis
import h5py

I get this error:

from ._conv import register_converters as _register_converters
  File "h5py\h5r.pxd", line 21, in init h5py._conv
  File "h5py_objects.pxd", line 12, in init h5py.h5r
  File "h5py_objects.pyx", line 1, in init h5py._objects
ImportError: DLL load failed: No se encontró el proceso especificado.

I am using QGIS 3.2 and H5py 2.8.0 with Python36 interpreter contained in QGIS installation folder. H5py was installed through the command:

python -m pip install h5py

I tried other versions of H5py (2.6.0, 2.7.1) but I had the same problem.

Is there anyone who knows how to solve this dll interaction problem?


texts:
I'm using GeoPandas to create a categorical choropleth of my data. I managed to get the correct colormap applied to the polygon fill but there are tiny slivers between the polygons that I would like to remove. What is the appropriate way to do this? I thought about the following options:

Apply the same color (same as used for fill) to each polygon line
Make the lines transparent
Set the line width to 0

This is my code so far:

fig, ax = plt.subplots()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin , ymax)
gdf.plot(ax=ax,column=column,cmap=cmap,categorical=True,linewidth=0)
plt.axis('off')
fig.savefig(output_filename+".pdf",dpi=300)
fig.savefig(output_filename+".jpg",dpi=300)

Gdf is a geodataframe with multiPolygons. The result looks like this. As you can see, there are weird slivers between the polygons. The shapefile does not contain any gaps.


texts:
I'm quite new to GDAL and to the remote sensing world so sorry in advance if I say something silly :-)
What I would like to do is to assign from scratch the EPSG 3857 System Of Reference to a Bing tile (a not georeferenced jpg) downloaded from the Bing server.
I've already computed the coordinates [x, y] this tile should have in the 3857 SoR and create a world file from them, but how can I now project my tile to 3857 by using this .tfw?

Thank you so much in advance!

Martina


texts: Diagonal white lines in base maps in Openlayers 3

texts:
What is the reason that Openlayers 3 (3.20.1 to be exact) renders diagonal white lines on base maps such as Bing, Google etc under Chrome? Things are OK under Firefox and Edge. (Chrome version is: Version 56.0.2924.87 (64-bit))

Please check the below image:


texts:
I have rasterdataset in gdb. I create layer from this rasterdataset.

    private static void CreateLayer()
    {
        ...
        IRasterLayer rlMosaicDataset = new RasterLayer();
        rlMosaicDataset.CreateFromDataset(rasterDataset);// name in gdb

        IRasterRGBRenderer stretchRenderer = new RasterRGBRenderer();
        IRasterRenderer rasterRenderer = (IRasterRenderer)stretchRenderer;
        IRaster raster = rasterDataset.CreateDefaultRaster();
        rasterRenderer.Raster = raster;
        rasterRenderer.ResamplingType = rstResamplingTypes.RSP_BilinearInterpolation;  // smooth
        rasterRenderer.Update();
        IRasterStretch stretchType = (IRasterStretch)rasterRenderer;
        stretchType.StretchType =
            esriRasterStretchTypesEnum.esriRasterStretch_ESRI;
        stretchType.StandardDeviationsParam = 2; // what is?
        rlMosaicDataset.Renderer = rasterRenderer;

        rlMosaicDataset.MaximumScale = 0.1;
        rlMosaicDataset.MinimumScale = 500000000;
        ...
    }

This layer is black-and-white imagery. I'd like that its color was green and red.
(In ArcMap I can do: Layer Properties -> Symbology -> Stretched -> Color Ramp -> here I can choose strip with green, yellow, orange, red colors)


texts: Debugging Global name ''QgsMapLayerComboBox' is not defined?

texts:
I try to develop little plugin for Qgis. I used custom widgets in QT that comes with QGIS 2.18 and loaded QgisMapLayerComboBox. But when I reload my plugin in QGis it gives me an error:

  2018-01-14T00:06:10   1   Traceback (most recent call last):
              File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python\qgis\utils.py", line 342, in startPlugin
                plugins[packageName].initGui()
              File "C:/Users/Michal.DESKTOP-M9QN5BH/.qgis2/python/plugins\DMS2Point\DMS2Point.py",
  line 172, in initGui
                parent=self.iface.mainWindow())
              File "C:/Users/Michal.DESKTOP-M9QN5BH/.qgis2/python/plugins\DMS2Point\DMS2Point.py",
  line 139, in add_action
                self.dlg = DMS2PointDialog()
              File "C:/Users/Michal.DESKTOP-M9QN5BH/.qgis2/python/plugins\DMS2Point\DMS2Point_dialog.py",
  line 30, in init
                self.setupUi(self)
              File "C:/Users/Michal.DESKTOP-M9QN5BH/.qgis2/python/plugins\DMS2Point\DMS2Point_dialog_base.py",
  line 48, in setupUi
                self.mMapLayerComboBox_2 = QgsMapLayerComboBox(self.widget)             NameError: global name 'QgsMapLayerComboBox' is not defined

I don't understand the last line. Code is compiled and everything, step by step made according to tutorials 'How to make a Qgis plugin...'. This error does no appear when I use 'normal' QT ComboBoxes. I wanted to use Qgis Custom Widgets because they might be easier in use and modification.


texts:
I'm looking for automated/scripted ways to generate PNG images based on GeoJSON multi-polygons.

I have multiple polygon isochrones for a given starting point, e.g. 10 minutes, 15 minutes, 30 minutes, etc. These are generated in GeoJSON format. Are there any techniques for automating the generation of an output like this?

i.e. I want to be able to generate this type of colourful PNG maps automatically, as I have over 1000 to do.


texts: Simple software for mapping residential address in Singapore

texts:
I am new to GIS and I recently have a project where we will like to map clinic addresses vs patients' address, to show the distance between these locations.

In addition, at each clinic address, a number of patients would have visited for either chronic or acute visit. I will like to be able to illustrate the different type (acute or chronic) of visits, as well as the number of visits (9 acute vs 5 chronic for Clinic XXX).

Lastly, we will need to do this on a regular basis (meaning we will have a continuous feed of address that we will key in periodically), and also to layer the different maps (options to select, unselect layers).

Is there any simple freeware or software out there that can readily do this? My colleague has recently learnt ArcGIS, but I feel that it's an overkill.


texts: How to read Arc/Info e00 vector files (no binaries) with Python?

texts:
I can easily read e00 files with QGIS and it seems as if in older versions of ogr there was driver support for the arc/info e00 file format. Nowadays that seems to be deprecated, as running ogrinfo --formats doesn't list arc/info.

I also found this useful link, a general description of the file format. Based on that information I would be able to write my own program to put all the pieces together. But hey, this must have been done before, right?

Can somebody point me into the right direction?


texts: Find regions of high density of points of the same class

texts:
I have a dataset consisting of geo points of two classes. For simplicity we can call one class "successes" and the second one "failures".

What I want to do is to find regions (areas, polygons) where ratio of failures to successes is relatively high.

If I had one class, I would do something like DBSCAN or OPTICS to find clusters.

I thought about separating whole region into smaller pre-defined sectors, computing ratio of failures to successes in each of them, and then to do LISA (Local Indicators of Spatial Association) clustering on them, but then my result will depend on the way I separate region into sectors.

How can I approach this problem?

Also, what if I wanted to divide whole region into sectors of several classes ("high probability of failures", "medium probability of failures" etc)?

I generally use Python and PostGIS but open to other tools. I dont have ArcGIS license.


texts: QQIS composer: Using bold and normal text styling in a label at the same time?

texts:
I am writing text explanations for map symbols and need to write three lines for each symbol -  basically explaining the same thing in three different languages.

I'd like to bold the main language - but when I select bold from the composer 'font' menu, it bolds all three lines - not just the one I've highlighted.

Could it be there is no way to mix stylings, font sizes and so on in a single composer label? It would be very handy option if its possible to have.


texts:
I searched all the internet and I didn't find a solution to my problem. I want to convert a JTS geometry to a GML3.1 geometry. Geo tools says that have support but the have not (please check encode method from GML class, where the GML3 version is not checked). Their examples on official page are not working at all because ... guess what ... exactly, encoding for GML3 is not implemented (I also tried with 13-RC1, not only with 12.2). Any other suggestion to do this ? Or I will end up by implementing something on my own that will generate the GML3.1 xml ?

PS: I am just a little desperate. :)