texts: Refreshing display for newly added labels with ArcObjects from Python?

texts:
I've added a shapefile to an MXD and Edited the symbolization of it with ArcObjects.  Now I'm trying to edit the labeling with ArcObjects.  I've successfully been able to edit the default class in the IAnnotateLayerPropertiesCollection2 object.  I've also successfully added two more label classes to that collection object.

The default class that I edited is displaying correctly, but the extra two Label Classes aren't displaying after the script is run.  When I manually open the properties of the layer and open the properties for each label class that I added and choose OK or Apply, then they show up.  So I know that they are being added properly, but how do I get them to display through the script?

I've tried MxDoc.UpdateContets(), MxDoc.ActiveView.Refresh(), arcpy.RefreshActiveView() and none of them work for this.

I'm using Python to handle the ArcObjects, but if you know how to do it in VB or C#, I can interpolate it.

def ApplyLabels_arcobjects(layersource, pFLayer, pStyleGallery, carto):
    # If there are no labels in the custom symbology object
    if not layersource.symbology.labels:
        return

    # Get Annotation Properties Collection
    geofeaturelayer = Snippets.CType(pFLayer, carto.IGeoFeatureLayer)
    pAnnoPropColl = geofeaturelayer.AnnotationProperties
    pAnnoPropColl2 = Snippets.CType(pAnnoPropColl, carto.IAnnotateLayerPropertiesCollection2)

    # Loop through each label in the custom symbology object
    i = 0
    for labelsource in layersource.symbology.labels:
        print("This is label number {0}".format(i+1)) # TODO Remove, for debugging

        # Get the premade Label Style from the style file
        pStyleGallery.LoadStyle(labelsource.style_filename, labelsource.label_classname)
        gallerystyle = getLabelSymbol(labelsource, pStyleGallery, carto)
        pLabelStyle = Snippets.CType(gallerystyle, carto.ILabelStyle2)

        # If I can get a Annotation Property that is already there, great:
        if i <= pAnnoPropColl2.Count - 1:
            print("Editing the one already there")  # TODO Remove, for debugging
            pAnnoLayerProp = pAnnoPropColl2.QueryItem(i)[0]

            # Class Name
            pAnnoLayerProp.Class = labelsource.name

            # Sql string for Where expression
            pAnnoLayerProp.WhereClause = labelsource.sql_string

            # Set Sympol and Placement Properties
            pLabelEngineProp = Snippets.CType(pAnnoLayerProp, carto.ILabelEngineLayerProperties2)
            pLabelEngineProp.Symbol = pLabelStyle.Symbol
            pLabelEngineProp.OverposterLayerProperties = pLabelStyle.OverposterLayerProperties

        # If I have to add an Annotation Property to the Collection, then:
        else:
            print("Adding new Label Engine")# TODO Remove, for debugging
            # New Label Engine
            pLabelEngineProp = Snippets.NewObj(carto.MaplexLabelEngineLayerProperties, carto.ILabelEngineLayerProperties2)

            # Apply Symbol
            pLabelEngineProp.Symbol = pLabelStyle.Symbol

            # Apply Placement Properties
            pLabelEngineProp.OverposterLayerProperties = pLabelStyle.OverposterLayerProperties

            # Cast
            pAnnoLayerProp = Snippets.CType(pLabelEngineProp, carto.IAnnotateLayerProperties)
            # Class Name
            pAnnoLayerProp.Class = labelsource.name
            # Sql string for Where expression
            pAnnoLayerProp.WhereClause = labelsource.sql_string
            # Display Annotation
            pAnnoLayerProp.DisplayAnnotation = True
            # Add it to the collection
            pAnnoPropColl2.Add(pAnnoLayerProp)

        # Turn on Annotation for the Map Layer
        geofeaturelayer.DisplayAnnotation = True

        i += 1


texts: LayerSwitcher() not working with Group Layers from Geoserver

texts:
i have used LayerSwitcher() before and it has worked fine but I am now trying to use it with the group layers in Geoserver and the map blanks out.

tiled = new OpenLayers.Layer.WMS(
                    "Geoserver layers - Tiled", "http://niles.iplantcollaborative.org:80/geoserver/wms",
                    {
                        LAYERS: 'Shapefile.py',
                        STYLES: '',
                        format: format,
                        tiled: !pureCoverage,
                        tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
                    },
                    {
                        buffer: 0,
                        displayOutsideMaxExtent: true,
                        isBaseLayer: true,
                        displayInLayerSwitcher: true
                    }
                );
var imagery = new OpenLayers.Layer.WMS(
                "Global Imagery",
                "http://maps.opengeo.org/geowebcache/service/wms",
                {layers: "bluemarble"}
            );
                map.addLayers([tiled, imagery]);

                // build up all controls
                map.addControl(new OpenLayers.Control.PanZoomBar({
                    position: new OpenLayers.Pixel(2, 15)
                }));

                map.addControl(new OpenLayers.Control.Navigation());
                map.addControl(new OpenLayers.Control.Scale($('scale')));
                map.addControl(new OpenLayers.Control.MousePosition({element: $('location')}));
                map.addControl(new OpenLayers.Control.LayerSwitcher());
                map.zoomToExtent(bounds);

The moment I comment out the LayerSwitcher(), the map is rendered fine.

Also, I tried adding a new WMS, imagery (in the code above), and if I try to do the following,

map.addLayers([imagery, tiled]);

I am only able to see a broken image. Can we not add a layer through code with group layers?

Where am I going wrong?


texts: Using GRASS GIS/QGIS to georeference TIFFs in 3D space, vertically?

texts:
I am a geophysicist working for a small shallow surface geophysics firm. We create 2D slices of the ground showing data collected from various surveys. These are generally cross sections (TIFFs) of uneven shapes.

This is what i want to do using QGIS/GRASS:

.

Is there a way to use GRASS GIS/QGIS to georeference these TIFFs in 3D space, vertically? If so is there a guide anywhere?

Would I have to use QGIS as part of workflow using another program (like sketchup)?

(the image is created in Discover, which I want to move away as we start to intergrate into QGIS.)


texts: Modifying a projection ArcMap 10.0 Arcpy

texts:
I have the need to project data from wgs84 to UTM 15N NAD1927 FOOT_US.

I would like to have arcpy do this for me.  I can project the data using one of ArcMaps predifined prj files in the Coordinates folder..but how would I do the "Modify" piece?

UTM 15N NAD1927 is defined in the prj file as "Meters" but I would like to do the arcpy.project_management using a prj file that has been modified either in ArcMap and then saved to disc or my another method...

I found this but am a little unsure.

Thanks for any insight or solution to using a modified prj file when doing a arcpy.project_management on shapefiles.


texts: How to select fields returned by GeoServer WMS GetFeatureInfo request?

texts:
I'm running simple local mapping service based on PostGIS/GeoServer/OpenLayers stack.

Following opengeo tutorial I've tied click events in OpenLayers to WMS GetFeatureInfo call using:

    map.events.register('click', map, function (e) {
        var url = "http://localhost:1979/geoserver/wms"
          + "?REQUEST=GetFeatureInfo"
          + "&EXCEPTIONS=application/vnd.ogc.se_xml"
          + "&BBOX=" + map.getExtent().toBBOX()
          + "&X=" + e.xy.x
          + "&Y=" + e.xy.y
          + "&INFO_FORMAT=text/html"
          + "&QUERY_LAYERS=MYDATA:MYLAYER"
          + "&LAYERS=MYDATA:MYLAYER"
          + "&FEATURE_COUNT=5"
          + "&SRS=EPSG:3857"
          + "&STYLES="
          + "&WIDTH=" + map.size.w
          + "&HEIGHT=" + map.size.h;
        window.open(url,
          "getfeatureinfo",
          "location=0,status=0,scrollbars=1,width=800,height=125"
        );
      });

I've used GeoServer's FreeMarker templates to format the response correctly. The problem I have is that the response I get returns all of the fields from the layer.

Is there a way to select the ones that are relevant to the user?


texts: Suddenly unable to save ArcGIS Online WebApp?

texts:
I'm putting together a quick webapp for a client using ArcGIS Online's webapp builder. However, it all of a sudden will not allow me to save any of my changes to the webapp. The save button and all of the options in the expandable menu are grayed out. As far as I know, it is like this even nobody has it open (although I believe this shouldn't be a problem anyway).

I have been unable to find an answer through ESRI's site. Anyone have any ideas?


texts: Debugging errors of SAGA interpolation algorithms in QGIS Processing?

texts:
I encounter recurrent problems with the geoalgorithms that perform interpolations.
No matter what technique I use I end with the following errors. As stated in the OGR tab QGIS 2.0 is unable to open a file in C drive.
I have very little IT background so the message error is meaningless.

Is there anybody able to decrypt the message and solve the problem?

here this is the algorithm i used


texts: Downloading Sentinel-2 data from Google Cloud Storage with Python?

texts:
I'm trying to download Sentinel-2 data from the Google Cloud Storage and basically adapted the FeLS-module (great work, by the way!). As it is done in the module, I also download the index.csv first (latest version here, but be careful... it's huge!) to search for scenes that match my requirements. But after several successful tests, I figured that there are scenes available within the bucket that are not listed in the csv.

I don't know why and have already contacted Google about this, so now I'm looking for another solution: is there a way to establish a connection to the tile I'm looking for and then list all available subfolders? For example, I would like to do something like this:

import urllib2
from bs4 import BeautifulSoup

url = 'https://console.cloud.google.com/storage/browser/gcp-public-data-sentinel-2/tiles/39/R/YH/'
page = urllib2.urlopen(url)
soup = BeautifulSoup('html.parser')
# the following has been working for another link (https://landsatonaws.com/L8/001/003/), ...
# ... this is just to make clear what I want to do:
table = soup('table')[0].find_all('td')
scenes = [table[i].string for i in range(0, len(table), 3)]

For the given Google Cloud Storage url, I never get something back except for an HTTPError: HTTP Error 404: Not Found.

If I could get rid of searching the index.csv, I would truly get all scenes that are available. Is this possible somehow?


texts: How to set resolution according to a given scale in OpenLayers

texts:
This is a followup to this question.

What I know is:

given the papersize, I can calculate how many pixels for a given DPI fill a millimeter (dpi * 1/25.4) =pixelpermilimeter
the number of pixels needed to represent one meter could be calculated via getPointResolution()*proj.getMetersPerUnit() =pixeltomapmeter
the scale is given as 1000 * pixelpermilimeter * pixeltomapmeter

From the issue mentioned above, I learned, that with a given papersize and a given screensize, I could set the resolution to fit the paper as follows:

map.getView().setResolution(resolution / scaling);

where resolution is the map resolution and scaling as the quotient of one map dimension in relation to the paper dimension.

What I am looking for is a way to determine the correct mapsize / resolution /  scaling such, that a given scale is printed on paper.

Say the scale is 1:10000, I need the center of the map printed such, that 1 cm represents 100m.

Does that make sense?

I think, I have the dots, but I am not able to connect them in the right way.


