texts: Rescale Raster Band With GDAL Python Bindings

texts:
I am trying to figure out how to rescale a raster band's pixel value using the GDAL python bindings.  I am looking to emulate the -scale switch in gdal_translate.exe; rescaling the pixel values in my raster to 0-255.

EDIT:  Some goofing around today has produced this, seems to work thus far.

from osgeo import gdal
import numpy
import math

def singleBandToRGB(source,output):
    src_ds = gdal.Open(source)
    band = src_ds.GetRasterBand(1)
    min = band.GetMinimum()
    max = band.GetMaximum()
    data = src_ds.ReadAsArray(0,0,src_ds.RasterXSize,src_ds.RasterYSize)
    data = ((255-0)*((data-min)/(max-min)))+0
    data = data.astype(numpy.uint8)
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(output,src_ds.RasterXSize,src_ds.RasterYSize,3,gdal.GDT_Byte)
    for i in range(1,4):
        dst_ds.GetRasterBand(i).WriteArray(data)
        dst_ds.GetRasterBand(i).ComputeStatistics(False)
    dst_ds.SetProjection(src_ds.GetProjection())
    dst_ds.SetGeoTransform(src_ds.GetGeoTransform())

    # housekeeping
    driver = None
    dst_ds = None
    src_ds = None


texts:
I would like to use lidR to produce canopy height rasters. I was following along in the 'rasterizing perfect CHM' user guide and now I'm wondering how to export this raster to a GeoTiff or .img. In particular, I think I'm going to loop through each tile (of thousands) and produce a CHM using this line of code...

chm6 = grid_canopy(las, 0.25, pitfree(c(0,2,5,10,15), c(0,1), subcircle = 0.2))

Of course, I'll change the pixel size, etc. for my use case.


texts: Reproject full swath data with VRT files

texts:
I am attempting to follow the following link to reproject swath data to a grid. http://geoinformaticstutorial.blogspot.com/2014/06/reading-and-map-projecting-amsr-2-data.html

In this link it is suggested to implement .vrt files which link to other vrt files detailing the latitude and longitude arrays of interest which are stored in netcdf files.

My specific files I have created are as follows and are approximated for the purpose of this question:

lat_Copy.vrt

NETCDF:"C:\Path\to\latlon\netcdf\file\file.nc"://lat_array

lon_Copy.vrt

NETCDF:"C:\Path\to\latlon\netcdf\file\file.nc"://lon_array

data_Copy

"C:\Path\to\lat_Copy.vrt"
"C:\Path\to\lon_Copy.vrt"
NETCDF:"C:\Path\to\data\netcdf\file\data.nc"://data_array

And then I run gdal with the following command:

gdalwarp -geoloc -t_srs EPSG:4326 "C:\path\to\data_Copy.vrt" test01.tif

Unfortunately the following error occurs when running this command:

ERROR 4: `C:\path\to\data_Copy.vrt' not recognized as a supported file format.

I am not sure where I have gone wrong in the file format. I have attempted adding in the 1's and 0's after each line as suggested in the blog post but nothing has come of it and only the same error appears.

I do have alternate .vrt files for the lat lon datasets that I have made but still produce the same error as above but I will post them if requested.

EDIT: I was able to solve my conundrum by copying the file format at this link in the answer provided by fer.marino Map project a raster having separate latitude and longitude raster bands


texts: GPX-Tracks are not displayed correctly by many applications/web applications

texts:
I have a problem regarding GPX tracks that were created from Shapefiles using QGIS.

The shapefiles were created by selecting all parts of a route from a PostGIS database.
The polylines were merged with QGIS and then exported as a GPX file by using the 'Save as' option.

The resulting GPX files look like this:

<?xml version="1.0"?>
<gpx version="1.1" creator="GDAL 1.11.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogr="http://osgeo.org/gdal" xmlns="http://www.topografix.com/GPX/1/1" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
<metadata><bounds minlat="51.533449555703989" minlon="6.983701320292362" maxlat="51.554475633460669" maxlon="7.010408625662222"/></metadata>

<trk>
<extensions>
<ogr:OBJECTID>31751.0</ogr:OBJECTID>
</extensions>
<trkseg>
<trkpt lat="51.54792023947288" lon="6.986589223144231">
</trkpt>
<trkpt lat="51.547997675144856" lon="6.986820893755352">
</trkpt>
<trkpt lat="51.548138136988818" lon="6.986910766789045">
</trkpt>
<trkpt lat="51.548352505219341" lon="6.987013118209317">
</trkpt>
<trkpt lat="51.549051234512156" lon="6.987331861261017">
</trkpt>
</trkseg>
<trkseg>
<trkpt lat="51.54792023947288" lon="6.986589223144231">
</trkpt>
<trkpt lat="51.547840418735099" lon="6.986046522387352">
</trkpt>
<trkpt lat="51.547799254904852" lon="6.985851100932796">
</trkpt>
<trkpt lat="51.547750380482739" lon="6.985724171010038">
</trkpt>

And so on. There are 8 segments of the track in total.

The issue is that while some web applications are able to display this track without a problem (like http://www.gpsvisualizer.com/)

many others are not (http://www.maplorer.com/view_gpx.html )

I have already tried adding timestamps to the points, but that didn't change anything.

I assume that the issue has something to do with the segments because if I merge all segments of the track using the GPX_Editor from Sourceforge the route gets distorted exactly like it does in the maplorer example.
No matter what application i use, the track is either displayed correctly as in the first example or falsely like in the second example.


texts: Basic instructions for loading OS Terrain 50 data

texts:
Are there any basic instructions for using the OS Terrain data available online?

I do not have an issue extracting - it is all sat in a folder ready to go.

What I cannot work out how to do is to load the data in QGIS and then to change the settings so that I can use it to do a viewshed analysis.

Please note that I don't know how to use scripts or anything complicated! I'm using Wein on a Windows PC and have downloaded the data from OS as shapefiles and GML.


texts: Converting Smallworld data into other GIS formats?

texts:
Is there a tool for converting Smallworld data into a variety of other GIS formats?

We have a number of Smallworld instances and need to create a central repository for the data that can be used by a mobile application.


texts:
I have a SHP file which have UTF8-encoded char fields. It also has a corresponding CPG file with the value 65001. The language byte (#29) in the DBF is currently 0.

How do I get it to keep the encoding in the output, when I use ogr2ogr to select a subset of the records? Output is also SHP. At the moment it changes it into $57 = codepage 1252.


texts:
In Geoserver I have a layer based on SQL view. I enabled WMS direct integration in GeoWebcache and I enabled caching in layer.

In Google Maps javascript API I add a layer in this way. Layer is projected, but caching not happens:

var wmsLayer =
new google.maps.ImageMapType({
getTileUrl: function (coord, zoom) {

    var s = Math.pow(2, zoom);
    var twidth = 256;
    var theight = 256;

    //define boundingbox
    var gBl = cmap.getProjection().fromPointToLatLng(
        new google.maps.Point(coord.x * twidth / s, (coord.y + 1) * theight / s)); //zuidwest coördinaat
    var gTr = cmap.getProjection().fromPointToLatLng(
        new google.maps.Point((coord.x + 1) * twidth / s, coord.y * theight / s)); //noordoost coördinaat

    var bbox = gBl.lng() + "," + gBl.lat() + "," + gTr.lng() + "," + gTr.lat();

    //basics WMS URL
    var url = "http://geoserver.xxx.be/geoserver/rcs/wms?";

    url += "&service=WMS";
    url += "&version=1.1.0";
    url += "&request=GetMap";
    url += "&layers=rcs:rcs-sql-server-source-with-parameters-public";
    url += "&styles=";
    url += "&format=image/png";
    url += "&TRANSPARENT=TRUE";
    url += "&srs=EPSG:4326";
    url += "&bbox=" + bbox;
    url += "&width=256";
    url += "&height=256";
    url += "&customheader=" + new Date().getTime();
    url += "&viewparams=kind:bike;groupid:1";
    url += "&tiled=true";

    return url;
},

tileSize: new google.maps.Size(256, 256),
opacity: 0.85,
isPng: true
});


texts:
I am always after improving my mapinfo outputs and I wondered what tips you have to improve the general look of an output.

I suppose this comes from having a good selection of both raster or vector layers.

I have the OS layers. Any tips?


texts: Using original Landsat multispectral image or pan-sharpened image in land cover classification?

texts: I am currently creating land cover maps of a number of national parks across Asia from Landsat 8 imagery and also later from SPOT imagery. I plan to use supervised classifications in ERDAS IMAGINE and we hope to create accurate maps showing individual types of vegetation. I am wondering if it would be better to perform this classification on the original multispectral image or the pan-sharpened image. Obviously, it would be preferential to use the pan-sharpened image as the spatial resolution would be higher, but I'm afraid that there may be some loss of spectral information in the resampling, especially when computing NDVI values. I am not sure whether this loss of information is worth compromising for the increased spatial resolution. I am fairly new to using this technique so just hoped someone could give me some advice!


texts: Allow null values from Value Table filter list in Python toolbox parameter

texts:
I'm setting up parameters in a Python Toolbox and have a Value Table with three columns. Two of the columns have value lists for the user to pick from. I would like for the user to be able to choose one, but leave the other Null. However, the tool will throw and error message that Null is not part of the value list and the tool will not run. If I try adding Null to the value list, the parameters break and do not show when the tool is opened. Is there a way to a user to pick a value from one column in a value table parameter, but leave the another blank?

Here is my code for my value table parameter:

target_elements = arcpy.Parameter(
    displayName = "Target Element(s) Name",
    name = "element_name",
    datatype = "GPValueTable",
    parameterType = "Optional",
    multiValue = True,
    direction = "Input")
target_elements.columns = [['GPString','Target Element Type'],['GPString','Target Element'],['GPString','EO ID (if known)']]
target_elements.filters[0].list = ["Insect","Lepidoptera","Other Invertebrate","Plant","Vertebrate Animal","Community"]
target_elements.filters[1].list = et_list


texts: Error: The requested Style can not be used with this layer

texts:
I'm trying to calculate the area of a polygon using the code from this page and I'm getting the following error

  The requested Style can not be used with this layer.  The style specifies an attribute of geometry and the layer is workspace:layer

Here's my code:

<?xml version="1.0" encoding="ISO-8859-1"?>
<sld:StyledLayerDescriptor version="1.0.0"
xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
    xmlns="http://www.opengis.net/sld"
    xmlns:sld="http://www.opengis.net/sld"
    xmlns:ogc="http://www.opengis.net/ogc"
    xmlns:gml="http://www.opengis.net/ogc"
    xmlns:xlink="http://www.w3.org/1999/xlink"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
    <sld:NamedLayer>
        <sld:Name>style_osm_polygon</sld:Name>
        <sld:UserStyle>
            <sld:Title>'Polygon'</sld:Title>

            <sld:FeatureTypeStyle>
                <sld:Rule>
                    <ogc:Filter>
                        <ogc:And>
                            <ogc:PropertyIsEqualTo>
                                <ogc:PropertyName>natural</ogc:PropertyName>
                                <ogc:Literal>water</ogc:Literal>
                            </ogc:PropertyIsEqualTo>
                            <ogc:PropertyIsGreaterThan>
                                <ogc:Function name="area">
                                    <ogc:PropertyName>geometry</ogc:PropertyName>
                                </ogc:Function>
                                <ogc:Literal>1000000</ogc:Literal>
                            </ogc:PropertyIsGreaterThan>
                        </ogc:And>
                    </ogc:Filter>
                    <sld:MaxScaleDenominator>7500000</sld:MaxScaleDenominator>
                    <sld:PolygonSymbolizer>
                        <sld:Fill>
                            <sld:CssParameter name="fill">#777777</sld:CssParameter>
                        </sld:Fill>
                    </sld:PolygonSymbolizer>
                </sld:Rule>
            </sld:FeatureTypeStyle>

        </sld:UserStyle>
    </sld:NamedLayer>
</sld:StyledLayerDescriptor>

I've tried capital GEOMETRY. Still get the same error.


texts: Raster clip in ArcMap changing the pixel size

texts: I have a masked raster layer (ESRI grid format) of a large area that I'm trying to clip in ArcMap 10.3. Cell size is key for further analysis, thus I need to maintain the cell size of the clipped layer to 30 x 30 (the original masked layer is 30 x 30). However, when I clip the layer, the cell size for both X and Y reduces to 29.99806155. Unfortunately, this minute change does impact my analysis, and I don't know (a) why this is occurring, and (b) how to prevent it. I've tried adjusting the processing extent and cell size to prevent this change, but to no avail. Any advice?


texts: PostGIS: ST_Union fails with valid polygons?

texts:
I am working with this dataset. I have imported it into a Postgres 9.6 database with shp2pgsql as follows:

shp2pgsql -I -s 5070:5070 us_govt/PADUS1_4Fee.shp usgovt | psql -d mydb

I have made all the polygons valid:

UPDATE usgovt
 SET geom=ST_Multi(ST_CollectionExtract(ST_MakeValid(geom), 3))
 WHERE NOT ST_IsValid(geom);

Then simplified them, and re-validated them:

UPDATE usgovt SET geom_simpler=ST_Multi(ST_SimplifyPreserveTopology(geom,0.5));

To confirm, all the simplified polygons are valid (the following produces zero):

SELECT COUNT(*) FROM usgovt WHERE NOT ST_IsValid(geom_simpler);

However, trying an ST_Union on a group of the polygons fails:

# SELECT ST_Union(geom_simpler) AS geom from usgovt WHERE gap_sts='2';
ERROR:  GEOSUnaryUnion: TopologyException: found non-noded intersection between LINESTRING (-3.36702e+06 4.94515e+06, -3.36699e+06 4.94512e+06) and LINESTRING (-3.36699e+06 4.94512e+06, -3.36702e+06 4.94515e+06) at -3367018.7763304636 4945150.2479274161

Why could this be failing, and how can I debug it?

UPDATE: Tried buffering with SELECT ST_Union(ST_Buffer(geom_simpler, 0.000001)) AS geom from usgovt WHERE gap_sts='2';, but that fails with: ERROR:  GEOSUnaryUnion: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 3216270.6962009948 -12370.16660008162 at 3216270.6962009948 -12370.16660008162.

UPDATE 2: I've tried creating my own version of ST_Union with better exception handling, as suggested in this answer to a related question:

CREATE OR REPLACE FUNCTION safe_union(geom_a geometry, geom_b geometry)
RETURNS geometry AS
$$
BEGIN
    RETURN ST_Union(geom_a, geom_b);
    EXCEPTION
        WHEN OTHERS THEN
            BEGIN
                RETURN ST_Union(ST_Buffer(geom_a, 0.00001), ST_Buffer(geom_b, 0.00001));
                EXCEPTION
                    WHEN OTHERS THEN
                        RETURN ST_GeomFromText('POLYGON EMPTY');
    END;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;

It returns CREATE FUNCTION but when I actually try to use it:

SELECT safe_union(geom_simpler) AS geom from usgovt WHERE gap_sts='2';

I get HINT:  No function matches the given name and argument types. You might need to add explicit type casts..

Also, I'm not sure how safe this method will be, if it's unioning tens of thousands of polygons, because I don't understand the internals of how ST_Union works.

Is it possible that it could successfully union 5000 polygons, then get to the 5001st, hit an error, and silently replace the union of the previous 5000 with a null polygon, ultimately producing an extremely inaccurate result?


texts: Add Table View to Data View in ArcMap gives Assertion Error

texts:
I'm trying to create a Python AddIn for ArcMap 10.6 and the first process of the tool will be to add a tableview to a map document. At this point I just want python code to bring up rows from a table specified by a where clause when creating the table view so you can see this tableview in the data view display of the map document.

Here's my code so far:

import arcpy

# construct where clause
val = "\'%2%\'"
f1 = "DIST_NAME"
WC = """{} LIKE {}""".format(f1, val)

# create tableview
cList = arcpy.MakeTableView_management("CensusTract", "cList", WC)

# add tableview to dataframe
mxd = arcpy.mapping.MapDocument(r"C:\TableView.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
arcpy.mapping.AddTableView(df, cList)

And here is the error that I'm getting:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.6\arcpy\arcpy\mapping.py", line 114, in AddTableView
    assert isinstance(add_table, TableView)
AssertionError


texts: ArcGIS Mosaic to New Raster returning single value

texts:
I am using Mosaic to New Raster to combine two rasters that overlap. RASTERONE contains two different values, 0.4 and 1.1 with some areas of no data. RASTERTWO contains a single value, 0.03, in every cell where RASTERONE contains no data.
RASTERONE metadata:
Number of bands: 1
Cell size (X, Y): 50, 50
Format: Tiff
Pixel Type: floating point
Pixel Depth: 32 Bit

RASTERTWO metadata:
As for RASTERONE save
Format: GRID

I am combining them using Mosaic to New Raster with the following settings:
INPUTS RASTERS:
RASTERONE (entered into the dropdown first) then RASTERTWO
Raster Dataset Name with Extension: myraster.tif
Pixel Type: 32_BIT_FLOAT
Cellsize: 50
Number of Bands: 1
Mosaic Operator: First
Mosaic Colormap Mode: First

This returns a raster with a single value across the entire area. The raster's data range shows correctly in the legend: 0.03 -> 1.1. However, the actual raster is 0.03 across the entire area.
Switching the Mosaic Operator setting from "First" to "Last", without changing the order that rasters are entered into the dropdown, gives the same output: a raster of 0.03 across the entire area.
I have used this same process in the past with success.
What are the possible causes of this problem?


texts:
I'm new to QGIS.  I'm interested in geocoding addressed and have seen from research that the MMQGIS plugin is a great resource.  I do not see this as a preloaded plugin (most references online indicate that it generally is).

I downloaded a MMQGIS folder of the current plugin and tried to drop that under contents in the QGIS content folder but that was unsuccessful.  Any assistance or information (even just knowing why it isn't a preloaded plug in) would be great.


texts: Radiometric calibration of WorldView-2 imagery using ArcGIS Pro

texts:
I am attempting to work with worldview-2 imagery that I have downloaded and I am trying to radiometrically calibrate the raw imagery following the methodology outlined here. Since these are my first attempts, I wanted to make sure that I am doing it correctly.

In the first equation of page 9, one multiplies the image pixel value by the absolute radiometric calibration factor (Kband) and divides the result by the effective bandwidth (deltaK for each band). These values are contained in the .IMD file that comes with the downloaded imagery.

As an example, if a blue band pixel value is 4, and the .IMD data is

BEGIN_GROUP = BAND_B
    ...
    absCalFactor = 1.783568e-02;
    effectiveBandwidth = 5.430000e-02;
    TDILevel = 10;
END_GROUP = BAND_B

Then the resulting value would be (4 * 1.783568e-02) / 5.430000e-02 = 1.31

Is this the correct interpretation of the data and the correct method of calibration?


texts: Is there a way to get create label classes similar to symbology classes?

texts: Is there any way to get the same classes of label size as by the symbols? For example: i devided towns according to their population to 5 classes (as graduate symbolology), and want to label each class by different size of fonts. Thanks.


texts:
I use QGIS 2.18.26.
I have a layer with 10 million location points over which I lay a 1m by 1m raster.
Then I want to count the number of points in each 1m² of this raster.
When I launch the count I get an error without a log file.

This has worked for 1 million location points with a 5m by 5m raster.

How can i make this work for the 10 million points within each 1m²?


texts:
I am having some problems exporting from Excel to Arc.

Exporting directly from .xlsx to Arc

columns such as Name and latitude disappears all together

Exporting from .xlsx to CSV, then erasing formatting, then exporting to Arc

large numbers such as ones with 9 or more digits fail to display

Does anyone know a failsafe procedure to follow when exporting tables from excel to ArcMap 10.2 so that ArcMap retains all values?


texts: Zoom to a selected feature with Openlayers 3

texts:
I have a function that allows me to change color to a selected feature in my layer's source. Within the same function I would like also to change the map view by zooming in to the selected feature extent.

Since I have different layers (with different sources) I have grouped first my layers and my sources

layers_group = [];
layers_group.push(layer_1, layer_2);

source_group = [];
source_group.push(source_1,source_2);

My function is as follows:

(function($) {
//when clicking the select
document.getElementById("show_feature").addEventListener('click', function () {
    resetSelectElement(region_select);
//first step it loads the vector layer and sets the style to transparent
    var i;
    for (i = 0; i < layers_group.length; ++i){
        layers_group[i].setStyle(transparentStyle);
        layers_group[i].setVisible(true);
    };
//second step when the select changes (i.e. I make a choice)
    $("#show_feature").change(function() {
//it starts the function to change styles according to the selection made
        var selectedFeature = $("#show_feature").val();
        var index;
        for (index = 0; index < source_group.length; ++index){
            source_group[index].forEachFeature(function(feat){

    //if substring(selected text) exist within feature property('Region')
    //should return any value other than -1
                if (feat.get('Region').indexOf(selectedFeature)!=-1) {
                        //sets the style on the feature if it exists
                        feat.setStyle(selectStyle);
                        //should zoom to the selected feature polygon
                        var extent = feat.getGeometry().getExtent();
                        map.getView().fitExtent(extent, map.getSize());

                } else {
                //and if it doesn't exist switch back to a transparent style - not visible
                feat.setStyle(transparentStyle);
                }
            });
        };

    });
});
})( jQuery );

Everything works well except for the zooming, as I get a "TypeError: map.getView is not a function" and no zooming is performed.

Could somebody help me in solving this?


texts: Ortelius and Nicolosi projection in QGIS

texts:
I am trying to define ortelius oval projection and Nicolosi globular projection in QGIS using proj4 foramt as described in documentation but QGIS is unable to project into these projection even if they are supported in proj4 library... I need those projection for further analysis of old maps.

Does anybody know what could be the problem?


texts: Analysing Raster, Shapefile and DBF together using ArcGIS Pro

texts:
Note at the bottom for a possible solution, though it's not very elegant!
I’m at a bit of a brick wall – I think this probably has a really simple solution, but I’m not quite sure where to start. Using ArcGIS Pro.
Problem: I can’t figure out how to calculate the area of various categories in a raster file, based on country boundaries in a shapefile.
Description of raster:
Environmental data based on 15 categories of land cover. The attribute table has three fields: Objectid, Value and Count.
Value stores a category label, 1-15. I have a separate dbf which maps these to land cover (e.g. large forest, small forest, lake, sea) and Type of cover (e.g. forest, water), but I’m not sure if I should be bringing that in at this stage. I have one version where I’ve joined it to the raster att table, and another version where they’re still separate.
Description of shapefile:
This stores various boundaries. The attribute table has 6 fields: FID, Shape, Country, Province, Region, Area.
What I want to achieve:
There are five countries. For each country, I want to know the area of forest, the area of water etc.
Rough Illustration
Dbf table:

Raster attribute table:

Shapefile att table:

Goal:

Whether “Forest” and “Water” actually contain the km2 area, or just a cell count that I can convert based on cell size – I don’t mind. It would also be OK to still have Description (large forest, small forest) rather than grouped Type (forest), since I will be analysing everything in Excel anyway (sweet home, Excel).
I thought Zonal Statistics might be a good place to start, but I don’t think that will really achieve what I’m after, since I want to split out the data in the raster. I used that previously with the same shapefile, trying to get the population across each country. But I suppose that was a case of getting less detailed (summing across provinces), whereas now I want to get more detailed (split the “Count” field by country).
Can anyone point me in the right direction to get started?
Edit: I may have found a way to do this. First I extracted each country to a separate layer (select, right click, make layer from selected features). Then I used extract as mask on my main raster. I did this three times, once for each country. Now I have three raster files, each one corresponding to a unique country. I think this will do the job, but I can't help but feel that there's a tidier way.


texts:
I'm new to the environment QGIS, and my question and problem is when changing a layer coming from http://ide.cime.es/menorca/wms/base_referencia/

I want to load the 1-> 86 layer and at first everything works well, the problem comes when scaling to work at street level -- when I get to the scale of 1: 332,000 you can no longer see all the map.

Do you know why it may be that you no longer see anything?

Google translated from:

  Soy nuevo en el entorno de QGIS, y mi duda y problema es que una vez
  cargo una capa con el servidor WMS
  http://ide.cime.es/menorca/wms/base_referencia/

  quiero cargar la capa 1->86 al principio todo bien, el problema viene
  cuando cambio de escala para poder trabajar a nivel de calle y cuando
  llego a la escala 1:332.000 ya deja de verse todo el mapa.

  Sabéis a que puede deberse que ya no se vea nada.


texts: ArcGIS resampling rasters by maximum pixel value?

texts:
I need to resample a binary image of linear features like roads, rivers to a coarser resolution, conceptually combining multiple pixels into a single one. 1 indicates presence of the linear feature, and 0 otherwise.

To this end, I was wondering if there is an existing function for resampling a raster such that if any one of the smaller pixel is a one, the output resampled cell is one (, and zero otherwise)? In other words, I want the resampled cell to be the maximum of its constituent cells.

In ArcToolbox, the resample tool only have: NEAREST, BILINEAR, CUBIC, and MAJORITY. None fits the above need.

Is there a way to do this in ArcGIS (or other environments like PostGIS)?


texts: Calculating area with Europe Albers Equal Area in QGIS

texts: I am trying to calculate the area for a number of polygons in EPSG:102013, having changed it from EPSG:3857. When using the field calculator $area, the entire column was returned as nan. Has the transformation to 102013 not worked properly and if so, how can I force the data into 102013?


texts: Exporting a PDF from QGIS Print Composer for Avenza Pdf Maps

texts:
Alright I actually did this successfully a few months back following these steps http://www.avenza.com/resources/blog/2016/12/06/exporting-qgis-avenza-maps.

Now when I bring the map into Avenza, I get a "Map Not Referenced"

Now I have no clue what I am doing wrong.  I have tried and tried.  I have unistalled and reinstalled.  Tried to export as image then convert format.

A strange thing to note, I can not select reference map under Print Composer-> Export Settings -> Reference Map.


texts: What are the standards to implement scale bar?

texts:
I am working on .osm data files, and now I want to implement scale bar for available data.

Logically I have created the scale pixel is to geo location. I want to know whether there is any standard to display this scale bar.


texts: Exporting data from Google Earth to QGIS

texts: I am trying to export Google Earth searches as a kml file to QGIS. How come all existing locations of a search (e.g. Walmart) in the United States only shows select locations. I'm doing a project on food deserts and want to include all major chains and their locations in the U.S.


texts: Calculate field not working in Modelbuilder using ArcGIS Pro

texts:
I am trying to do a very simple calculation in ModelBuilder, where I use Add Field to generate a column with for volume data. I use Calculate field where I select the volume column as "field name", python as "expression type" and use the following code $feature.POLY_AREA * 0.5.
When I run it I get this error:

ERROR 000539:   File "", line 1
$feature.POLY_AREA * 0.5
^ SyntaxError: invalid syntax  Failed to execute (Calculate Field (5)).

I did this in ArcMap and it worked without any issue.


texts: Converting zero values in raster using average of neighboring cells in ArcMap?

texts:
I have an image

and I want to fill the empty spaces using the average of the neighbouring cells.  I setNULL these values and followed Fill in NoData gaps in raster using ArcGIS for Desktop?

There was no change despite the code running.

I tried mosaicking with another raster and replaced the empty bit with Zeros

.

How can then can I replace the zeros with the averages of the neighbours?


texts: Selecting isolated polygons which are not donuts but may overlap using ArcGIS Desktop?

texts:
I'm exploring a few different options for this process and I'm unsure what would be most effective. I'm trying to find all isolated, non-intercepting polygons in the same layer (not to merge). Most share boundaries, and a few overlap. Area varies so it unfortunately cannot be the deciding factor for selection (e.g. selecting by max area). I'm using ArcMap 9.3.1.

I created a 10 meter buffer around the polygons, erased the center of my buffer features, then select all polygons that intersect the buffer, thinking I'll just do a switch selection. The problem is that all polygons get selected, I guess because the selection is too tolerant. I don't want merely tangent/touching boundaries, I need them to overlap within the buffer (in effect, avoiding selection of those polygons each buffered feature is created from). The "are within" of course selected nothing. I wonder if dissolve can be used.

edit: as per request, here's a pic representing what I mean by isolated - no other polygons border it from any side.


texts: Data goes missing in Esri file geodatabase

texts:
We currently operate a system whereby survey data is collected in the field, exported in shapefile format, edited and then appended into a file geodatabase using the Load function in ArcCatalog. This all works nicely.

However, when I create new features within geodatabase feature classes using ArcGIS Desktop: Editor -> Start Editing -> Create New Features, or copy and paste features between feature classes during an ArcGIS Desktop editing session, they have disappeared the next time I access the geodatabase.

I have a feeling that this is because the newly created or pasted feature classes have no attribute data.  None of the fields are required apart from FID/OID. Can anyone confirm or say otherwise?


texts: Why does ArcGIS Field Calculator expression unexpectedly return a binary value (0 or 1)?

texts:
I am running simple field calculator expressions and experiencing some of the silliest results I have ever seen.

I am working on census data (living as a shape file) and trying to perform a simple [fieldName]/[fieldName] calculation on a newly created field.

ArcGIS happens to enjoy playing around with me... and returns values 0 and 1 in no particular order/reason that I can really justify.

Both input fields from the field calculator expression are of double data types. I have created a new field of double data type as well. However, the problem persists with other data types, including long/short integer.

Is this a bug? Has anyone else experienced this? Shortest route to create + update a new field other than moving between SPSS/Excel and ArcGIS all day long?


texts: Finding Ring Smaller than a certain tolerance

texts:
In Simplifying a shapefile for world timezones

http://efele.net/maps/tz/world/

all generalize/simplify algorithms fail on 5 countries/states, shown in this shape file

http://enclick.com/tz_world_handpicked.rar

The failure of the simplifying scripts is due to small rings below the size of the tolerance used for simplifying. The entire multi-polygon is deleted for each of these countries

How can I find and delete rings below a certain size with QGIS?


texts:
I'm new to GEE. I'm currently working on NDVI data analysis. So this is my code for MODIS Ndvi. I want to know how to divide all the images there inside ImageCollection by 10000 to get real NDVI values.

How do I get this work done?

{
var ndvi = ee.ImageCollection('MODIS/006/MOD13Q1');
var modisNDVI = ndvi.select('NDVI');
var collection = ee.ImageCollection(modisNDVI.filterDate('2017-07-01', '2018-12-31'));
var clipped = collection.mean().clip(fc)
}


texts:
I created a small dataset in google sheets. The columns are date, latitude and longitude, and size. I added the csv layer to my map, and then added the layer to timemanager plugin. for some reason, when i turn timemanager on, the layer no longer appears. any ideas why this might be happening?

I am using qgis 2.12.3 lyon.


texts: Is it possible to achieve sub-meter accuracies with a smartphone, without using an external GNSS-receiver?

texts:
I am currently working with QField on my tablet (Samsung Galaxy Tab S5e) and my smartphone (Samsung Galaxy A40). Since the positioning accuracy on my devices is rather low (absolute best case ~3m on an open field with clear sky), I am looking for lowcost ways to improve it in real time. Due to the high costs, I don't want to use an external GNSS receiver.
I've read that with RTK networks accuracies can be improved to sub-meter levels if you have raw measurements from your GNSS receiver, use a Ntrip-client and have access to correction data.
This discussion: GNSS raw measurements on Android: Does real time/post processing Differential techniques possible?, and this article and the corresponding video: https://www.gpsworld.com/how-to-achieve-1-meter-accuracy-in-android/; https://www.youtube.com/watch?v=vywGgSrGODU
gave me an overview of the topic, but I am still missing some things:
In the mentioned article/video from 2018 it says that there might be apps in a few years which will offer high-accuracy positioning with your smartphone. Does anyone know of an app that has been developed in the meantime and that achieves sub-meter accuracies with the GNSS-data received from a smartphone?


texts: Converting results of reduceRegion in GEE to something that can be passed as visualization parameter

texts:
I'm trying to dynamically set my visualization parameters based on the range of pixel values in my layer of interest. This seems like something that should be easy to do, since visualizing any given dataset requires you to understand the range of values that you're working with, and yet until recently, the only way I could figure out how to do this was by printing out the 1st and 99th percentile pixel values (obtained using reduceRegion), and then setting that by hand as my min and max in my visualization  parameters.
Now, I'm using .getInfo() to convert the results of my reduceRegion percentile function to numbers that can be passed to my viz params. However, GEE warns against using .getInfo(), and indeed it is working very slowly and crashing my webpage fairly frequently. Does anyone have a superior solution to dynamically set viz params based on the range of pixel values in a band that doesn't involve calling .getInfo()?
Here's my code, which works, but slowly and sometimes crashes.
// Load an input region: CA
var ca_counties = ee.FeatureCollection('TIGER/2016/Counties').filter(ee.Filter.eq('STATEFP', '06'));
var ca = ca_counties.union(); //making into one state instead of many counties

//Loading AOD data
var start_date = '2019-03-01';
var end_date = '2019-07-01';
var aod = ee.ImageCollection("MODIS/006/MCD19A2_GRANULES");

function get_state_mean(image_collection, band, state_borders, start_date, end_date) {
  var aod_state = ee.ImageCollection(image_collection)
        .filterDate(ee.String(start_date), ee.String(end_date))
        .select(ee.String(band))
        .reduce(ee.Reducer.mean())
        .clip(ee.FeatureCollection(state_borders))
  return aod_state;
}

// Get statewide mean AOD by pixel
var aod_ca = get_state_mean(aod, 'Optical_Depth_047', ca, start_date, end_date);

// Get scale (in meters) information from Optical Depth band
var aod_scale = aod_ca.projection().nominalScale();
print('Optical_Depth_047_mean', aod_scale); // ee.Number

//Getting viz parameters--min and max
var percentClip = aod_ca.reduceRegion({
  reducer: ee.Reducer.percentile([1, 99]),
  geometry:ca,
  scale: aod_scale,
  maxPixels: 3e7
});
var keys = percentClip.keys();

var palette = ['#00abff','#95bd00','e5e700','#cc2900','#a400c3'];
var viz = {palette: palette,
  min: ee.Number(percentClip.get(keys.get(0))).getInfo(), // <<< THE QUESTIONABLE PART
  max: ee.Number(percentClip.get(keys.get(1))).getInfo()  // <<< THE QUESTIONABLE PART
}

Map.addLayer(aod_ca, viz, 'ca aod');


texts: Opening PDF and other attachments in QField

texts:
I want to open a PDF file in QField, but I don't know how to configure a project to be created in QGIS.
Can you tell me how to do this?


texts: QGIS interpolation gives same value across the grid, what am I doing wrong?

texts:
Following this question, I've tried to use interpolation in QGIS to get a sense of the topography of  a map based on about 3500 track points with elevation. The elevation ranges between 213 and 353 meters with a 312 meters average.

I then try to use IDW with the following parameters:

It takes a few moment to process and then I get a raster with absolutely no variation.  The max is the average value and I've got a min of NaN. Of course, the rendering is all of the same color since there is no variance.

By asking for an ASC file instead of a TIF, I can confirm that it looks like it's creating a grid with the same value (average elevation) all over..

What am I doing wrong?


texts: What product does "EDN ArcGIS Enterprise Advanced" key unlock?

texts: I have a key to unlock something called "EDN ArcGIS Enterprise Advanced"   I tried to unlock ArcEditor with the key - seemed like the obvious choice, but that didn't work.  The key starts with first three letters "ECP".  It seems like it unlocks a Server, but I thought the server license cost $30K and I don't think we paid for that license, but it is available to me in the EDN portal customers.esri.com


texts: How to hide more than one QGIS attribute field at a time?

texts:
Using QGIS 2.14.1, I use the Layer Properties > Fields > Edit Widget to hide the display of an attribute in the attribute table, as shown in this screenshot:

Is it possible to hide multiple attributes all at once without the time-consuming process of going through the Edit Widget one field at a time?


texts: Calculating a "flag" based on two images and add it as a new band in Google Earth Engine

texts:
I have two images, high1 and low1. Pixels in high1 represent high areas. Pixels in low1 represent low areas.
I would like to create a new image that has a band with a flag for low and high areas. I would like the flag for low areas to be values of 10 and the flag for high areas to be values of 100. This new band would consist of all the pixels of high1 and low 1 and each pixel in that band has a value of 100 or 10.
Does anyone know how to accomplish this in Google Earth Engine or the GEE package in r (rgee)? I think this needs to be done through an image function.


texts: I am new to QGIS, I have been using ArcGIS for most of my maps.  I have two layers that are not overlapping properly. One is a shp file from a GIS clearinghouse and another from a delimited text source that I created myself.  I am using the same CRS for the delimited text source as the shp file from the clearinghouse, and I checked "enable on the fly CRS transformation" to see if that would solve my problem as well.


texts:
I am using gdal2tiles to convert a vrt to tiles but getting the following message:

  ERROR 4: 'Depth (Max).vrt' not recognized as a supported file format.
  Usage: [.....]
  gdal2tiles.py: error: It is not possible to open the input file 'Depth (Max).vrt'

Here's my code located in a batch file:

@python "%OSGEO4W_ROOT%\bin\gdal2tiles.py" -s EPSG:2278 -z 12 -w none -k "Depth <Max>.vrt" depth_max

I am running this from the OSGeo4W Shell. I had this working at one point so I suspect it is just a syntax problem with the quotes?


texts:
I have two Interstate labels on my map which appear to be different sizes even though they are the same font size and have the same marker text background size. The labels have been converted into annotation.

I am using ArcMap 10.2. What could be the problem?


texts: I am attempting to edit a linear feature layer and need to drag vertices at node endpoints to connect across slight gaps.  In ArcMAp, it was possible to lock onto a vertex and drag it to a new position.  Does QGIS offer similar functionality? I am using Edit> Reshape Feature.  Is this correct?


texts: Mask and crop a raster to multiple polygons in R

texts:
Is it possible to mask and/or crop a raster to multiple polygons, with the output being multiple raster's (i.e., a raster stack)?
Here's a reproducible example:
library(tidyverse)
library(raster)
library(sf)
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)

#  Create a data frame containing the centroids of a polygon
(df <- data.frame(
  x = c(180000, 180500),
  y = c(331700, 332200)))

# create a point shapefile from the dataframe, then convert to polygon
polygon <- df %>%
  st_as_sf(coords = c("x", "y"), crs = crs(r)) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POINT") %>%
  st_buffer(dist=200)

# plot the raster and polygon
plot(r)
plot(polygon, add=TRUE)

# mask and crop the raster
        m <- r %>% mask(polygon) %>%
          crop(polygon)

The above code returns a single raster. How to return multiple rasters for each polygon?


texts: Okay I can't figure out for the life of me how to log a live GPS track in qfield. I can add points and draw polygons, both of which are useful, but how does one track? I am trying to use qfield to collect field data for natural resources inventories so I need to track stream courses, trails, stone walls, etc., while keeping property boundaries and other features (shapefiles) visible. Anyone know how to record a track?


texts: VBA to VB.NET (Visual Studio 10) Migration - ArcMap Focus

texts:
I will be migrating several VBA applications currently in ArcMap 9.3.1 to Visual Studio 10 in the near future, and wanted some feedback on the general migration steps listed below.  I wanted to make sure that I did not miss any critical steps and that the order of the steps is correct.

1.)  Run the ArcGIS Code Converter tool to update old Object names.

2.)  Create a Visual Basic project and add all the mxd specific modules, and classes to the new VB project.  Re-create any forms or command buttons that existed in VBA.

3.)  Run Code Advisor for Visual Basic 6.0.

4.)  Open the VB project in Visual Studio 10 and run the Visual Basic Upgrade Wizard.

5.)  Manually update any common dialogs, MsgBox, and Win32 API code sections.

6.)  Add ESRI references and libraries.

7.)  Register new Visual Studio dll.

Here are some website references that I have used so far.

Migrating VBA/VB6 ArcObjects Applications to .NET
http://resources.arcgis.com/gallery/video/arcgis-desktop/details?entryID=A721AC1F-1422-2418-341B-DFD761945F47

How to migrate from VB6 to VB .NET (wiring events)
http://resources.esri.com/help/9.3/arcgisengine/dotnet/ca0956fb-b3d8-48a0-8c85-3559ae02dc0f.htm

ArcGIS Application Migration: Using the Visual Studio 2008 Upgrade Wizard
http://resources.arcgis.com/gallery/video/arcgis-desktop/details?entryID=93EC2C14-1422-2418-880E-0CFF2043C7FD

Thank You


texts: Reprojection from WGS to Slovenian Gauss Kruger in QGIS

texts:
I have gpx file (track) in WGS84 (EPSG: 4326) coordinate system from Garmin GPSMAP 62 which is right projected in BaseCamp and MapSource. When I save track as *.GPX and import it to QGIS and use Settings/Project properties/Coordinate Reference System/ and set as EPSG:3787 I get a wrong projected track (cca 300 m western of right position).
Thanks!


texts: Run GDAL on android and perform operations

texts: Could one run gdal on android devices to perform some GIS operation, like viewshed, contours, etc?


texts: Viewing/displaying fields with high character lengths in ArcMap?

texts:
I am trying to import an Excel spreadsheet with a field that has more than 255 characters.

Do I have to convert the Excel table to a file geodatabase table first in order to view this field in ArcMap?

How can I view this field in ArcMap without using a hyperlink to the Excel table and without writing a script?


texts: Does ArcGIS have something similar to Attachment Widget Type of QGIS?

texts:
I have shapefile and I would like to be able to add the path to a chosen file to its attribute in ArcGIS. It is a pretty straightforward task in QGIS.  I simply changed the widget type of the field to Attachment as you can see in the images. It allows me to add the path in the edit mode as can be seen in the GIF.
I am looking for a similar solution in ArcGIS. The closest I can found is Enabling attachment. However, I do not want to add an attachment, I am only interested in the path to file.
Any idea how I can achieve this in ArcGIS? Arcpy script is also okay.


texts: Removing intersect of one polygon from another using R?


texts: Python Shapely - filling small gaps between multi polygons

texts:
I have a set of rectangular polygons, which currently form a directional path.
The library I am using is Python Shapely / GeoPandas.

The issue I am having is as the path changes direction, there are small gaps and overlaps between the polygons.

I could iterate through all my polygon objects and calculate the angle in between each intersecting point and generate a new set of polygons but I figure there is a much better way.

I am fairly close - using the shapely buffer function, as seen in the following code snippet. I am able to remove the gaps, but it ends up merging my polygons into 1x. But I require them to remain 2x seperate polygons.

polygons = MultiPolygon([Polygon([(0, 0.5), (0, 2.5), (1, 2.5),(1,0.5)]),Polygon([(0, 0), (1, 1), (2,0),(1,-1)])])
polygons = polygons.buffer(2).buffer(-2)


texts:
I have a tiff raster in my PostGIS 2.0 database. I want to get the value of the raster at a  specified lat/lon point.

I have:

postgresql-9.1
postgis-2.0-svn with raster support
a raster imported into PostGIS its SRID is 3035


texts:
I'm working with Sentinel 2 satellite and I have sampled on the map 4 blue points indicating an irrigated areas, and other 4 red points for non-irrigated ones.

Then, I have considered 4 collections, one for each month, starting by May

#May
collection_mag = (ee.ImageCollection('COPERNICUS/S2')
              .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10))
              .filterDate('2019-05-01', '2016-05-31')
              .filterBounds(area))

#June
collection_giu = (ee.ImageCollection('COPERNICUS/S2')
              .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10))
              .filterDate('2019-06-01', '2019-06-30')
              .filterBounds(area))

#July
collection_lul = (ee.ImageCollection('COPERNICUS/S2')
              .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10))
              .filterDate('2019-07-01', '2019-07-31')
              .filterBounds(area))

#August
collection_ago = (ee.ImageCollection('COPERNICUS/S2')
              .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 10))
              .filterDate('2019-08-01', '2019-08-31')
              .filterBounds(area))

For each collection, I have calculated the NDVI, then reduced the collection into one image

def NDVI(image):
    return image.expression('float(b("B8") - b("B4")) / (b("B8") + b("B4"))')

img_mag = collection_mag.map(NDVI).mean()
img_giu = collection_giu.map(NDVI).mean()
img_lul = collection_lul.map(NDVI).mean()
img_ago = collection_ago.map(NDVI).mean()

This is the result by now (printing only the last month collection) - :

I want to plot a chart indicating the NDVI value, month by month, of the sampled points. Can anyone help me?


texts:
I would like to know if it is possible to link a text to a leaflet map.
I would like to highlight a country in my map when the mouse is hover a word in my text.
My map is created with leaflet library.


texts:
Somebody know how can I calculate the x,y,z,m parameters from a MultiLineStringZM in Postgis?

I would like to know x,y,z and m from the points what are made the polyline.

When I try calculate it by this way "*st_x,st_y, st_z, st_m*" postgis go up an error message:
ERROR:  Argument to M() must be a point

Any idea?

Thanks,

JM


texts: Understanding concept of Spatial Reference System?

texts:
I am trying to understand the concept behind the Spatial Reference System.

Why is it used and what is so specific about the SRID number?

I went through many definitions and almost all went over my head.

I thought, it was a reference to draw maps from a real world place, expressed in map, but then I inserted a point in PostGIS, with SRID 300 and the same point with SRID 350. I found the two points showing the same place, in QGIS viewer.


texts: Downloading image from Google Earth Engine

texts:
I would like to download an image from GEE to ArcGIS, with all information (coordinate, pixel value). I am trying but I couldn't.
How can I to do that?
I got this message:
Error: Image.clipToBoundsAndScale, argument 'input': Invalid type.
  Expected type: Image<unknown bands>. Actual type: ImageCollection.

This is the example of what I want to download (image)

This is my code:
var point = ee.Geometry.Point([-46.633286, -23.550510]);

// generate a buffer distance list
var distance_list = ee.List.sequence({
    start: 1000,
    end: 50000,
    step: 1000
});

var buffer_func = function (distance) {
    // new feature containing the buffer
    return ee.Feature(point.buffer(distance), {});
}

var buffers = ee.FeatureCollection(distance_list.map(buffer_func));

var clipToCol = function (image) {
    return image.clip(buffers)
}; /////////////////////////////

var AOD = ee.ImageCollection('MODIS/006/MCD19A2_GRANULES')
    .select('Optical_Depth_055')
    .filterDate('2019-05-13', '2019-05-28').map(clipToCol);
//.filterBounds(RMSP);

var viz = {
    min: 0,
    max: 350,
    palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};
Map.addLayer(AOD.mean(), viz, 'Optical Depth 055');

var outline = ee.Image().byte().paint({
    featureCollection: buffers, //////////////////////////
    color: 1,
    width: 1
});

Map.addLayer(outline, {
    palette: ['black']
}, 'buffers'); /////////////////////
Map.setCenter(-46.63203, -23.55221, 9);

Export.image.toDrive({
    image: AOD,
    description: 'image',
    scale: 1000,
    region: geometry
});


texts: Regenerate an isochrone through a widget in CARTOframes

texts:
Is there a way to re-generate an isochrone in a CARTOframes map from some kind of user interaction?
For example, modify the minute's value:
isochrones_gdf, _ = Isolines().isochrones(stadiums_df, [15*60], mode='car')

to something like
isochrones_gdf, _ = Isolines().isochrones(stadiums_df, [myDynamicValue], mode='car')

I need to create something similar to this using CartoFrames


texts: Expressing raster rotation with GDAL geotransform (python)

texts:
I want to generate a raster which is rotated with respect to its coordinate system. As I understand this requires modifying the geotransform array. I have not been able to find much information about this (though I did find a related SE question in Rotating rasters using GDAL Geotransform (and ArcGIS Desktop), and a thread on the GDAL-dev mailing list here, apparently started by the same person).

Most sites seem to only provide information on the structure of the unrotated geotransform. Based on the linked post & mailing list thread, I understand that the rotated geotransform should have the following elements:

[0] x-coordinate of upper-left raster corner
[1] cos(θ) * x_resolution
[2] -sin(θ) * x_resolution
[3] y-coordinate of upper-left raster corner
[4] sin(θ) * y_resolution
[5] cos(θ) * y_resolution

However, I am not sure if this is correct.

I have tried using python to save a rotated raster using the above as a guideline for my geotransform, but the output is not what I expected. I defined my geotransform as:

# Define raster array to save
array = np.arange(0,25).reshape(5,5)

# Define rotation & geotransform
rot = rot*math.pi/180
gt = (0,math.cos(rot)*x_res, -math.sin(rot)*x_res, 0, math.sin(rot)*y_res, math.cos(rot)*y_res)

# Save raster
save_raster(outfile, array, srs, gt)

And got the following result:

In a sense the raster has been 'rotated', but it appears my array has been resampled (which I do not want) while the raster cells are still aligned with the coordinate system x/y axes.

Am I just misunderstanding how the geotransform array works? Is there even a way to do what I am trying to do?


texts: How many levels for administration divide?

texts:
Excuse my ignorance, I just want to know how many levels for boundaries divide in a global level?

What's Admin 4 btw?


texts: Join Excel attributes to shape attribute

texts:
I am wondering if there is a way to join the attributes from Excel tables to shapefiles in QGIS. I have abbreviations in the shapefile and want to add the associated complete names from Excel file by automatic join. The abbreviation is also given in an Excel column.

Once there was a tool "join attribute" that is missing now. I tried to join by using properties/ join but then it is only possible to choose shapefiles.


texts: Implementing French language in PyQGIS Plugin

texts:
I build a plugin with PyQGIS.
If the user changes the language of QGIS interface then automatically my Plugin dialog and messages need to change that particular language (we just looking into French only).
Is it possible? If yes, how?


texts:
I am trying to add a widget to a webmap that will allow users to report wildlife spottings (birds, coyotes, etc). The client (City staff) would like to have to be able to review any reports prior to them being published online. We do want the users to be able to attach a photo to their report.

What is the easiest way to do this? Is it with the Report Feature widget? Is there a way to modify what menus come up with this (I don't need the severity menu)? Is there a way to allow photos to be added to the reports?


texts: Dissolving adjacent polygons and adding attributes

texts:
I am trying to merge together adjacent polygons and have the attributes added together for the merged features.

The layer is of certain census blocks across a state (thousands of features) and I would like to merge the adjacent census blocks and aggregate the population data of the merged features.

When I run the dissolve tool, I can get the adjacent polygons to merge, but the attributes are added for the whole table, rather than just the merged polygons.  Is there any way to have the feature attributes added together to make sense with the new geometry?

Another way I tried to do this was to use the dissolved layer and spatially join it to the census block population data.  For some reason, this is not accurate as the population totals are not equal.


texts: Filed calculator python code for auto increment not working

texts:
I have a bit of python code for the field calculator for numbering records from a given number by a given interval. It works most of the time, but sometimes it'll number out of order. It will use the right interval and end at the expected number, but not in the order I want. I've checked if the FID and ObjectID are all in numerical order, which they are. The only thing I can think of that may be causing an issue is that the field is a Text field. I've posted the little code snip below.

rec=0
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = pStart
 else:
  rec = rec + pInterval
 return rec


texts: How to measure similarity of approximately congruent polylines

texts:
My initial question was posted here. The answer of whuber provides a nice solution for the case when polylines are congruent.

However, the problem is that my polylines are approximately congruent (an example of the small difference of polylines is marked in black):

The objective is to measure the similarity between two polylines, assuming that small deviations (around 5%) are allowed and must be ignored. So, the similarity between above provided polylines must be close to 100%.

The question is how to measure similarity between such polylines using R? My approach is explained here, however it is very computationally expensive, and the solution of whuber assumes that polylines are congruent, not approximately congruent.


texts:
I do not understand how the operator works some times.

Here is an example that works,

"SELECT ST_Distance(ST_GeomFromText('POINT(-3.16496271127842 55.9262620593642)'),ST_GeomFromWKB(a.geometry)),ST_AsText(ST_Centroid(a.geometry)) FROM <table> As a;"

and here ia an example that fails,I use a simple execute method like this,

cc = conn.cursor()
a = ('POINT(-3.16496271127842 55.9262620593642)')
b = "SELECT ST_Distance(ST_GeomFromText(%s),ST_GeomFromWKB(a.geometry)),ST_AsText(ST_Centroid(a.geometry)) FROM <table> As a;"
cc.execute(b,a)

sometimes, the operator works fine and sometimes it does not. I have read documentation about the operators and how to use python with SQL but I have not found WHY it acts like this in a random kind of way.

Does anyone have any hints and tips?


texts:
What is the recognised naming style for field names in Geographic Information Systems (i.e. in shapefile layers)? Most datasets that I use with QGIS seem to use all caps (e.g. 'FIELDNAME') and underscores for word breaks ('FIELD_NAME'), but these have mostly been created with ArcGIS, so I was wondering whether this was because of a limitation to upper case in that package. What about other options, such as 'Fieldname', 'FieldName', 'Field-name'? Are any of these okay, or frowned upon?

What about conventions for filenames?


texts: Improve my arcpy script for simple querying

texts:
I am quite new to coding python in arcmap so I would like to get my code reviewed and hopefully improved.

I have two Shapefiles, one contain polygons of cities borders and the other is points representing city center of gravity. There is only one point per city but a polygon can coverage several cities.

My goal is that each point will have the objectid of the polygon containing it.

My code complexity is O(N^2) so its little disturbting and probably can be much better but I still dont understand good enough how to go over feature class and selecting each row at a step without looping and quering.

10X advance.

import arcpy

# fcYesuvimPolygon holds the location of the polygon feature class.
fcYesuvimPolygon = "h:/acadwork/Ortho/yesuvim_rec_final.shp"
# fcYesuvimPoint holds the location of the point feature class.
fcYesuvimPoint = "h:/acadwork/Ortho/yesuvim_name_point.shp"

# create layer for SHP.
arcpy.MakeFeatureLayer_management(fcYesuvimPoint,"fcYesuvimPoint_lyr")
arcpy.MakeFeatureLayer_management(fcYesuvimPolygon,"fcYesuvimPolygon_lyr")

# Go over the point SHP.
with arcpy.da.UpdateCursor(fcYesuvimPoint,["OID@","PolygonID"]) as cursorPoint:
    for rowPoint in cursorPoint:
        # query for selecting each row.
        qry = "FID = " + str(rowPoint[0])
        # select by attribute of each row FID
        arcpy.SelectLayerByAttribute_management("fcYesuvimPoint_lyr", "NEW_SELECTION", qry)
        # select polygons intersecting each point
        arcpy.SelectLayerByLocation_management("fcYesuvimPolygon_lyr", "INTERSECT", "fcYesuvimPoint_lyr", "", "NEW_SELECTION")
        # Go over selected polygon (one for each point)
        with arcpy.da.SearchCursor("fcYesuvimPolygon_lyr","OBJECTID") as cursorPolygon:
            for rowPolygon in cursorPolygon:
                # insert the objectid of polygon to point SHP field names Polygon ID
                rowPoint[1] = rowPolygon[0]
                cursorPoint.updateRow(rowPoint)


texts: QGIS 3x problem with rescaling vector (shapefile) layer

texts:
I would like to rescale my shapefile layer in a particular location on my map in QGIS.
Unfortunately every time my layer is rescaled in completely different location.
According to this query
Resizing shapefile features proportionally in QGIS
I can do it with Affine transformation.
Unfortunately when I run the affine transformation and fill it with my XY coordinates as well as the scale factor I am getting my layer away from my working area.

In the X and Y I put the coordinates from Gridreferencefinder.com
In the scale factor I put accordingly 1.3 for X and Y and this is working. The worst thing is the layer itself, which appear at completely different coordinates (bounded black).
I found some solution here:
Computing parameters for QGIS Affine transformation?
but it looks like it's not the quickest one.
The problem is exactly the same when using SAGA vector transformation and GRASS v.transform.
Is there any way to fix it?
Or possibly can I scale my layer manually using i.e. advanced digitizing toolbar?


texts: How do I use env parameter variables in an SLD colormap?

texts:
We have data being represented by the colormap copied below. I need to use env variables to change the hard coded color steps into whatever's input to the url parameters. I know you can use<ogc:Function name="env"> to pull variable values out of the url but I'm not really sure how to use that to replace the hard coded values. The colors aren't changing, just the quantity values. The getMap url would have something like &env=step1:1.0;step2:1.5; etc...

  <RasterSymbolizer>
    <Opacity>1</Opacity>
    <ColorMap type="ramp" extended="true">
      <ColorMapEntry color="#000000" quantity="0" opacity="0" />
      <ColorMapEntry color="#0000FF" quantity="0.75" opacity="1" />
      <ColorMapEntry color="#00FF00" quantity="1.25" opacity="1" />
      <ColorMapEntry color="#FFFF00" quantity="1.75" opacity="1" />
      <ColorMapEntry color="#FF6600" quantity="2.25" opacity="1" />
      <ColorMapEntry color="#FF0000" quantity="3.00" opacity="1" />
    </ColorMap>
  </RasterSymbolizer>

I tried doing what http://docs.geoserver.org/latest/en/user/styling/sld-reference/rastersymbolizer.html suggested by using

<ColorMapEntry color="#0000FF" quantity="${env('step1',1)}"/>

but that didn't work. Geoserver throws

  IllegalArgumentexception: Argument "quantityString" should not be null.

Which makes sense.

So the CQL expression didn't seem to work as I hoped. How do I get that to work or use the ogc function?

Using Geoserver 2.7.


texts: How to call a stored pgsql function from Geoserver?

texts:
I created a stored function in PostgreSQL : myfunction(a text, b integer).
I want to create an SQL View in Geoserver using this query : Select * from myfunction('%nom%', %likes%).

My problem is that Geoserver don't recognize my function. And thus when I click on "refresh" nothing is displayed.


texts: Hitting Memory Limits of PostGIS ST_ConcaveHull / Incomplete Shrink Wrap

texts:
I am trying to use PostGIS 2.0 on Postgres 8.4 to do some mapping. I am working with the NOAA Shapefiles, specifically their US States & Territories file located here:

http://www.nws.noaa.gov/geodata/catalog/national/html/us_state.htm

The app I'm building  needs to do polygon in polygon tests, and point in polygon tests (PiP). The shapefiles as-is are too detailed, and the queries are taking too long. I'm trying to get PostGIS to simplify the complex state shapes, I want to achieve a high quality shrink wrap of each state, so that rivers and bays are included in the states but small coastal islands or other tiny landmasses are ignored . My goal is to have two tables: both the states raw geometries table for high quality output situations, and also a separate states_simple table where I can use pre-compressed geometries to speed PiP tests.

I have finally figured out how to properly implement ST_Concavehull, but I'm running into limits of either the server hardware, or the configuration of postgres.

Here is an image to illustrate the issue I'm running into:

Processing "square" states like Oklahoma using ConcaveHull process quickly and without any issue. States that have SOME coastal islands / complexity like Georgia process a little slower, but reasonably fast and come out in a way that will be useful for PiP test (i.e. tight shrink wrap):

However, some states like New Jersey are exhibiting some of the same problems Florida has using ConcaveHull:

The top and bottom of New Jersey are EXACTLY what I want, but the large part of PA including Philladelphia, and the inclusion of Staten Island basically make this geometry useless for my purposes.

So right now the query I'm using to retrieve these polygons is as follows (Example is FL):

SELECT ST_AsGeoJSON(ST_ConcaveHull(ST_Union(geom),.4)) as geom  FROM states WHERE id = 11;

tl;dr:

If I turn down the tolerance from .4 to something like .6, the state comes out like a giant triangle. .4 seems to be the floor limit of ConcaveHull, anything lower than that and queries run for 3+ minutes without returning results.

I've been following this article for tips and steps to increase performance of Postgres and PostGIS:

http://postgis.net/workshops/postgis-intro/tuning.html

Changes I've made to postgresql.conf:

shared_buffers = 256MB
work_mem = 16MB
wal_buffers = 1MB
checkpoint_segments = 6
random_page_cost = 2.0

The server all this is running on is a Turnkey Linux LAPP Server running on an Amazon EC2 Medium M1 Instance (64 bit, 1cpu, 3.75gb RAM)

Am I going about this all wrong? Is there another freely available set of "smoother" CONUS states I can use? I'm theorizing that the relative size (read: huge) of the state geometries is what's gumming up PostGIS/Postgres... just not really sure how to diagnose or monitor that beyond running

ps aux | grep ^ postgres

and watching how long queries take to run.

My main questions:

How can I monitor the resource usage of this process better?
Am I using the right PostGIS functions to accomplish this goal?
Is it a matter of turning down the tolerance of ConcaveHull, and increasing system limits?


texts: Reducing Land Cover Data Sets in Earth Engine

texts:
I am trying to use earth engine to reduce the USGS land cover dataset and calculate the frequency of specific land cover types in several regions. The output that I am getting creates a single dictionary with a list of land cover types and values for each county (region).

I am curious if I can break this dictionary up to allow a table to be exported to google drive with each land cover type having a unique column?

Here is my code:

// Load a FeatureCollection of Maine Counties.
var Counties = ee.FeatureCollection('ft:1S4EB6319wWW2sWQDPhDvmSBIVrD3iEmCLYB7nMM').filter(ee.Filter.eq('StateName', 'Maine'));

//Imports NLCD Land Cover Data
var LandCover2011 = ee.Image('USGS/NLCD/NLCD2011')

// Clip the image to the polygon geometry
var LandCover2011 = LandCover2011.clip(Counties);

// Extract the landcover band
var landcover = LandCover2011.select('landcover');

// Print out the frequency of landcover occurrence for each county
var frequency = landcover.reduceRegions({
  collection: Counties,
  reducer:ee.Reducer.frequencyHistogram(),
  scale:30
});

//Prints Feature collection
print(frequency);

//Exports table to drive
//Export.table.toDrive(frequency, "MaineLandCoverData");


texts: Are there any open source Python modules compatible with arcpy that will allow complex SQL queries?

texts: I am trying to handle string (non-geometry) duplicates in a large attribute table of a shapefile without having to convert to gdb and back or connect to a postgres db. I need to be able to use a count function for my query. I'd prefer not creating an output file, but rather simply updating the file itself. I've looked at the cursors but I cannot seem to accomplish what I need to with them. I essentially need to get a count for duplicates, do a comparison, and then accordingly append values of one field with the values of another field.


texts: Adding attribute table to raster in PostGIS?

texts:
I would like to add attributes to a raster stored in a postgis database. I use the following command to insert the raster in the database:

raster2pgsql -s 4236 -I -C -M /home/user/mysrc/outlabeled.tif public.test | psql -d template_postgis

The raster outlabeled.tif is in fact a one-band raster with a label per group of pixels (object). Then I have computed .txt files with attributes that I'd like to add to the raster (like an attribute table in a shapefile). Is it possible to do it (and how)? So far I've the impression that the only way to do it is actually to add a raster band (and not text files). If it's not possible, would an attribute table added with gdal before the importation do the job?


texts: Applying base heights to raster is giving me different results

texts: I have applied the DEM's base heights to both the raster image, the raster DEM, and a vector file.  The DEM and vector file match up but the satellite image does not match up at all.  Any suggestions to what could be causing this?


texts: Adding multiple users to an enterprise database in one go

texts:
I couldn’t figure out how to add multiple users to an Enterprise database in one go.

Please consider the following situation:

1 ) I have an  enterprise database named “ Test02”

2) Then I try to add multiple users to it in one go, these users are stored before in the SQL server.

3) Now the problem that I can’t add them all to the database and I need to add each user alone and set a password to it.

Is there any method to see all users in the SQL server and add (all or some ) of them to the database.


texts: Creating table containing all filenames (and possibly metadata) in File Geodatabase?

texts:
I need to create a metadata-table for all files in a File Geodatabase.

Is it possible to automatically be updated with new filenames, and possibly even metadata (as created in the ArcGIS metadata editor)?


texts: Dissolving polygons based on multiple attributes with Python (shapely and fiona)?

texts:
I would like to dissolve a polygon shapefile based on two different fields, using Python and open source libraries.

I'm just able to dissolve the shapefile using one single field, with Gene's answer to Dissolving polygons based on attributes with Python (shapely, fiona)?

Here's the code he suggests, that should be the starting point for what I need:

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
with fiona.open('cb_2013_us_county_20m.shp') as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open('dissolve.shp', 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'STATEFP' field
        e = sorted(input, key=lambda k: k['properties']['STATEFP'])
        # group by the 'STATEFP' field
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['STATEFP']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})


texts: Need to supply a projected polygon shapefile in latitude and longitude?

texts: I am supplying a polygon with multiple features to be uploaded as a geofence layer in a Web AVM. The dataset is projected (MGA94). For some reason they can't work with metres and have requested the dataset in latitude and longitude. How can I convert or "unproject" my dataset? I have tried setting the coordinate system of the dataframe in ArcMap to a geographic coordinate system (GDA94) and exporting the data that way but the value of the coordinates doesn't change.


texts:
I have downloaded an osm file from Geofabrik download server and opened this .osm.pbf-file in QGIS. So there are 5 classes - lines, multiline strings, multipolygons, other_relations, and points.

I wanted to add some additional number of points and give them attributes in this vector layer. And finally, save it as .osm.pbf again, As I need the .osm.pbf for osmosis.

Is there a method to do the same? As editing in that layer is not happening.


texts: Access SQL database in postgres for routing

texts:
i am trying to access sql database in pgadminIII. Using postgres_fdw I am trying to connect. But I am struck in the middle..Following are the steps I followed

CREATE EXTENSION postgres_fdw;

CREATE SERVER ws253 FOREIGN DATA WRAPPER postgres_fdw OPTIONS (host '192.168.0.253', dbname 'allahabad', port '1433');

CREATE USER MAPPING FOR postgres SERVER ws253 OPTIONS (user 'sa', password 'track');

CREATE FOREIGN TABLE routing
(
 gid serial NOT NULL,
 name character varying(255),
 district character varying(255),
 state character varying(255),
 the_geom geometry,
 name_hindi character varying,
 name_hindi_display character varying(255)
)
SERVER ws253;

What is the next step I have to do?


texts: Saving fields after a raster join with python arcgis

texts:
I am working on a raster processing script in which i join a categorical raster with a table that includes some extra data on each category.  After much hassle, I have the join working, but I can't seem to save the raster with the fields included.  I would like to be able to use the lookup tool by field name later on my script.

Here is my code.

AddJoin_management("joinedRaster", "VALUE", "joinTable", "joinField")
F = Raster("joinedRaster")
F.save("Ras_j")

If I try to do:

fieldList = ListFields(F)
for field in fieldList:
     print("{0} is a type of {1} with a length of {2}"
        .format(field.name, field.type, field.length))

with F, then I see that join has been completed successfully.  But if I try to do:

a = Raster("Ras_j")

Then I only see the fields from the original raster dataset, without the joined fuels.  Any advice on how to make sure that the join persists?


texts: Matching sde database port to FME Reader?

texts:
I had been trying to connect a MSSQL Server/SDE 9.3 to an FME (2013) Reader for 2 days until I noticed from an earlier fme project that each sde Database has a unique port number. I thought the sde connection was as standard 5151 but it seems that my database was 5158.

Does every database have a seperate port number? Or is this somehow defined per sde connection?

How can I get a list of which database is connected to which port?


texts:
I joined two tables to a shapefile and the output table has a field from each. Both fields represent the same thing so I'd like to combine them into one. You can see a clip of my fields below, where most of nulls are populated in the adjacent field.

Simply enough, I just want the final output to read as (based on the fields above)

<null>
Complete
Verified
Verified
Verified
Verified
Construction
Construction
Construction
<null>
Verified
Approved
Complete
Verified
<null>
Complete
etc

I was thinking about using field calculator to say "If field 1 is None, field 1 = field 2" but that either doesn't work or I was getting the syntax wrong. Also, I have a few sets of three fields that I need to merge so that wouldn't work for that unless I calculated twice?


texts: GeoServer cascading WMS, FeatureInfo is lost

texts:
I'm trying to cascade a WMS service in GeoServer (http://grb.agiv.be/geodiensten/raadpleegdiensten/GRB/wms?REQUEST=GetCapabilities&SERVICE=WMS).

I can view the cascaded layers just fine, but whenever I try to get FeatureInfo from a cascaded layer I get an error message saying that the layer isn't queryable.

I can get the FeatureInfo from the original service, so I know that there is actual FeatureInfo on the original layers.

I've combed through every available WMS setting in GeoServer but only found the 'Queryable' option in the 'WMS Settings' panel in the Publishing tab of the cascaded layers, which is obviously enabled.

Is there anything I'm missing, or is it just not possible to get FeatureInfo from a cascaded WMS service?


texts: QGIS 3.4 - create mbtiles from QGIS server

texts:
We are trying to develop a workflow to generate .mbtiles (for offline tiled map viewing) from a QGIS project.

A few years back, in QGIS 2.4, we went with QGIS -> multiple tiffs -> multiple GeoTIFFs -> VRT -> tiles.  That introduced some image quality loss.  See that question here.

Even back then, it seemed the real way to get around image quality loss would be to use QGIS server and download the tiles that it serves up.  It was over my head and I did not make any attempts along those lines.

Now, QGIS 3.4 serves WMTS so the serve-and-download method seems to be more feasible.  I set up WAMPServer (Windows 10) and QGIS Server (3.4), and I am getting pretty close but am not having luck finding a good way to grab all the tiles served up by QGIS server.

I can grab the entire map image of a simple locally-served map using WMS:

http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?MAP=D:\Documents\GIS\Projects\serverTest.qgz&SERVICE=WMS&REQUEST=GetMap&WIDTH=400&HEIGHT=200&CRS=EPSG:4326&LAYERS=Community_Boundaries

I can grab one tile using WMTS:

http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?MAP=D:\Documents\GIS\Projects\serverTest.qgz&SERVICE=WMTS&REQUEST=GetTile&Version=1.0.0&Style=default&TileMatrixSet=EPSG:4326&TileMatrix=1&Format=image/png&Layer=CommBnd&TileRow=1&TileCol=1

(At the moment, that is serving up a solid black 256x256 image even though the full WMS map looks correct with a few town boundaries, but, I will figure that out separately.)

What's the best way to grab a full tile set or mbtiles file based on the WMTS single-tile request above?


texts:
I would like to get the average value of specific attributes. I have a table with columns like:

ID, P1, E1, Age1, P2, E2, Age2, P3, E3, AverageAgeP, AverageAgeE

Sometimes there is a value for P, sometimes there is value for E - otherwise it's NULL. I would like to calculate the average value for all cells of P using the value of the Age-field defined in this group. The same I would like to do with E.

The difficulties I see, is to select all P which are not NULL, summarize the age belonging to it and divide by the number of values used - of course excluding the NULL-values.

As I have more then 3000 records I need to do it automatically using ArcGIS 10.3 . Maybe using Python and loops...


texts: Save Mapbox-Vector-Tiles (PBF) on indexeddb and use them offline

texts:
I'm developing a PWA(Progressive web app) where I need to retrieve PBF files from a POSTGIS database, then save them on an offline layer, to work as a offline-first app, right now I'm using a PGRestApi

PGRestAPI:https://github.com/spatialdev/PGRestAPI

I'm trying to save this tiles using leaflet offline layer:
Leaflet offline layer
Is there a more efficient way to save those pbf vector tiles on indexeddb and then use them on a leaflet map?


texts:
I have a point cloud (.ply format) that I got after using Structure from motion techniques (with MicMac). I would like to calculate its volume with r.report in GRASS.
Will r.report work on a ply file?
If not, how can I convert my file into something that will work with r.report?


texts: Can I ignore very low levels of spatial autocorrelation in my data?

texts:
I have >1 million postcodes clustered in ~7000 middle super output areas for England and Wales which I'm analysing using MLwiN. I have used a multilevel model with random effects and have clustered the postcode observations at the middle super output area.

My concern is that when analysing the middle super output areas using a Queens 1st order contiguity weights matrix (in GeoDa), Moran's I came to 0.00324916 (p=0.013) (see graph below). While there is some spatial autocorrelation, this seems to be a very, very small amount - is it small enough for me to ignore?

If I couldn't ignore this then I would have to cluster at the regional level leaving me with only 32 clusters. As I'm using MCMC estimation this would dramatically increase the computational intensity of my model runs (leading to days per run instead of hours), hence I would quite like to avoid this if possible.

.


texts: Why won't this PostGIS query display in QGIS

texts:
I am having trouble rendering the following PostGIS query in qgis1.803.

SELECT
s.musym,
st_intersection(s.the_geom,f.the_geom) as geo,
nextval('"soilsplittersequence"'::regclass) as k
From public."SoilsSpatial" as s JOIN public."Fields" as f
ON ST_Intersects(s.the_geom, f.the_geom)
where f."FarmNumber"=37

When I run this query in the db manager module it returns 500 ish results, but I can't display it.   I also can't display this query if I use it to create a view in my postgis database.

The only way I can get it to display is if I make the query into a create table.

I'm assuming that the problem is related to my use of a sequence to generate a unique id.  The reason that I'm using this technique is that there is a many to many relationship between the two tables.

Thanks for the advice.


texts:
I am trying to learn how to develop QGIS plugins but I struggle to get going by some simple principles.

As a test I am just trying to connect a function to a button. So when pressed the function run. Where to put the function and connections are not obvious to me.

This example give me an error message TypeError: test_print() takes no arguments (2 given)
I cannot see how this give test_print 2 arguments?

def run(self):

    # Run test_print method when pushButton is clicked
    self.dlg.pushButton.clicked.connect(self.test_print)

    # show the dialog
    self.dlg.show()

def test_print():
    print 'It works!'

Does anyone has any general advice how to wire up buttons and functions properly?

EDIT: Working test example

This worked.

def run(self):

        def test_print():
            QMessageBox.about(None,'test','It works!')

        # Run test_print method when pushButton is clicked
        self.dlg.pushButton.clicked.connect(test_print)

Now, how can i run the test_print method if its not inside def run(self)?


texts: Cannot change layer name in GeoServer 2.18

texts:
I install the GeoServer 2.18 in windows server 2012 r2, and with tomcat 9, JRE 1.8.
it works fine.
but I can not change the layer title/name.
when I want to change the name, and I will get the error message

java.lang.reflect.UndeclaredThrowableException

I find some way to try to fix, like use the Xmax=512M,
Could anyone help me to solve this question?
Thanks


texts: How to query points that fall outside their designated regions?

texts:
I have a shapefile with 250 points. Each point has a region attribute. In another shapefile I have all the possible regions with their names (it is polyline geometry).

Unfortunately the coordinates of the points are sometimes wrong.

I would like to query which points don't lie in the region specified by their region attribute.


texts: CREATE TABLE AS with SELECT clause throws error in QGIS DB Manager

texts:
I am not able to run the CREATE TABLE AS command in QGIS DB Manager (Win 10, QGIS 3.4.5-Madeira. I am getting the same error in QGIS 2.18 (Win 8).

Query execution error on CREATE TEMP VIEW _tview AS CREATE TABLE VV as SELECT 55 AS FF: 1 - near "CREATE": syntax error

Here is the error message in a screenshot:

Actually, I see even the following does not work on my Win 10.

   CREATE TABLE CUSTOMERS(
   ID   INT              NOT NULL,
   AGE  INT              NOT NULL,
   SALARY   DECIMAL (18, 2),
   PRIMARY KEY (ID)
);

Am I missing something in syntax, or the DB manager is missing something?


texts: Generate new coordinates with distance + azimuth

texts:
I currently have the following Excel data for a list of animals:

Lambert 93 coordinates of the observer
Azimuth of the animal
Distance between the observer and the animal

I would like to get the Lambert 93 coordinates of the animals.
Is there a way to do this with ArcGis 10.3?


texts:
I have the following scenario:

1) I need to create GPX files for use in testing of a mobile device.

2) The tracks to be described within will span large distances over long time periods (i.e. it's not practical to just go out and record a track).

3) The people performing the testing are not very GIS knowledgeable and will struggle with things like manual file conversion, etc.

I have in mind some simple application (it could easily be a QGIS python plugin) in which you set certain parameters (time, direction & speed of track "stages") and then generate a GPX file accordingly.

But before I throw myself into something like that, does anyone know of a simpler way to do this?

I've looked into various options (creating a tour in kml and then converting into gpx, manually cut-and-paste the xml together by hand, etc.) but it all seems like a lot of hard work for what I would have thought was quite a common use case...


texts: Select multiple features in OpenLayers 3

texts:
I would like to use a click interaction to select multiple features (I have two vector layers) and show their names inside a infobox. What is wrong?

var selectClick = new ol.interaction.Select({
 condition: ol.events.condition.Click,
 toggleCondition: ol.events.condition.shiftKeyOnly,
});
var select = selectClick;
map.addInteraction(select);
var selectedFeaturesClick = select.getFeatures();
select.on('add', function(e) {
    var info = [];
    selectedFeaturesClick.push(feature);
    info.push(feature.get('DIVISION'));

    if (info.length> 0) {
    infoBox2.innerHTML = info.join(', ');
    }
});

I also tried with this following lines but still not working well. I would like to use a click interaction to select multiple features and show their division names inside a infobox2. Using the following code, if I add a new feature (using the Shift key) the infobox2 is updated with the additional feature name (DIVISION) but it is repeated again in the infobox2 when I try to remove this feature from the selection (with the alt key). I would like to deselect a feature and have its division name removed from the infobox2. How is possible?

var infoBox2 = document.getElementById('info2');
var selectClick = new ol.interaction.Select({
 addCondition: ol.events.condition.shiftKeyOnly ,
 toggleCondition: ol.events.condition.never,
 removeCondition: ol.events.condition.altKeyOnly,
});
 select = selectClick;
 map.addInteraction(select);
var features = select.getFeatures();
var info = [];

var displayFeatureInfoClick = function(pixel) {
    var features = [];
    map.forEachFeatureAtPixel(pixel, function(feature, layer) {
    features.push(feature);
    });

        if ( features.length > 0) {
            for (var i = 0, ii = features.length; i < ii; ++i) {
            info.push(features[i].get('DIVISION'));
            }
      infoBox2.innerHTML = info.join(', ');
    } else {
      infoBox2.innerHTML = '&nbsp;';
    }
  };

map.on('click', function(evt) {
displayFeatureInfoClick(evt.pixel);
});


texts: Why would Tools toolbar of ArcMap be partially frozen?

texts:
Whilst using ArcMap one of my students 'froze' bits of the Tools Toolbar e.g. Zoom In/Out and Go To XY.

I couldn't figure out why.

In the end we reloaded the last saved file and it was fine.

Any ideas on what we did?


texts: Updating Web Map Service (WMS) by default when data amended in Oracle source tables?

texts:
I have a WMS which source data is held in Oracle.

Every year an amendment to the view in Oracle is run to update records.

Does the web map service need to be re-published once these amendments are made or are these updates to the data automatically reflected in the WMS?


texts: Setting a Background Color with RVBT field with python

texts:
I'm trying to set the color of a background of a layer. To do this i created an attribute "Color" with the RVBT numbers associate to each feature. As in the picture, i want to set the color with those features to set differents color, but i don't know how to do this with python.

Here is an image :

Here is my code (I set red color for example) :

background_color = QgsTextBackgroundSettings()
background_color.setFillColor(QColor('red'))
background_color.setEnabled(True)
layer_settings  = QgsPalLayerSettings()
text_format = QgsTextFormat()
text_format.setFont(QFont("MS Shell Dlg 2"))
text_format.setSize(12.5)
text_format.setColor(QColor("black"))
text_format.setFont(QFont("MS Shell Dlg 2",11,QFont.Bold))
text_format.setBackground(background_color)
layer_settings.setFormat(text_format)
layer_settings.fieldName = "Valeur"
layer_settings.enabled = True
layer_settings = QgsVectorLayerSimpleLabeling(layer_settings)
layer.setLabelsEnabled(True)
layer.setLabeling(layer_settings)
layer.triggerRepaint()
qgis.utils.iface.layerTreeView().refreshLayerSymbology(layer.id())


texts: Features color based on layer's color in QGIS

texts:
Is there a way to synchronize the selection color with a layer's color in QGIS?.
Following this post https://gis.stackexchange.com/questions/111784/displaying-only-selected-features-on-map-in-qgis#=  I am able to change the order and apparently color of selected features with a based-rule symbology, however the selection color does not take rule based color only width and line type.


texts: Using Append Parameter in the Add Locations Network Analyst Tool in ModelBuilder?

texts:
I am trying to get ModelBuilder to create a route for each row in a .csv file. Each row has one start and stop locations (Point A and Point B). I need to add both start and stop locations to a network layer to create a route. I then need to delete these stops, and iterate back through, creating a route for each record between each start and stop point. In the Add Locations tool help, I read that the Append Parameter might be what I need to use. In the tool dialog box, there is just a checkbox, but in the Python reference it says

•   APPEND — Adds the new network analysis objects to the existing set of objects in the selected sublayer.

•   CLEAR — Deletes the existing network analysis objects and replaces them with the new objects.

I think I need to do both Append and Clear actions to get what I am trying to achieve but it is just a checkbox in ModelBuilder.

Will I need to use Python code to be able to append, then clear, stops in an iterator model?

I am not very strong writing Python Code and was hoping to use ModelBuilder to get what I need.

I added my model builder file for clarification. The problem lies in the Add Locations tool. It only accepts one pair of Lat/Lon coordinates


texts: Why does not equal computing line in `OpenLayers` with `ArcGIS` when i use measure control?

texts:
My map projection is: "EPSG:900913" ,Display projection is "EPSG:4326" and units is "m", however does not equal the length is calculated in OpenLayers with ArcGIS, in fact Openlayers show the length more than the actual length when i use measure control.

For example:

the calculation length in ArcGIS is:26,072.443575

and in OpenLayers is: 3.0757.

How to resolve this problem?


texts:
I know this question has been asked quite a few times before but I couldn't find an answer that suited my situation at all.

So this is what happens. Currently, I am using ArcGIS to map the sewer connections, which have a dataset attached to them (size, length etc.) Once a week I export these to a personal geodatabase called "Sewage.mdb". This file is used by QGIS (which several others use at work) and shows the updated points with the associated detail. All the other people at work do is open the QGIS file shortcut on their desktop and it loads the existing qgis model data from the .mdb file.

Now I am trying to update everyone from QGIS 1.80 to 2.8.2, when I have trialed this I get an error that my sewer connections are going to cause an error because it can't read from the .mdb file. There I get the error, "Sewage.mdb is not a valid layer and can not be added to the map"

TL;DR
ArcMap 10 export > Sewage.mdb > QGIS loads Sewage.mdb > error Sewage.mdb is not a valid layer and can not be added to the map.

Also I have tried to add the layername after to see inside of the .mdb file with no luck
"|layername=DCCC_Assets_DBO_SewerNetwork_Junctions"


texts: GeomorphicFloodArea Plugin is not working

texts:
I'm trying to get the plugin GeomorphicFloodArea to work but cannot.

I tried using QGIS 3.10 but got an error message, so I downloaded version 2.18.28 and it still will not work. I get an error message saying:

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

Not sure what else to try at this point.


texts:
Does anyone have an idea of how you could allow the user to draw a polygon on a resium map?

Resium is a react library using cesium. I have an app that has entities placed on the map and I'd like to allow the user to select an area on the map by drawing a polygon using their mouse and be able to filter out entities not in the polygon the user created.


texts: Seeking Demographic Information of U.S. Census Blocks

texts: I am struggling to find the blockwise (not blockgroup or tract) data like employment status, household income, commuting distance, list of place of interest, industry information. By using American Factfinder I get easily the updated data of these factors by tract and 2010 version of blockgroup. ACS Employment Status Data by Block Group, 2006-2010 is given here.


texts: Is it possible to do a visibility analysis and create a 3d - cone as a result in ArcGIS?

texts:
I have a digital surface model of a region and the coordinates and height of an observer point and need to do a 3d - visibility analysis for this point.

However, as a result I would need a threedimensional 'visibility cone' (probably a multipart feature?) which depicts the areas that are visible, not only on the ground but also in the air (similar to the area confined by AZ1/AZ2, R1/R2, V1/V2 in this image.
I would also need to be able to set the various parameters like aperture angle, range etc.

Is doing something like this possible in ArcGIS?

Source: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/GUID-6E9426D5-6274-4CBA-9BE6-E71DE9598F7D-web.gif


texts: Error: Multiple Define with Arcgis Javascript API

texts:
I'm working on displaying GeoJSON data on ArcGIS map using ArcGIS javascript API.

Below is my code,

 <!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta http-equiv="X-UA-Compatible" content="IE=7,IE=9">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
    <title>Geoj</title>

    <link rel="stylesheet" href="https://js.arcgis.com/3.9/js/esri/css/esri.css">
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.1.1/css/bootstrap.min.css" >
    <link rel="stylesheet" href="https://esri.github.io/bootstrap-map-js/src/css/bootstrapmap.css">

I'm getting the following error:

 Error: multipleDefine
    at Error (native)
    at c (https://js.arcgis.com/3.9/init.js:15:378)
    at Ka (https://js.arcgis.com/3.9/init.js:36:286)
    at https://js.arcgis.com/3.9/init.js:37:216
    at f (https://js.arcgis.com/3.9/init.js:15:291)
    at Ya (https://js.arcgis.com/3.9/init.js:37:194)
    at h (https://js.arcgis.com/3.9/init.js:34:301)
    at HTMLScriptElement.<anonymous> (https://js.arcgis.com/3.9/init.js:39:128)

I'm not able to figure it out. Could anyone help me out on this?


texts: When does a QGIS .qix spatial index of a shapefile need to be re-created?

texts:
I understand the basics of spatial indexes but would like to understand when exactly a spatial index would need to be re-created or updated?

The reason I ask is because when a folder was moved to a new drive on my computer rendering slowed down considerably for the shapefiles contained within it. Before moving they all had a spatial index created. After the move and after I realised the rendering had slowed I created new ones not really knowing if it would help but it seemed to work.

I use QGIS 2.18.20