texts: Counting categorized points in polygons using PostGIS

texts:
I want to do an operation i PostGIS similar to the count points in polygon in QGIS.
I have

A multipolygon table (plan.polygon) consisting of several areas
A point table (basemap.points)

in a PostGIS database.
The point layer is categorized by the column code_text.
I want to count the number of each category for each of areas resulting in a multipolygon table with rows areas and columns with the count of each category.
Below you can se the areas (red lines) and the points with categorization (red, yellow, green dots).

I've been trying the following resulting in a query processing forever.
SELECT b.gid, b.the_geom,
    (SELECT COUNT(v.code_text)
    FROM basemap.points v, plan.polygon x
    WHERE entity = 120 AND ST_within(v.the_geom, x.the_geom)
GROUP BY a.code_text)
AS count_category_a,
    (SELECT count(y.code_text)
    FROM basemap.points y, plan.polygona z
    WHERE entity = 130 AND ST_within(y.the_geom, z.the_geom)
GROUP BY a.code_text)
AS count_category_b
FROM basemap.points a, plan.polygons b;


texts: How to design a web application using georeferenced aerial imagery

texts:
I'm working on a web application that monitors the state of plants in their respective fields, of which we are taking georeferenced aerial imagery with a drone.

Unfortunately I'm new to this technology, but what I would like to do is to present those fields via a web application displaying points on those maps each of which will represent a plant.

Do I need to store all the aerial images on the local server and install some kind of provider locally or should I load them on a gis service and access them via an API and add my points on the map received?


texts: Rendering polygons outside latitudes [85,-85] in OpenLayers 3

texts:
I have some polygons I am drawing on an OpenLayers 3 map.

Some of these are above latitude 85, and used to be rendered fine with OpenLayers 3 v3.15.1: https://jsfiddle.net/adlerhn/98mn4w5z/

However, in the latest version of OpenLayers (3.20.0), there is clipping in place, and polygons outside the extent of the map are not rendered: https://jsfiddle.net/adlerhn/xt5m1o6d/

The reason I am trying to render outside of the bounds of a map is because I am using ol3-cesium to render into a 3D globe, and in there you can see latitudes between 85 and 90, which would otherwise appear as a gap on the top.

I have tried adding the extent property to the vector layer, but I haven't got any result.

var extent = ol.proj.get("EPSG:4326").getExtent();

var vectorLayer = new ol.layer.Vector({
    map: this.map,
    source: new ol.source.Vector({
        features: [new ol.Feature({
            geometry: geometry
        })]
    }),
    extent: extent,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'rgba(0, 128, 255, 0.5)'
        })
      })
});

How can I get OpenLayers to render outside the extent of the map?


texts: How does ModelBuilder compare with Vensim and Stella for environmental modeling?

texts: Does anyone have any advice or has anyone experienced any limitations in ArcGIS software when trying to use model builder for any type of environmental phenomenon? I know that Vensim and Stella can be used with more detail for just about anything when modeling an environmental phenomenon. What kind of limitations might there be with ArcGIS software?


texts: Converted Vector to Raster file is Black and White in colour (GDAL Rasterize)

texts:
def main1(shapefile):

    #making the shapefile as an object.
    input_shp = ogr.Open(shapefile)

    #getting layer information of shapefile.
    shp_layer = input_shp.GetLayer()

    #pixel_size determines the size of the new raster.
    #pixel_size is proportional to size of shapefile.
    pixel_size = 1

    #get extent values to set size of output raster.
    x_min, x_max, y_min, y_max = shp_layer.GetExtent()

    #calculate size/resolution of the raster.
    x_res = int((x_max - x_min) / pixel_size)
    y_res = int((y_max - y_min) / pixel_size)

    #get GeoTiff driver by
    image_type = 'GTiff'
    driver = gdal.GetDriverByName(image_type)

    #passing the filename, x and y direction resolution, no. of bands, new raster.
    new_raster = driver.Create(output_raster, x_res, y_res, 1, gdal.GDT_Byte)

    #transforms between pixel raster space to projection coordinate space.
    new_raster.SetGeoTransform((x_min, pixel_size, 0, y_min, 0, pixel_size))

    #get required raster band.
    band = new_raster.GetRasterBand(1)

    #assign no data value to empty cells.
    no_data_value = -9999
    band.SetNoDataValue(no_data_value)
    band.FlushCache()

    #main conversion method
    gdal.RasterizeLayer(new_raster, [1], shp_layer, burn_values=[1])

    #adding a spatial reference
    new_rasterSRS = osr.SpatialReference()
    new_rasterSRS.ImportFromEPSG(2975)
    new_raster.SetProjection(new_rasterSRS.ExportToWkt())
    return gdal.Open(output_raster)

This is my code for converting vector .shp to raster .tif. It converts the shapefiles into .tif files BUT the output is only black and white in colour.

This is my input shapefile.

This is the output vector.

So what changes do I make in the code so that my output file has similar colour?

This is the current symbology properties:

This is when I select multiband colour.


texts: Accessing non-spatial PostgreSQL tables using OGR Python bindings?

texts:
When connecting to PostgreSQL/PostGIS using the OGR Python bindings is it possible to get non-spatial table names (i.e. OGR layers) from the connection?

Currently, doing:

conn = ogr.Open("PG: #params")
for layer in conn:
  print layer.GetName()

Will only print the names of spatial tables, however I can access non-spatial tables by (for example) directly specifying the table name, e.g.:

layer = conn.GetLayerByName('non_spatial_table_name')

The OGR docs refer to the "PG_LIST_ALL_TABLES=YES" parameter, but I can't find how this links to API (let alone the Python SWIG bindings!).


texts:
I try to label a layer by using a rule based labeling. What I specifically try to do is to find an expression which says:

"If column 1 shows the attribut x then label layer "Test" by the attribute of column 2"

edit:
I'm using QGis Version 3.10.0-A Coruña.


texts: Free implementation of Polynomial Approximation with Exponential Kernel (PAEK) algorithm

texts:
I need a Free (i.e. Open Source) implementation of the Polynomial Approximation with Exponential Kernel (PAEK) algorithm, preferably in C, C++, Python, Julia or R.

The algorithm is one of the methods ArcGIS offers for line smoothing and described in Bodansky, Eugene; Gribov, Alexander; and Pilouk, Morakot, "Smoothing and Compression of Lines Obtained by Raster-to-Vector Conversion", LNCS 2390, Springer, p. 256-265, 2002.

If I can't find a ready to use implementation, I will have a go at it myself.


texts: Which HTML standard supports the print composer of QGIS?


texts:
I haven't been able to find a solution with the source so am trying to work on the output file.  I have also asked this question on SO

I have a large (200MB) geoJSON file that has a lot of complex polygons and multipolygons.  A very truncated example is at https://gist.github.com/jinky32/81f61e1fc118822ba103?short_path=d16949b

As you can see this file is comprised of polygons and multipolygons that have a String property of either 1 or 2.  Below is an example of how these shapes look on mapshaper.org when highlighting a multipolygon of either value in the same tile (essentially c.90+% of this tile is made up of a multipolygon with one value or the other)

I do not need to differentiate between these different values and polygons / multipolygons with a String value of either 1 or 2 can be combined together which I hope will reduce the file size.

Can anyone advise how I can achieve this - preferably with a cli tool?


texts: Eliminating sliver polygons automatically

texts:
I have thousands of sliver polygons outside of my study area. What is the best way to eliminate them automatically? In ArcGIS or QGIS.
Those polygon islands are an example of the slivers.
Note: I can´t select by attributes (small areas) because then it would eliminate other polygons inside the study area.


texts: Joining two vector layers with polygons and then selecting one via other in QGIS?

texts:
This could very well be a duplicate question, but i can't find the answer in the site.

Using QGIS 3.4, I have two vector layers with polygons and same projections: one with all individual watersheds and the other with cumulated upstream watersheds. The first layer has a field that can be joined with the second, so that each watershed has a corresponding upstream area. The second layer is designed for this specific purpose.

I want to join the layers and then select any watershed, resulting in a selection of that watershed's cumulated upstream area. Exactly what the second layer is designed for.

I can join the layers in the properties tab, and the first layer's attribute table is populated with the second's data. However, when I select the first layers watersheds, nothing is selected in the second layer.

What is the next step I must take to achieve my selecting results?

The description of the second layer states that the polygons are overlapping each other.

Watershed layer: https://www.smhi.se/polopoly_fs/1.24469!/Menu/general/extGroup/attachmentColHold/mainCol1/file/aro_y_2012_2.zip

Upstream area layer: https://www.smhi.se/polopoly_fs/1.31173!/Svaro.zip

Both layers are in SWEREF99TM


texts: Geoserver GetFeatureInfo WMS attributes: how to hide geometry?

texts:
I know this has been discused before, but before somebody tell me about data base views or HTML templates, this is my concern:

Geoserver returns the geometry field from WMS getfeature info, enabling anyone to create scripts to download the data from the server as if it was WFS.

Is there a way to prevent this behaviour?


texts: Filtering fields from a PostgresSQL/PostGIS database layer in QGIS by Query builder

texts:
I incurring in an error while trying to make the following filter of fields in the layer properties Query builder in QGIS 3 from a PostGIS database.

select "field_A", "field_B2", "field_D" from my_database

what should be the correct sentence to filter fields in the Query builder?


texts: Use ArcPy for finding value throughout all fields and export to new layer

texts:
I need to find a way to find a specific value in the attribute table and then export all rows with said value to a new feature class.

Below is the code that I have managed to come up with (using other answers on this site). It gives the FID of the rows with the value searched but I can't find a way to select and export those records to a new file.

import arcpy
arcpy.env.overwriteOutput = True

# Set the environment
workspace = r"\\C:\\GIS\\"
arcpy.env.workspace = workspace
fc = "species.shp"

# Get fields in the attribute table
field_names = [f.name for f in arcpy.ListFields(fc, field_type='String')]
print(field_names)

# Loop through fields
for field in field_names:
    # Create a search cursor through rows in the feature class that looks for "In Danger" values and prints the record's FID number
    search = arcpy.da.SearchCursor(fc, "*")
    for row in search:
        if "In Danger" in row:
            print("In Danger value for FID {0}".format(row[0]))
            arcpy.Select_analysis(fc, 'new_fc', '{} like \'%{}%\''.format(field, 'In Danger'))


texts: Dissolving polygons based on spectral differences from QGIS?

texts:
I'm looking for a package or algorithm that allows to dissolve polygons based on their spectral characteristics and a threshold defined by the user in QGIS. Equivalent to what the "Spectral Differences" algorithm does in eCognition.

I did a segmentation using the Meanshift algorithm of the OTB in QGIS, based on a three band image of Sentinel 2 satellite. I want to analyze some oversegmented polygons, and dissolve them with their neighbours if their spectral differences are less than a certain threshold.

Does anyone knows of a tool that does this kind of work? or something similar at least? I haven't found anything in the web, but I suppose there might be someone who has already done it.


texts: How do you create a dynamic legend in ArcGIS 10x that summiaries the areas of a polygon based on a layout?

texts:
Hi I have a land use polygon layer with around 700 maps to produce. Am trying to create a dynamic legend that summarises the areas M2 based on the layout extent. The layout extents and maps are produced using data driven pages but I dont want to have to manually create 700 legends of summarised data.
I need the legend to display the area of each landuse type based on the map extent.
I have five land use types:
Grass
Shrubs
Water
Hard standing
Wildflower Meadow
Any help would be gratefully received.
Many thanks James


texts: Create a new raster TIFF file which is masked based on the GeoJSON file

texts:
I was trying to mask the building given its TIFF file and GeoJSON file  creating a new raster TIFF file. I don't know what is the problem, the result is all black, here is the code below:

from osgeo import gdal, ogr

def create_building_mask(rasterSrc, vectorSrc, npDistFileName='',
                        noDataValue=0, burn_values=1):

    source_ds = ogr.Open(vectorSrc)
    source_layer = source_ds.GetLayer()

    srcRas_ds = gdal.Open(rasterSrc)
    cols = srcRas_ds.RasterXSize
    rows = srcRas_ds.RasterYSize

    memdrv = gdal.GetDriverByName('GTiff')
    dst_ds = memdrv.Create(npDistFileName, cols, rows, 1, gdal.GDT_Byte)
    dst_ds.SetGeoTransform(srcRas_ds.GetGeoTransform())
    dst_ds.SetProjection(srcRas_ds.GetProjection())
    band = dst_ds.GetRasterBand(1)
    band.SetNoDataValue(noDataValue)

    gdal.RasterizeLayer(dst_ds, [1], source_layer, burn_values=[burn_values])

    dst_ds = 0

    return

create_building_mask("img.tif", "img.geojson", npDistFileName='mask.tif', noDataValue=0, burn_values=1)

The result of mask.tif is all black, instead of masked img.


texts:
I'm drawing marker line with user defined geometry of the marker (see example on the picture). Markers are placed with interval of screen (not map) units.  It's possible to do with svg-marker or font marker. Can't find way to do it with geometry generator.
Symbols on the picture are for example. I have a lot of user defined symbols, exported from other GIS. Most of them are not so easy to draw.


texts:
I am trying to create a kmz using open source software. I have been testing the orfeo toolbox (https://www.orfeo-toolbox.org/CookBook/Applications/app_KmzExport.html) but it is just creating a flat dataset.

Anyone have tested this one?

Or any other solution?


texts: Measure distance between multiple ordered points

texts:
I'm using QGis 2.18.25 in Windows 10. I have a point vector with multiple points that represent a location in a certain date, so they have an integer field with the numbered order for each point as in the image. I would like to be able to get the distances between point 1 and point 2, point 2 and point 3, and so on. Using the measuring tool is not really an option since I have other hundreds of point vectors and it would take too long.

I've tried using a distance matrix but the end results are way different than when I manually measured them, and the Points2One plugin only generates one single line that bends in each point while I need the distance of each segment between them.

Is there any other way to do this?


texts: Why is my large DDP document losing symbology in later pages of PDF export?

texts: The map document is 77 pages long and it seems that later pages are losing some symbols for features. The scale isn't changing, so it is not an issue of scale based displays. Any thoughts?


texts: Addresses within radius of vector layer in QGIS

texts:
I want a CSV file of all addresses that intersect with historical stream features.
Currently, I have 3 datasets (see my blog) for this in QGIS 2.16.3 (MacOSX):

a layer (.gdb) of historical stream data (blue)
a layer (.shp) of address data (green)
a layer (.shp) of
property data (teal)
an OSM layer (grey) for reference only

I assume I need a query that does the following:

Select all properties that intersect with vectors from stream-layer
then from the resulting properties
select all addresses that intersect with those properties

Can this be done in a single Query (e.g. using ST_INTERSECT) or would I need something like PyQGIS to do this and would someone be able to help me with the syntax for that query/script?


texts: How to "merge" multiple raster ASTER L1T bands?

texts: I want to generate a Land Surface Temperature map using ASTER products. I used the Semi-Automatic Classification plugin and followed this tutorial. However my study area is too large and required at least 3 images to cover the entire study area. I want to "merge" them so I can clip them using the "Clip multiple raster" by the same plugin. I am sorry for using the word "merge" for I don't know the correct term on concatenating them, and it seems might require some algorithm to "dissolve" those overlapping bands.


texts: Create a fishnet inside an a symmetric polygon using ArcGIS desktop 10.3

texts:
I try to create a fishnet with 5 columns and 3 rows to this polygon:

but i get this result:

whereas is wanted the grid will be adjusted to the polygon direction like this :


texts:
I have a pandas dataframe containing a Shapely geometry object. The shapely object (circle) was computed from the following SO https://stackoverflow.com/questions/27431528/find-the-intersection-between-two-geographical-data-points, where I have lat, long, and radius as input.
import pandas as pd
from shapely.geometry import Point
from shapely.ops import transform
import pyproj
from functools import partial

WGS84 = pyproj.Proj('epsg:4326')

def latlonbuffer(lat, lon, radius_m):
    proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj4str) # azimuthal equidistant
    project = partial(pyproj.transform, AEQD, WGS84)
    return transform(project, Point(0, 0).buffer(radius_m))

For example, using the above function, I created a table like below, with an additional centroid column
geometry, centroid
POLYGON((26.48306 50.09625, 26.47916 50.09604, ..)), ((26.48307336330026, 50.052005610561245))

But when I convert this table to geopandas and plot it on mplleaflet
import geopandas as gpd
import mplleaflet
import matplotlib.pyplot as plt

df_gpd = gpd.GeoDataFrame(df)
ax = df_gpd.plot(figsize=(20,20), color='r')

mplleaflet.display(fig=ax.figure)

it displays the polygon in a different location than expected.
Below is what's plotted

And this is the location of the expected circle centroid

I figured out that even though the Polygon shows the coordinates in (lat, long) format, when mplleaflet is plotting, it swaps the coordinates to (long, lat).
Is Polygon object holding the coordinates in (long, lat) format, and mplleaflet is flipping the coordinates to properly plot it on the map?
The result I'm expecting is a circle centered around the centroid (26.48307336330026, 50.052005610561245) which is in (lat, long) format.
tldr: I would like to know if mplleaflet flips the coordinates to (long, lat) when plotting Polygon object in Geopandas, or the Polygon object using the function latlonbuffer() results in (long, lat) format.
Update:
The problem was with pyproj flipping the coordinates while transforming. We can keep the Polygon object to the GIS standard (long, lat) by adding always_xy=True in the transform parameter.
import pandas as pd
from shapely.geometry import Point
from shapely.ops import transform
import pyproj
from functools import partial

WGS84 = pyproj.Proj('epsg:4326')

def latlonbuffer(lat, lon, radius_m):
    proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj4str) # azimuthal equidistant
    project = partial(pyproj.transform, AEQD, WGS84, always_xy=True)
    return transform(project, Point(0, 0).buffer(radius_m))

Using this yields
POLYGON((50.09625 26.48306, 50.09604 26.47916, ..))

which can be directly plotted using mplleaflet without flipping the coordinates.


texts: What is the DBSCH command being passed to CartoDB's invalidation service?

texts:
CartoDB's PostgreSQL extension (cartodb-postgresql) has a function defined here:

https://github.com/CartoDB/cartodb-postgresql/blob/2e665a56b4a7dbc3bd8e3d91ac18b52bb932663c/scripts-available/CDB_GhostTables.sql#L2-L50

which includes the instantiation of a Redis python client (line 32). Then a call is made to that client (line 41) which looks like this:

client.execute_command('DBSCH', db_name, username, event_name)

DBSCH is not a Redis command, and I cannot find any documentation of it elsewhere in the CartoDB repositories.

Does anyone know where to find the source of that command?


texts: FME area on area overlayer dropping components of some features

texts:
As it says in the title, is this expected behaviour?

What happens is that the dataset before the area on area overlayer has some features (in the image below, the large green parts are the same feature, there is a point where they touch). The lower part of the green feature is lost when the transformer is run. This shouldn't happen!

I have an example dataset and workspace, here: https://www.dropbox.com/s/n50h39yftp7nhtx/fme_bug.zip.

There are a number of other cases where this happens in my dataset. Any thoughts?


texts: Modify a polygon that WKT is POLYGON EMPTY

texts:
First of all, I'm not a GIS professional, just a developer using SQL spatial tools.  I am inserting data in to ArcGIS feature class to be modified by users within ArcGIS desktop.  For most cases, there is an existing polygon that I can create from an existing ArcGIS layer as the starting point.  The user then within ArcGIS desktop edits the polygon by modifying the vertices.

This all works well, with the exception of when there is not an existing polygon to use as a base.  Initially I would insert a null as the shape into the feature class, but that caused issues in the application.  Then I found the concept of POLYGON EMPTY and that worked well within the application.  But my GIS friends are telling me that users can't modify this "empty" polygon.  They say can't modify what doesn't exist.

This doesn't make sense to me, can someone validate I am being told truths.


texts: Faster zonal stats for equal rasters with Python

texts:
I'm trying to extract data out of many raster images with different bounding boxes. Using the great rasterstats.zonal_stats it's a good approach, but it takes a considerable amount of time, especially for a large list of geometries. I wrote a small solution that works really fine for a set of raster files with equal bounds and resolution. Nonetheless, extending this to more complex scenarios, where rasters do not necessarily are equal (see image below)

Case A: raster images with equal bounds and resolution:

Looking at the source code behind rasterstats, the code iterates through each of the geometries passed to zonal_stats and use rasterio.features.rasterize, to mask the raster to the selected geometry and make the needed calculations. Despite being a robust approach it might be overacting when extracting > 100 rasters with similar metadata.

When a raster has the same metadata: same resolution, same boundaries
(bounding box), height, width, etc. It is possible to do the following:

1. Create a crosswalk dictionary:

Store all the indices for each of the geometries using one of the rasters sources (again, here I have a list of identical rasters). Hence, the indices will be equal across all files.

with rasterio('path/to/raster/.tif') as src:
    with fiona.open(/path/to/shape, 'r') as shp:
        geoms = [feature['geometry'] for feature in shp]
        index = [feature['properties'][id_var] for feature in shp]

        crosswalk_dict = {}
        for idx, geom in zip(index, geoms):
            geom_rasterize = rasterize([(geom, 1)],
                                       out_shape=src.shape,
                                       transform=src.transform,
                                       all_touched=True,
                                       fill=0,
                                       dtype='uint8')

             crosswalk_dict[idx] = np.where(geom_rasterize == 1)

2. Extract data

Loop through all the files and geometries in the crosswalk dictionary and calculate the mean using the source data for each of the files.

for path_raster in list_paths:
    with rasterio.open('/path/to/raster/.tif') as src:
     r_array = src.read(1)
     r_array[r_array == src.nodata] = np.nan

     mean_dict = {}
     for key, value in crosswalk_dict.items():
      mean_dict[key] = np.nanmean(r_array[value])

      df = pd.DataFrame.from_dict(mean_dict,
        orient='index',
        columns=['mean'])

This process outputs similar results that rasterstats.zonal_stats means and it is quite fast since the rasterization is done only once. Nonetheless, in the image shown above, images have no similar meta. Each of the images has its own affine and boundaries, although they have the same resolution.

Case B: raster images with different bounds:

In this case, I made a similar approach. All the functions in brackets are in this gist. I didn't write them here for space readability issues, they're quite verbose.

Calculate an envelope boundary (the blue box in the image) (draw_envelope)
Give all files the bounds from the envelope built in step 1, but keeping the same data (transform_raster_bounds)
Now with "equal" files, create a crosswalk list using either a meta dict, or a file to create the dictionary (get_crosswalk)
Extract!

This is fine -I believe-, but somehow the results are incorrect, or not comparable to the ones made by rasterstats. I know this is a huge question, but I'm out of ideas on what can be wrong with the transform_raster_bounds function. Also, is this an overkill? There's a way to optimize rasterstats?


texts: Normalize OSM geometry to 1x1 coordinate system

texts:
I'm extracting a couple of geometries with overpass turbo and need to "normalize" their shape from lat-lon to [[0-1][0-1]] coordinates.

The first thing I've tried is probably pretty wrong: get highest & smallest lat-lon values in the geometry and remap them to 0-1. This is the JS code:

function (geometry) {
    let biggestLat = geometry[0][0];
    let biggestLong = geometry[0][1];

    for (let i = 0; i < geometry.length; i++) {
        if (biggestLat < geometry[i][0]) {
            biggestLat = geometry[i][0];
        }

        if (biggestLong < geometry[i][1]) {
            biggestLong = geometry[i][1];
        }
    }

    let smallestLat = geometry[0][0];
    let smallestLong = geometry[0][1];

    for (let i = 0; i < geometry.length; i++) {
        if (smallestLat > geometry[i][0]) {
            smallestLat = geometry[i][0];
        }

        if (smallestLong > geometry[i][1]) {
            smallestLong = geometry[i][1];
        }
    }

    let coordinatesArray = [];

    for (let i = 0; i < geometry.length; i++) {

        let lat = this.map_range(geometry[i][0], smallestLat, biggestLat, 0, 1);
        let lon = this.map_range(geometry[i][1], smallestLong, biggestLong, 0, 1);
        coordinatesArray.push([lat, lon]);

    }

    return coordinatesArray
}

And the result is obviously pretty off:

I've been reading about the different projection methods and that I need to decide for one, but I don't really understand which I should use in this case. Mercator seems to be wrong?

To summarize my aim: I want to scale any Geometry i get from osm to the same size in a x,y coordinate system. Can anyone explain me how to do it correctly or what framework to use to do the calculation? I've been using turf.js in other projects.


texts: Resolving Java script Alert that This page was unable to display Google Maps element?

texts:
Whenever I try to load the Google Hybrid layer on QGIS 1.8.0 and 2.0.1, I get this error:

  This page was unable to display a Google Maps element. Please contact the site administrator. If you are the administrator of this site, please check the JavaScript console or check the following page for troubleshooting: http://g.co/mapsJSApiErrors

How do I correct the error so that the map can load?


texts: I am working with the latest stable version QGIS 3.4.9 Madeira. Unfortunately, I cannot find the window relating to projected coordinate system properties (see the picture), so I cannot change the coordinates of the central meridian. Can you help me, please?


texts: Label QGIS points layer with auto incrementing number count from 1

texts:
I have imported a CSV layer to my map and want to label the points with a numeric count from 1 upwards, counting just the visible points.

I have almost found a way to do this using the HouseNumbering plugin to automatically create a number count in a field (having converted the table from CSV to SHP file)

However I need to create maps of several different area views which each need to be numbered from 1, so will need multiple different fields with different numberings in the attribute table, and then swap labels depending on which view I'm in - which is messy

Is there a way to create an incremented count 'on the fly' without it needing to be in the attribute table?

The closest I can find is to label with $id which doesnt' work as the numbers don't start from 1 on my specific map view
Or $uuid which labels the points with a 38 characters unique string (not very attractive!)


texts: Proper labeling when multiple maps on one sheet of paper

texts:
Suppose I want to display 4 maps on one sheet of paper, all with the same scale and same direction for North. Assume they all have the same projection and datum.

Is it customary to place a unique scale bar and north symbol for each individual map? If not, how do you choose where to place it?


texts:
In QGIS, you can apply CMYK to text, using the CMYK sliders palette. However, I can only find an RGB colour picker for lines and fills.

Is there a way to apply CMYK to these elements?


texts: How to create a legend for quantitative visualization with increasing icon size?

texts:
I am trying to find a way to do a map and legend like this:

I am using QGIS 2.6 Brighton. The only way was to do a pie chart but I think it's not possible to create classes. I can create circles with a categorized style but I couldn't figure out how to create such a legend.


texts: GeoServer REST API (geoserver-restconfig) create coverage store and mosaic error

texts:
I am working on Ubuntu 18 and using Kartoza GeoServer docker v:2.15. I try to push Landsat band 4 (for sample I have two images epsg_4326) with Python library geoserver-restconfig to GeoServer with coverage store with below code

cat.create_coveragestore( 'name112', workspace=None, path='file:ndvi_4326/1.tif', type='GeoTIFF',create_layer=False, layer_name=None, source_name=None, upload_data=False,contet_type="image/tiff",overwrite=False)

when create layer is 'False' I can create store but if it is 'true' I get an error

  'FailedRequestError: Failed to create coverage/layer 1 for : name12, 400'

I got same error when I was trying with create_imagemosaic function.

cat.create_imagemosaic('band4','file:ndvi_4326',configure='all',workspace='ndvi')

~/Desktop/django_tutorial/djangoenv/lib/python3.6/site-packages/geoserver/catalog.py in create_imagemosaic(self, name, data, configure, workspace, overwrite, charset)
    532             resp = self.http_request(url, method='put', data=upload_data, headers=headers)
    533             if resp.status_code != 201:
--> 534                 raise FailedRequestError('Failed to create ImageMosaic {} : {}, {}'.format(name, resp.status_code, resp.text))
    535             self._cache.clear()
    536         finally:
   FailedRequestError: Failed to create ImageMosaic band22 : 400,

I think, when trying to create layer I get an error and I checked to code I couldn't find anything. Any suggestions for this problem?

Edit-1
I tried image mosaic code in Windows machine and it worked. What is the source of error? GeoServer setting? or Docker settings?

Also I opened github issue
https://github.com/GeoNode/geoserver-restconfig/issues/3


texts: Are multiple separate feature classes more efficient than repeatedly deleting all features from same feature class?

texts:
Quick Version: Is it more efficient, in terms of database performance, to have multiple separate temporary feature classes or one feature class that you repeatedly empty out and fill?

I am working with ArcGIS 10.1, Windows 7, using personal geodatabases.

I have a python toolbox tool that selects survey cells based on their proximity to a seed-cell.  Once a boundary is determined, based on distance from the last selected set of cells, the final selected set is buffered to form the boundary.  I have been using a single feature class (BoundaryHopper)to hold the output of the buffer and then appending that output to a second feature class that holds all the boundaries for a project (DerivedBoundaries).  BoundaryHopper would be emptied out at the beginning of each loop iteration to avoid repeatedly appending the same boundaries to DerivedBoundaries  ]

Recently I switched the tool over to using  separate outputs for each buffer output (TempBoundary1,2,3...) and appending each individual feature class to DerivedBoundaries.  That solved some issues I was having from trying to overwrite the same feature class all the time.

So, my question is: Are there any significant performance benefits of one technique over the other? Does one technique create less 'bloat' in the database or result in faster processing?


texts:
How do you rasterize a polygon layer into a raster so all cells covered by polygons (i.e. not NoData) will be set to a constant value?

I do this to set cell according to attribute values:

gdal.RasterizeLayer(rasterHdl, [1], featureLayer, options=["ATTRIBUTE=%s" % attributeName])

... but how about burning in a constant?


texts: Snapping an annotation feature class (anno class)

texts:
Is there any work around for being able to snap an annotation feature class (non-simple feature type) to the vertex points of another layer en masse?

Similar to how the 'Editing toolbox > snap' tool or 'Data Management > Feature Class > Integrate' tool works for 'simple' feature types.


texts: Does anyone knows if it is possible to manually modify SAGA GIS colour palettes, for example using ColorBrewer references?


texts:
I've finally been able to get two rasters with the same pixel size and perfect overlap (resample + snap raster). Now I need them to have the same extent, giving the value of "-9999" in pixels not present and eliminating those that are out of the raster used as a template.

As you see in the image what I need is to fit the "light green layer" to the "dark green layer".


texts: Transforming a .txt file into a raster in R

texts:
I have a .txt file which consists of 160 rows and 320 columns. It covers the geographic region from 34 to 41.95 in latitude and 19 to 34.95 in longitude. Its spatial resolution is 0.05 x 0.05. I use the following code to read the .txt file:

a1 = as.matrix(read.table("/home/...", header=F, sep = "\t",   as.is=TRUE))
rast1 = raster(a1)
extent(rast1) = c(34,41.95,19,34.95)
projection(rast1) = CRS("+proj=longlat +datum=WGS84")

Which gives me the following result:

class       : RasterLayer
dimensions  : 160, 320, 51200  (nrow, ncol, ncell)
resolution  : 0.02484375, 0.0996875  (x, y)
extent      : 34, 41.95, 19, 34.95  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : 0, 245.2738  (min, max)

The resolution given above is not correct. When I try to resolve this issue by running:

res(rast1) = c(0.05, 0.05)

I get the following result:

class       : RasterLayer
dimensions  : 319, 159, 50721  (nrow, ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : 34, 41.95, 19, 34.95  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

In that case the number of rows and columns in reversed and diminished by 1! Also, all the cells have NA's as values.

Could you please explain to me why this happens and how can I properly read my .txt file and transform it to raster?


texts: Irregular XYZ point data to raster (TIF) - LASTools alternative?

texts:
I am trying to convert my irregular spaced point data (xyz) to a DEM and then, to a raster (geotiff). Because of the irregular spacing I need a triangulation first, right (so, no spaces are in the final raster)?

So the way is: XYZ -> Triangulate -> Interpolate To Raster. I have a lot of points (minimum 150.000), my raster size should be 0.25 x 0.25m or smaller. I need a way to do this with open source / free tools!

What I've tried so far:

LASTools (las2dem) - works great, is extremely fast <- that would be exactly what I need, but is not free.
FUSION (TINSurfaceCreate then DTM2ASCII) - works well, but slows down extremely as the raster size shrinks (1 x 1m and 0.5 x 0.5m is okay, everything below is really slow).
SAGA Tools / gdal_grid (via QGIS) -> extremely slow as well.

Any other suggestions I should try?


texts: Leaflet issue: Marker and polygon not agreeing on locations?

texts:
I've got a leaflet map, to which I added a GeoJSON object (a polygon, specifically) and a number of Markers. The markers should all be inside the polygon, but they're not being displayed that way:

In the example below, the GeoJSON is:

{"type":"Polygon","coordinates":[[[29.53125,49.921875],[29.53125,50.625],[30.9375,50.625],[30.9375,49.921875],[29.53125,49.921875]]]}

And the green-highlighted marker's coordinates are:
[30.9332103729, 50.3615112305]

As you can see that marker is showing up well outside the outline of the polygon. The problem is considerably worse at wider zoom levels -- zoomed all the way in, the green marker makes its way just barely inside the polygon outline.

Any idea what could be causing this?

UPDATE:
Turns out the issue appears to be the Markers, not the GeoJSON... when I zoom in on the map, the markers move. In this example, they all move to the northwest as I zoom in.


texts:
ESRI recently put out a new tool called the 'Activity Dashboard' that checks for new maps created, map uses, etc. Sadly, this tool only works for those maps created through the ArcGIS online system.

One of the problems I have as the new 'Web Map' guy for my organization is that multiple maps have been created, all of them as (basically) one-offs to fix a specific problem, and no one has a master list of them or the services they are running on. This left me with a very annoying situation. I 'thought' I knew about them all, but I find that I am regularly learning about old maps as they break. The recent address locator that ESRI killed off brought SEVERAL maps to my attention that I wasn't even aware of.

I am wondering if there are any tools that have been developed that track maps and/or services that have been created to give someone in my position a better chance of managing this system.


texts:
I'm hoping to connect to a Esri WMS image layer through QGIS by querying/filtering the WMS so only specific layers are shown - specifically, trying to query out individual satellite scenes based on a project name. However I haven't been able to achieve this - whenever I try the whole WMS is drawn. Ultimately I want to create a bunch of XML's of individual projects to share.
I have constructed the URL as per the help - both encoded and not but it does not work. Example below -
<qgsWMSConnections version="1.0">
    <wms smoothPixmapTransform="false" ignoreGetFeatureInfoURI="false" name="ProjectName" ignoreGetMapURI="false" ignoreAxisOrientation="false" url="https://sat-image.com/arcgis/services/Satellite/Service_Name/WMSServer?Query=Name%3D%27ProjectName%27" invertAxisOrientation="false" username="" referer="" dpiMode="7" password=""/>
</qgsWMSConnections>

Using QGIS 3.1 and ArcGIS Server 10.7.


texts: How to assign binary value to points inside polygon

texts: I have three different polygon layers with some polygons that overlap. I have points that fall both inside and outside of these polygons.  I want to assign a value of 1 if a point falls inside any polygon and a value of 0 if a point is not inside any polygon.  How to do this?


texts: I am producing maps using OSGB references. However, I am using a 500m square grid, not 1km square. I would like to show the grid reference in a 3 digit format, the final digit being in a smaller font than the first two digits, e.g. the Easting 63 showing as 630, with the 0 showing as a smaller font, the next Easting is 635 with the 5 in a smaller font, the next Easting as 640 and so on. I am clueless on scripts. Sorry. I have searched on this site and on Youtube but can't find what I'm looking for in a format I can understand. I have found a piece of code I can utilise to produce a three digit format, left(right(@grid_number,5),2)but can't work out how to make the last digit smaller.


texts:
I am having problems with Python Shell importing the module "qgis.core".
When I type "import qgis.core" the Shell ("idle.pyw") gives me this error:

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import qgis.core
ImportError: DLL load failed: Impossibile trovare il modulo specificato (Trad "Impossible to find the specified module").

I have already set the environment to point the right folders following "PyQGIS cookbook" instructions.
In my case, the paths are:

PYTHOPATH=C:\"QGIS_path"\apps\qgis\python;
Path=C:\"QGIS_path"\apps\qgis\bin.

Am I missing something? What are supposed to be the right libraries the Shell is unable to find?

I am working with windows 7 64 bit, QGIS Lisboa (version 1.8), Python 2.7.2 (the one that is included in the QGIS package).


texts: How to import raster without any positional information to Qgis

texts:
I have PNG image exported from python script without any meta-data. I know coordinates of top-left corner and bottom-right corner and projection (I guess it's Mercator). How can I import this image to correct position?

To be sure, I know about Georeferencer in menu. However it seems like bad idea to resample original raster image. (I have to choose corner points manually, then enter target coordinates, and it's still not working.) Is there a way where I just enter 4 numbers and projection, and raster will appear on correct position in map? Or is it possible to somehow create metadata for this image so it will be imported automatically to correct position?

Edit: Ok, it is working when I manually set corners to correct pixel values and use Linear transformation. But it still seems like a bad solution. Adding some meta-data to image would be nicer.


texts:
I'm trying to come up with the Intensity value for the Cells dataset (Diggle 1983). This dataset has coordinates from 0 to 1, probably percent, there is no actual unit to it as I think it was taken from bugs or something.

What tool should I use to calculate the intensity of the layer? I looked at Point density but that produces a 'raster' and I am looking for just a single value.

I'm new to ArcGIS.


texts: Populating a field in an attribute table with a consecutive number sequence?

texts:
I know that this is probably a very basic question but I can't seem to find an answer through Google searches.

I would like to create a field in an attribute table (in ArcMap 10.2) for a Number Label (i.e.. 1, 2, 3, 4, ... etc.) is there a way to calculate this field with the field calculator rather than typing in each individual number?


texts: How to classify agriculture employment data?

texts: I have employment data which are huge - about 500 data for each category, like agriculture field, fishing field etc. Now I want to categorise this data in to 7 ranges. Suppose <50=1, >50-100=2, etc. for producing an agriculture map of a country by showing these categories as a legend. I don't want to change my main GIS data. I only want to show this data into small categories in Map. How can I do it?


texts: Layer style is not rendered properly using PyQGIS

texts:
I'm trying to render a layer with the given .qml style into .png image using a standalone PyQGIS application. The issue is that not all the layers of the symbology are rendered, e.g. a polygon with the colour fill, pattern fill and the svg-marker outline will be rendered with only colour fill, but pattern fill and svg-marker outline will be missing. And during execution of the script I get this (seems to be Qt-error) message: QPainter::begin: Paint device returned engine == 0, type: 3. Tried to google this error but found nothing helpful. QGIS 2.14.

The code:

from qgis.core import *
from PyQt4.QtCore import QSize
from PyQt4.QtGui import QPainter, QImage, QColor

qgisApp = QgsApplication([], False)
QgsApplication.setPrefixPath('/usr', True)
QgsApplication.initQgis()
print QgsApplication.showSettings()

lyr = QgsVectorLayer('bbox.shp', 'bbox', 'ogr')
lyr_id = lyr.id()
QgsMapLayerRegistry.instance().addMapLayer(lyr)
lyr.loadNamedStyle( 'bbox.qml')

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

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

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

layers = QgsMapLayerRegistry.instance().mapLayers()
lst = layers.keys()
print 'lst', lst
renderer = QgsMapRenderer()
renderer.setLayerSet(lst)

extent = renderer.fullExtent()
rect = QgsRectangle(extent)
# print rect.toString()
rect.scale(1.1)
renderer.setExtent(rect)

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

# do the rendering
renderer.render(p)
p.end()

# save image
img.save('test.png', 'png')

QgsApplication.exitQgis()

A test layer and style can be downloaded here.


texts: Disable GWC cache emptying after style changes

texts:
GWC is removing all generated tiles for layer after any edit of this layers style.

Is there any possibility to turn off this feature?

I've got few changes to make on working published geoserver so if I could disable cache cleaning after changes I could make it, and preview results via WMS and after the job is done manually empty and reseed cache - this will shorten downtime of service.


texts: Displaying list of features included in polygon

texts:
I want to find the object that is part of the polygon and displaying it in a separate listbox.

I have two layers with polygons, the first polygon (Red) is the polygon in which you want to search for features. The second polygon layer (Green) is the polygon that you want to display in the listbox when you enter the first polygon.

A yellow marked polygon that should not be affected because it is included in the first (red) polygon not completely I have little experience in this, I would like at least part of the implementation code of this At the same time I have to choose from the listbox layer with a polygon (Green)

 public void button1_Click(object sender, EventArgs e)
        {
Tool1 dbconnect = new Tool1();

 IFeatureLayer pfeaturelayer = (IFeatureLayer)mxDocument.ActiveView.FocusMap.Layer[0];
                    IDataset pDS = (IDataset)pfeaturelayer.FeatureClass;
                    WellManagerPolygon wl = new WellManagerPolygon(pDS.Workspace);
                    Well pWell = wl.GetWellByIDPolygon(Convert.ToString(textBox2.Text), Convert.ToString(comboBox1.Text));

                    int mRange = Convert.ToInt32(textBox1.Text);
                    // range of x meters
                    ITopologicalOperator pTopo = (ITopologicalOperator)pWell.PolygonLocation;

                    IWorkspaceEdit pWorkspaceEdit = (IWorkspaceEdit)pDS.Workspace;
                    pWorkspaceEdit.StartEditing(true);
                    pWorkspaceEdit.StartEditOperation();

                    IFeatureWorkspace pWorkspace = (IFeatureWorkspace)pWorkspaceEdit;

                    IFeatureClass pWellRangeFC = pWorkspace.OpenFeatureClass("WellRange");
                    IFeature pFeature = pWellRangeFC.CreateFeature();

                    pFeature.set_Value(pFeature.Fields.FindField("WELLIDRANGE"), Convert.ToInt32(textBox2.Text));
                    pFeature.set_Value(pFeature.Fields.FindField("RANGE"), Convert.ToInt32(textBox1.Text));

                    pFeature.Shape = BarsRange;
                    pFeature.Store();

                    mxDocument.UpdateContents();
                    mxDocument.ActiveView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeography, null, null);

                    pWorkspaceEdit.StopEditing(true);
                    pWorkspaceEdit.StopEditOperation();


texts: How to Render KML "Hatched" Polygon Symbology in Google Earth?

texts:
After solving the issue of KML polygons that were only partially filled-in in Google Earth (here: Why are These KML Polygons Not Filled With Assigned Colour in Google Earth?) I'm now wondering how to make Google Earth display the hatched colouration fill of certain polygons.

Here is how they look in ArcMap.

Here you can see that Google Earth simply renders the hatched polygons as solid black.

How can I render these KML Polygons in Google Earth with their true "Hatched" Symbology?  (or is it even possible?)


texts:
I've joined a table to a SQL-QueryLayer like posted here
Using IRelQueryTable

When I set
geoFeatureLayer.FeatureClass = relQueryTable as IFeatureClass;
the Layer in ArcMap forgot that it is an Query-FeatureClass. So in the Source-Tab of the Layer-Properties-Dialog there is no entry for the DataType and in the Fields-Tab the Fields are listed full-qualified in the Alias-Names. But the Join seems to work because it is correctly listed.
Any idea?


texts: Replace error message in ArcObjects .net custom geoprocessing tool

texts:
I am writing some validation in the updateparameters method of my custom geoprocessing tool. I need to post an error message when specific conditions exists and I would like to include a standard esri error code with my message. I have the logic worked out and my error message posting but I must not be using the error code correctly because no mater what number I use for the error code, it displays as -1.

I would like it to display error code 427.

This is the simplified code setting the error message.

 public void UpdateMessages(IArray paramvalues, IGPEnvironmentManager pEnvMgr, IGPMessages Messages)
 {
     Messages.ReplaceError(5, 427, "You're doing it wrong.");
     return;
 }


texts: Specifying which parameters to be used when using arcpy.ImportToolbox function?

texts:
This tool I have designed is to create a buffer around Gas or Power lines then append the newly buffered data into the planning area. My custom tool's script is as follows:

# Import System Modules
import arcpy
import os

#Set Workspace
from arcpy import env
env.workspace = "C:/Code/Module7FProj"
arcpy.env.overwriteOutput = True

# Sets Parameters for Gas and Power line buffer for geoprocessing tool
GasPower_Lines_Buffer = arcpy.GetParameterAsText(0)
Output_GasPower_Buffer = arcpy.GetParameterAsText(1)
GasPower_Buffer_Distance = arcpy.GetParameterAsText(2)
Planning_Area = arcpy.GetParameterAsText(3)
Append_GasPower_Lines = arcpy.GetParameterAsText(4)

# create a 25 foot buffer from power and gas lines
arcpy.Buffer_analysis(GasPower_Lines_Buffer, Output_GasPower_Buffer,     GasPower_Buffer_Distance)

# append buffered power and gas lines with planning file
arcpy.Append_management(Output_GasPower_Buffer, Planning_Area,     Append_GasPower_Lines)

I need to create a separate script which imports this tool, sets the multiple inputs for "GasPower_Lines_Buffer" and sets the input for "Planning_Area".

So far I have:

import arcpy

arcpy.ImportToolbox("C:/Code/Module7FProj/FinalProj.tbx")
# below is where I need to set the parameters from the
# previous script. Alias of toolbox is "FinalProjToolbx"
arcpy.FinalProjToolbx.FinalProj


texts: Mapinfo 10 with new look does not show toggle state of button

texts: Is this only a problem to me or a problem to other's too. I have mapbasic application that runs good on mapinfo 11 and mapinfo 10 with old look but this toggle button does not work with mapinfo 10 (The icons dll files too won't work) . What might be the problem?


texts: extract values from raster function taking too long


texts: Adding field in attribute table using ArcPy

texts:
I have researched How to add attribute field to existing Shapefile via Python without ArcGIS? but could not implement it.

Can give me examples of how to add field name in attribute table?

I have to show a hyperlink in the attribute table , using an online query.


texts:
In QGIS, I have plotted buffers around points. I would like to calculate the area which falls under buffers without over counting due to buffer overlaps.

Is there a quick and easy way to do this?


texts: GEE filter image to month before specified date

texts:
I am using a for loop in GEE to loop over several years of data adding bands to an image (MODIS fire data). I am currently adding a band for the mean NDVI in the month the fire occurred and the previous month. However, I am struggling to get it to work for adding it to the previous month. I have used ee.Algorithms.if() so that the previous month will be 12 if it is currently 1 (so goes from January to December) and that the year becomes the previous year if the current month is January.
Make collection of NDVI with mean monthly values:
// currently values for every 16 days, take the monthly average
var month_mean = ee.List.sequence(0, 20*12).map(function(n) { // .sequence: number of years from starting year to present (Jan 2001 to Jan 2020)
  var start = ee.Date('2001-01-01').advance(n, 'month'); // Starting date
  var end = start.advance(1, 'month'); // Step by each iteration

  return ee.ImageCollection("MODIS/006/MOD13A2").select("NDVI")
        .filterDate(start, end)
        .mean()
        .set('system:time_start', start.millis());
});
//print(month_mean);

var meanMODINDVI = ee.ImageCollection(month_mean);

Set parameters for filtering image collections:
var startMonth = 1;
var endMonth = 12;
var startYear = 2018;
var endYear = 2019;
var d = 1;
var m = 1;
var y = 2018;

 for (y; y <= endYear; y++){
for(m; m <=endMonth; m++){

//// filter fire data
var fire = ee.Image(fireCCI.filterDate(y+'-'+m+'-'+d).mosaic());
///// NDVI at time of burn
var modisNDVI = ee.Image(meanMODINDVI.filterDate(y+'-'+m+'-'+d).mosaic());
/// NDVI month previous
var j = ee.Algorithms.If(m==1, 12, m-1);
var k = ee.Algorithms.If(m==1, y-1, y);

var modisNDVIprev = ee.Image(meanMODINDVI.filterDate(k+'-'+j+'-'+d).mosaic());

var fireImageAddDate = fire.addBands(modisNDVI)
                            .addBands(modisNDVIprev)
                            .setMulti({"day":d,"month":m,"year":y});
var fireImageWithDate = ee.Image(fireImageAddDate);
}}

If I print j and k they will print the correct values, but as they are a computed object I get this error with trying to run the code:
NDVI previous month: Layer error: DateRange: Bad date/time 'ee.ComputedObject({
  "type": "Invocation",
  "arguments": {
    "condition": false,
    "trueCase": 2017,
    "falseCase": 2018
  },
  "functionName": "If"
})-ee.Number({
  "type": "Invocation",
  "arguments": {
    "condition": false,
    "trueCase": 12,
    "falseCase": 5
  },
  "functionName": "If"
})-1'.

I have tried casting the values using ee.Number() and ee.String() but still get the same error messages.


texts: Clip of polygons does not work properly by invalid geometry

texts:
I want to clip a polygon shapefile by another polygon shapefile. Input polygon is the green one (first), clip polygon is the light blue one (second) and the result should look like the dark blue one (third). But the real result looks like the red one: there are some parts of the clip polygon which are missing (dark blue color under the red polygon). The projections are the same.

When I check the vector validity (Vector - Geometry Tools - Check geometry validity) for clip polygon, there appear 32 errors; probable cause of this wrong result.

I use QGIS version 2.14.2-Essen. What should I do to get the right result (dark blue one)?


texts: PostgreSQL/PostGIS Geom to JSON with Inner Join to another table

texts:
I am currently producing JSON from a database successfully using the following query:

SELECT row_to_json(fc) FROM ( SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
FROM (SELECT 'Feature' As type, ST_AsGeoJSON(lg.geom)::json As geometry,
row_to_json((SELECT l FROM (SELECT extent, activationid) As l)) As properties  FROM hat.projectsgeom  As lg
WHERE activationid = 'HOT-0001' AND agency = 'HOT' AND type = 'activation') As f )  As fc;

However, I am stuck when I try and use INNER JOIN to another table.

This is my attempt so far ..

ERROR:  relation "activations.activationid" does not exist

LINE 5: INNER JOIN ( select * from hat.activations.activationid) zz ...

SELECT row_to_json(fc) FROM ( SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features
FROM (SELECT 'Feature' As type, ST_AsGeoJSON(lg.geom)::json As geometry,
row_to_json((SELECT l FROM (SELECT extent, activationid) As l)) As properties FROM hat.projectsgeom  As lg
WHERE activationid = 'HOT-0001' AND agency = 'HOT' AND type = 'activation') As f )  As fc
INNER JOIN ( select * from hat.activations.activationid) zz ON (hat.activations.activationid = hat.projectsgeom.activationid );

I have tried a number of variations, I think this is close, but I cannot get a result ..

What do I need to change in the second query to get the Inner Join to work properly and return all fields form the second table?

UPDATE:

Based on Jakub's response, I have been able to get closer.

SELECT row_to_json(fc) FROM
    (
        SELECT 'FeatureCollection' As type, array_to_json(array_agg(f)) As features FROM
        (
            SELECT 'Feature' As type, ST_AsGeoJSON(lg.geom)::json As geometry, row_to_json
            (
                (
                    SELECT l FROM
                    (
                        SELECT *
                        FROM  hat.projectsgeom
                        INNER JOIN hat.activations ON hat.projectsgeom.activationid = hat.activations.activationid
                        WHERE hat.projectsgeom.agency = 'HOT' AND hat. projectsgeom.type = 'activation'
                    ) As l
                )
            ) As properties  FROM hat.projectsgeom As lg

        ) As f
    )  As fc;

Is returning the following error

ERROR:  more than one row returned by a subquery used as an expression

I am baffled on how to get this to work, with the Json.

The table projectsgeom has multiple rows that need to be INNER JOINED to table activationid and all the columns from both tables returned as the properties, using the WHERE Clause on table projectsgeom

What Query will work?


texts: Latitude Longtitude to X, Y, Z coordinate system using python

texts:
does anyone know how to go about converting lat, long, alt to x, y, z where you have an origin for x/y coordinate (specific place in London)? In other words, I have a lat long and alt data and would like to convert the to x,y,z with respect to a specific place, not the Earth center.

For example, lat long is (51.5081° N, 0.0759° W) tower of london with altitude 100 meters and I want to convert my gps lat long to x y z considering the tower of london as the origin.


texts: Storing Get Count Tool Results in Data table using ArcGIS ModelBuilder?

texts:
I'm trying to take counts of feature class.

In ModelBuilder there is option Get Count tool it is useful, but the result will appear in result window it shows counts.

How do I store counts in a data table?


texts: Geoserver SLD default style for feature?

texts:
I'm new with GeoServer.

I'm trying to style a Geoserver layer. The features come from a PostGIS table.

The feature are all lines. I'd like to style them based on attributes (of the table). I can do that. I created rules that say "if myfield = value1 draw a red line, if myfield = value2 draw a blue line, ..." and so on.

I'd like to write a default style for the line: "if none of the above rules is true than draw the line yellow".

I didn't find in the Geoserver SLD documentation anything to accomplish that. Or my way of reasoning is different from that of SLD.

In the docs there's a UserStyle is Default property but I think it's not what I'm looking for.

I Googled too to no avail.

Any suggestion?


texts: Using featureCollection to define geometry in Google Earth Engine

texts:
I'm starting to use the GEE Python API and I'm struggling with the translation of some instructions from JavaScript to Python.
I'm using the following code in JS:
var all_alerts  = ee.ImageCollection('projects/glad/alert/UpdResult')
var country     = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', 'CG'));
var alerts_2020 = all_alerts.select('conf20').mosaic().clip(country);

Export.image.toDrive({
  image:alerts_2020,
  description:'alerts_rdc_2020',
  scale: 30,
  region:country,
  maxPixels: 1e10
})

and it work like a charm.
Naive I'd like to make the same with the python API:
def get_alerts(country_code, year):
    country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(ee.Filter.eq('country_co', country_code))
    all_alerts  = ee.ImageCollection('projects/glad/alert/UpdResult')
    alerts = all_alerts.select(year).mosaic().clip(country);

    file_name = 'alerts_' + country_code + '_' + year
    task_config = {
        'image':alerts,
        'description':file_name,
        'scale': 30,
        'region':country,
        'maxPixels': 1e10
    }

    task = ee.batch.Export.image.toDrive(**task_config)
    task.start()

country_code = 'CG'
year= 'conf20'
get_alerts(country_code, year)

The task start on my GEE taskboard but rise the following error :

Error: GeometryConstructors.LineString, argument 'coordinates': Invalid type. Expected type: List. Actual type: FeatureCollection.

which I don't understand because it the exact same call I was making (or so I think).


texts: Derive individual land cover classes from .jpg map

texts: I have a .jpg image which is land cover classification image with 9 classes. I want to extract raster image (.img file) which contains the same classes as .jpg image (9 classes too). Is there any tool to be used for this purpose in ARCGIS, ERDAS and ENVI?


texts: Finding center of point layer in ArcGIS Pro

texts:
I have a point layer that has Lat/long coordinates. How to find the center of a given point layer. For more understanding, I attached a drawing.
How is it possible to do this?

I tried to use the tools "Mean center", "FeatureToPoint", but they give false results and find exactly not the center, can there be other tools or can you somehow find it using python, I have an attribute table that contains only the lat/long coordinates (-180 :180) (-90: 90). Maybe you can use some formulas?


texts: Connecting QGIS to query on MS Access database?

texts:
Is there a method to connect to an MS Access database file (.mdb) using QGIS?

In particular I need to connect to a query in it.

The database does not have spatial data.


texts: Error in reading project with delimited text layer using PyQGIS?

texts:
I have multiple QGIS maps with a KML vector layer and a delimited text layer on top (csv file). I want to save these maps as images using a script as follows in the python console to do this.

for project_name in projects:
    QgsProject.instance().setFileName(project_name + '.qgs')
    QgsProject.instance().read()
    # Save content of map canvas as image
    qgis.utils.iface.mapCanvas().saveAsImage(project_name + '.png')

However, the projects are not being read with the error regarding the delimited text layer which is:

  File cannot be opened or delimiter parameters are not valid.

This is not the case when I open the project manually in the QGIS UI. The projects loads correctly and also gets saved as an image.

Edit: I upgraded from QGIS 2.8 to 2.18 and now it throws a dialog box which is a 'Handle Bad layers' one. I don't understand why that is when all the file paths are correct and the project opens normally outside the console.


texts: Displaying static and dynamic content on top of Map using ArcGIS API for JavaScript?

texts:
Using ArcGIS Javascript API, we have a custom widget to allow user to select his role to analyse the data on the Map. I must display the "Role : XXXX" on top right of the Map. This should be visible to the user always. "Role: " is static & "XXXX" is dynamic as and when user changes his role. I tried using TextSymbol to display the content on the Map. TextSymbol requires "Graphic point" to be added to Graphics Layer.

  var t = new TextSymbol("Role: ", font, color);
  var p = new Point([-100,52]);
  var g = new Graphic(p,t);
  graphicsLayer.add(g);
  map.addLayer(graphicsLayer);

By giving the point in graphics, the text moves along with the map as the point becomes fixed.

How do I add text to Map without map point?


texts: Connect to Esri map service using c++ in stand alone desktop app

texts:
Can anyone tell me how to connect to an Esri map service using c++ ?

Here is my code. The Open call returns E_INVALIDARG.

IAGSServerConnectionFactoryPtr AGSConnectionFactory;
AGSConnectionFactory.CreateInstance(CLSID_AGSServerConnectionFactory);
IPropertySetPtr connectionProps;
connectionProps.CreateInstance(CLSID_PropertySet);

VARIANT TestURL = _variant_t("http://server.arcgisonline.com/arcgis/services/");
hr = connectionProps->SetProperty(BSTR("connectiontype"), _variant_t(esriAGSConnectionTypeInternet));
hr = connectionProps->SetProperty(BSTR("url"), TestURL);

IAGSServerConnectionPtr AGSConnection;
AGSConnection.CreateInstance(CLSID_AGSServerConnection);
hr = AGSConnectionFactory->Open(connectionProps, 0, &AGSConnection);


texts: Calculate Statistics along a Linear Feature

texts:
I have a linear feature class (road network) that is split into segments less than 2 miles in length.  Each segment is assigned a value. I want to calculate a distance weighted sum in some manner such as the following:

Calculate midpoint of each segment
Search X distance units along the linear feature the segment belongs to and sum all values within the search distance with distance weighting (fractional multiplier of each value based on defined distance bins)
Return this weighted sum to each segment as a field value

I know how to do 1 and 3, but I'm stumped on how to accomplish 2. Usual methods for calculating spatial statistics ignore geometry and search in a full circle. I want to only search along the linear feature the segment is a part of.


texts:
Using MapInfo 12.5, I want to be able to print a property at a set scale - say 1:1000 or 1:5000.

What is the easiest way to do it?

I want to be able to show a Scale Bar so that the user can simply measure length on the property with the scale bar.


texts: Looking for GIS dataset for Japan with attributes in English

texts: I have found some GIS data for Japan but the location attributes (prefecture, city names etc) are in Japanese. Anyone know where I could find this data in English, I could get it translated but it would take a long long time.


texts: Is there a way to undo the curve while editing or undo a point where you accidentally placed it? Sometimes I accidentally place a point where I didn't mean to while in edit mode for "add circular string" and I can't undo the point. The only way is to end the point and delete the curve. In ArcGIS you can click a point and undo immediately but I can't seem to find a way to do this in QGIS.


texts: Adding data at runtime into dataframes using ArcPy?

texts:
I am trying to create an empty MXD to act as a template.

I will have three dataframes for maps, each holding a different object from ArcSDE (1 raster object and 2 features).

I know how to query all of the data frames, but I am having trouble finding out how to add data to the dataframes.

Is there an online guide anywhere?


texts: Can GIS software be used to open .ulf files?

texts:
I was sent a .ulf file and I have been trying to open it using various software.

I read at http://fileinfo.com/extension/ulf that this file can be open using gis software.

  Data file formatted in Universal List Format (ULF); contains a list of
  elements similar to an .XML document; often used by Geographic
  Information Systems (GIS) applications to store coordinates and
  vertexes.

Can anyone confirm this and if so, how do I open the file using this software?


texts: Smoothing raster with binary cell values

texts:
Main question towards the end. Details follow now:

I have a slope raster

Region:            Bern, Switzerland
CRS:               EPSG:2056 - CH1903+ / LV95 - Projected
Extent:            2587757,1186900 : 2613583,1212068
Map units:         meters
Cell values:       degrees
Theoretical range: 0-90
Range:             0-81
Pixel size:        2x2m
File size:       ~ 500 MB

I ran the file through raster calc with the expression

( "filename" < 2 )

to obtain all cells with a slope of less than 2. I have now a raster with cell values of 0 and 1. Some areas of the raster are pretty homogenous, but some areas are veritable checkerboards where the slope of the original file is just around 2 and now alternates almost from cell to cell between 0 and 1.

I need to polygonize the raster. So I tried the GDAL function polygonize from the QGIS toolbox, but had to kill the task after running for 25 minutes.

How can I generalize or simplify a noisy binary layer to reduce processing time when polygonizing?


texts:
I'm quite a newbie to QGIS and GRASS.. I have installed QGIS 2.10.1 with GRASS 6.4.3 on my Windows 7 laptop and wish to open GRASS layers in QGIS. If I got it right, you can use the QGIS built-in browser for this. So, for example, I wish to see the "elements" vector layer:

I try to drag the folder (created with GRASS topology tools) to the Layers window. But nothing happens, the folder doesn't move with the mouse. Should it move, or is this perhaps some Windows issue? Or am I missing some very elementary skill with QGIS?


texts:
I need to use QGIS Web-client to publish a research project on web.
I setup the software on Windows and on Ubuntu Linux, changing the configuration as described in the README file and replacing the qgis_mapserv.fcgi with the server for the correct platform.

Now I have the client on Windows which isn't working at all. On Linux I can load the helloworld.qgs example, but I see the dreaded pink screen.
The 'natural' map is not loading, it keep going forever and without no one connection with the server.

The QGIS server is working correctly on Windows and not really in a good way on Linux, but it doesn't matter. I need a working configuration for one system, OS independent.


texts: How to find out what features /tags are included in an .osm or .pbf file?

texts:
I would like to load osm data into a postgis (postgresql) database. In
How to load OSM data into PostGIS preserving all Tags? there has been a solution posted:

Download the OSM files (in .osm format) and import using osm2pgsql.

When using osm2pgsql, there's a file called default.style that allows you to control
which OSM feature goes into which postgresql table, e.g.:

way        boundary     text         polygon
node,way   railway      text         linear

means that

    features tagged "boundary" go into the *planet_osm_polygon* table
    features tagged "railway" go into the *planet_osm_point* and *planet_osm_line*
        because the "railway" tag is common between railway stations (osm nodes) and
        railway tracks (osm ways).

Now I wonder how to find out which features are contained in the .osm or .pbf file to configure  "default.style"

Any input on how to make a good configuration file (default.style) would be appreciated. Thanks!


texts: Looping through multiple feature datasets and Intersect Files?

texts:
I am getting an error with a python script I have created.  I have 1 geodatabase with 3 feature datasets.  What I would like to do, is loop through 1 feature dataset (NB) to get a list of feature classes.  One by one, I would like to get the name of the feature class (ID_NB) and search for its corresponding feature class in another dataset (called CCH) so that I can perform the Intersect function on the two files and save the output in the 3rd feature dataset and call it ID_Ints.  **There may not always be a matching feature dataset between the two either.

The script works fine the first time through, but when it comes to the second feature class, it gives me an error saying "Dataset does not exist or is not supported".  I have a feeling it is because I have switched workspaces in a for loop, but I can't think of another way around this.

Anyone have any ideas?

#Define Environments
CCHws = "C://TEST.gdb//CCH"
PerMap = "C://TEST.gdb//NB"
opws = "C://TEST.gdb//opws"

#Set workspace
arcpy.env.workspace = PerMap
arcpy.env.overwriteOutput = True
NBList = arcpy.ListFeatureClasses()

#Start forloop, find first NB Feature Class and extract the ID (the first 5 characters)
for NB in NBList:
    NBname = NB[:5]
    NBname2 = NBname + "*"
    NBnameEx2 = "\"" + NBname2 + "\""

#Calculate Original Area of NB feature Class
    arcpy.AddField_management(NB, "NB_Area", "DOUBLE")
    arcpy.CalculateField_management(NB,"NB_Area","!shape.area@meters!","PYTHON")

#Set variables to use in Intersect
    CCH = str(CCHws) + "//" + str(NBname) + "_CCH"
    opresult = opws + "//" + str(NBname) + "_Ints"
    NBfile = str(PerMap) + "//" + str(NBname) + "_NB"

#Re-set workspace to look for matching ID with an '_CCH' in the CCH dataset
    arcpy.env.workspace = CCHws
    arcpy.env.overwriteOutput = True
    cchList = arcpy.ListFeatureClasses(NBname2)
    for cch in cchList:
        if len(cchList) > 0:
            arcpy.Intersect_analysis([NBfile,CCH],opresult)
        else:
            print NBname + " Does not have a matching CCH"


texts: Spatial clustering points on a spatial network in QGIS / PostGIS

texts:
I have some thousands of houses/addresses (points) and a road network and I would like to build clusters of defined minimum to maximum cluster size not using euclidean distance as a metric, but the distance between the addresses along the road network (where I'm not sure if I'm using the word 'cluster' in the correct context here).

I read this paper http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.6037&rep=rep1&type=pdf and had a look at PostGIS function ST_ClusterDBSCAN, which is unsuitable for my task because DBSCAN algorithm is not intented to limit cluster size nor is it possible to define a non-euclidean metric for the input data.

K-Means approach might not be suitable, because the numbers of clusters is not defined a-priori due to varying count of addresses within a cluster.

Everything seems to point me to hierachical clustering (I had alook at How to cluster points into clusters of a maximum diameter in PostGIS?), but by now I could not find a solution taking the metric along a spatial network into account.

For better understanding ref. to https://youtu.be/iQvpPl-J7b4?t=328 , the guys do something similar to what I'm after.

Are there any open solutions to this?


texts: Contour Lines exported as raster in Illustrator

texts: I am trying to export contour lines from arcgis to illustrator to finish the map. I have tried to export the map with "vectorize layers with bitmap markers/fills" option, but it still opens as a raster in illustrator


texts: how to create a grid of points within an irregular polygon with the maximum number of points?

texts:
I usually use qgis and arcgis.
I need to create a grid of points within a polygon with the maximum number of points. I also need to set the grid parameters: x distance, y distance, offset from the border of the polygon.
I would like to develop grids as in the examples (the first example is a grid where the rows are equal two by two; the second example is a grid where the y axis is inclined)
I need this grid to fit inside rectangle or rhombus(image3).


texts: Specifying band combinations when loading LANDSAT MTL files in ENVI IDL?

texts:
I'm somewhat new using ENVI and have a very limited background in programming. I was asked to compile 12 RGB and NIR composites of LANDSAT scenes within a specific timeframe and would like to use ENVI's IDL to not only create the composites with script but batch process each rather than do the 24 composites manually.

I have been looking around everywhere and the only somewhat helpful page has been https://www.harrisgeospatial.com/docs/ENVILayerStackRaster.html. However, I don't know how to call on specific bands from the MTL file to stack like I normally would when loading an MTL file into ENVI. All I can find is how to call on the OLI or TIRS bands as a group (https://www.harrisgeospatial.com/docs/MetaspectralDatasets.html).

Once I figure that out then I'm onto the task of figuring out how to batch the process.
A little guidance would be helpful. I know I'm rushing into this without the proper foundation to be using IDL.


texts: How do I calculate the percentage of progress using raster data using Python?

texts: How do I calculate the progress amount of two or more raster data being processed? I have seen examples of vector data but it becomes difficult to calculate the progress when no vector data are involved. I have only seen solutions for vector data here it easy by looping around the features in a vector file


texts: Applying Canadian Postal Code areas in shapefile to Google Earth?

texts:
So for someone who just downloaded and opened Google Earth and then searched for how to view postal codes on google earth this is all very confusing.. I downloaded the Crowdsourced Canadian Postal Code Polygon File Shapefile format As of September 22, 2016. but I don't know what to do with it.

Can someone help me?


texts: Exporting attributes from QGIS to .dxf/dwg

texts: I have a QGIS point layer with attributes in table. I want to export a layer with some of the attributes as a CAD block with attributes. Is it possible to do something like this?


texts: PyQGIS: Inline function to replace NULL values with 0 in all fields (not coalesce)?

texts:
Now I know there are several ways to tackle this problem but so far none of the following worked for me.
I got a huge attribute table out of various spatial joins so pretty much all of the fields have NULL values in some rows. This table will be processed and new fields are going to be calculated programatically so the solution involving the coalesce function is out of the race because the expressions for the following calculations are complex enough and I don't want to do all that typing. So here are my failed attempts of fiddeling together a function that i can use inline:

def replaceNulls(layer):
   context = QgsExpressionContext()
   context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))
   with edit(layer):
      for field in layer.fields().names():
          for feat in layer.getFeatures():
              expr = QgsExpression('if({f} is null,0,{f})'.format(f=field))
              feat[field] = expr.evaluate(context)
              layer.updateFeature(feat)

Using an expression from this answer replaces ALL the fields with 0 and some with field - independent of their field type.

def replaceNulls(layer):
   context = QgsExpressionContext()
   context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))
   with edit(layer):
      for field in layer.fields().names():
          for feat in layer.getFeatures():
              expr = QgsExpression('case when {f} is null then 0 else {f} end'.format(f=field))
              feat[field] = expr.evaluate(context)
              layer.updateFeature(feat)

This function replaces all values that are NOT NULL with 0 for some reason all though it is said here that:

  The expression syntax is exactly the same as in the field calculator GUI. So you can test them in the GIU first, before copying them into your PyQGIS script.

If that is the case what am I missing out on? The same source also mentions a more pythonic way but how to if feat[field] is None: when the NULL value is still a QVariant?

TL;DR: replace(NULL,0). How?


texts:
I have two National Imagery Transmission Format Standard (NITF) files loaded into QGIS 3.10 and they should be the same coordinate system but for some reason I cannot figure why they are not aligned.

Steps to recreate:

1) Loaded the following NITFs files into QGIS. One image is synthetic aperture radar and the other is electro-optical imagery.

2) Layers do not align:

3) I loaded both files into GDAL in python and printed out the result of GetProjection() and got the same result for both:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

Additional info:

Dataset is from this data source: https://github.com/AFRL-RY/data-unicorn-2008

Edit:

This is the .prj file I get when use QGIS after loading the images (Raster>Projections>Extract Project) rather than use python GDAL's getProjection()

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

I tried changing to the different CRS's as suggested the misaligned SAR image either stayed the same or disappeared (probably incorrectly projected somewhere else on the globe). Some of the EPSGs were not present in the options in QGIS.

I set a basemap and we see that the EO image is aligned but the SAR seems rotated.


texts: Advice on learning Remote Sensing software?

texts: Can you please advice me on which Remote Sensing software ? Which one is popular ? Which one is the least ? Which one is the easy to learn curve ? Pros and Cons welcome ?


texts: Adding barriers to a Network Analyst layer programmatically

texts:
Right now, I have code that allows the user to indicate on a map, which contains a Network Dataset layer, stops that are needed to create a route.  I need to add functionality to the interface that allows the user to also create barriers as well as stops.

The way this is normally done in ArcMap is as follows:

Create the route layer by clicking on Network Analyst and New
Route
Make the Network Analyst window visible
Highlight the Restriction layer inside the Line Barriers layer
Click on the button for Create Network Location Tool
Draw the barrier on the map.

I already have implemented the ability to create stops by using ControlsNetworkAnalystCreateLocationToolClass().

My question is, how can I also implement the capability of creating barriers as well as stops into my code?


texts: Editing layers in GeoServer with web application

texts: I wanted to edit layers in GeoServer where they are store in PostGIS database. I have refereed http://www.gistutor.com/openlayers/22-advanced-openlayers-tutorials/47-openlayers-wfs-t-using-a-geoserver-hosted-postgis-layer.html link for that, which works. But now i also wanted to edit previous feature and also add attributes of feature that has been added (note: from this link edited feature cannot be edited again once web application is refresh and attributes cannot be added) . I also wanted to correct edited mistake without entirely removing the feature. I also refereed links like Ref Link 1 and Ref Link 2. Please help me out for same. From this i wanted to create an web base application for editing layers from GeoServer . If any other opensource software if you can suggest for this web application please suggest. Please help me out for same.


texts: Label in the middle of linestring in OpenLayers

texts:
Is there anyone can put label in the middle of the linestring dynamically? I've used labelAlign option. I've changed this value from bc, cc, cm, etc. But it doesn't work at all. There's also option for x offset and y offset. But I don't know how to calculate the right position so the label can show at the middle of linestring. Not at the first line string.

Do you have any idea how to do it?

Thanks,


texts: Showing Linestring from SQL statement with CartoDB

texts:
I have cartodb sql statement to return what polygons are intersected by a linestring, but I'd like to have the linestring visualized over the polygons on the results. Also not add the linestring to the database.

SELECT country, permit_inf, permit_req, the_geom_webmercator
FROM airspaceshape_1
WHERE
  ST_Intersects(
    the_geom,
    ST_SetSRID(
      ST_GeomFromGeoJSON('{ "type": "LineString","coordinates": [ [27.456667,5.065], [30.456667,-1.088] ] }'),4326)
  )


texts: What is expression for finding perimeter in ArcPy and Field Calculator?

texts:
What is the correct format for finding perimeter as an expression for arcpy.
I was able to use expression !shape.Area! to find the area of my polygon.
However when I use the expression !shape.Perimeter! I get invalid expression.
I've tried the following expressions;

!shape.Perimeter!@KILOMETERS!

!shape.geodesicPerimeter!

!shape.geodesicPerimeter@KILOMETERS!

These do not work. I was able to find the perimeter using the attribute table and right clicking and going to calculate geometry. Therefore I know its possible, just I require the syntax for my python code.

My python code that works for area looks like this.

arcpy.CalculateField_management("Iceberg","Area","!shape.geodesicArea@SQUAREKILOMETERS!","PYTHON_9.3","#")

However my python code for perimeter is not working and it looks like this;

arcpy.CalculateField_management("Iceberg","Area","!shape.geodesicperimeter@KILOMETERS!","PYTHON_9.3","#")

I must be missing something, and I've looked at some examples, and have also done, X and Y coordinates, all work, but I cannot get my perimeter to work :(


texts: Blocked loading mixed active content error when loading named map as layer in CARTO?

texts:
I am unable to get my layers to draw on my Carto map using https. My map is created I can see the basemap, but the layers do not draw. The layers load locally, but they won't when I host the site on my server. Instead, I get the error "Blocked loading mixed active content: “http://[user].carto.com/api/....":

I have followed the instructions at this answer  and CARTO https-support

1) I used https when importing both the JS library and the CSS files:

<link rel="stylesheet" href="https://cartodb-libs.global.ssl.fastly.net/cartodb.js/v3/3.15/themes/css/cartodb.css" />
<script src="https://cartodb-libs.global.ssl.fastly.net/cartodb.js/v3/3.15/cartodb.js"></script>

2) I am using createLayer, so I set the { https: true } parameter. It looks like this:

cartodb.createLayer(map, {
    user_name: '[user]',
    type: 'namedmap',
    https: true,
    named_map: {
        name: "linessimple",
        layers: [{
            layer_name: "t"
        }]
    }
})
.addTo(map);

I even tried adding the alert that is suggested in the help documentation. The alert says the layer was loaded, but then the browser blocks it.  Tested in FireFox, Chrome, and Internet Explorer. Same problem.