texts: spplot R and order of panels for Lattice plots

texts:
I have a shapefile with various attributes. I want to plot the shapefile using `spplot'.

I first read the shapefile

A <- readOGR('myshapefile.shp')

I am using the following command:

p <- spplot(A, c('lulc', 'lulc_2', 'lulc_3',
                 'lulc_4', 'lulc_5', 'lulc_6', 'lulc_7'),
            col.regions=myColors,
            col='transparent',
            names.attr = c(1998, 1999, 2002, 2003, 2009, 2014, 2017))

Surprisingly, the map order is from bottom to top and left to right in 3 rows.

I want them to appear from top to bottom and left to right order i.e. 1998 should appear first and 2017 should appear last.


texts:
We have a legacy tile server with tiles georeferenced to EPSG:4326 accessed as follows through openlayers:-

var OLDSERVICE =    new ol.layer.Tile({
                    title: 'OLD SERVICE',
                    type: 'base',
                    source: new ol.source.XYZ({
                                url: 'http://OLD_URL/maps/api/nls/{z}/{x}/{y}.jpg',
                                maxZoom: 15,
                                tilePixelRatio: 1})
                    });

The current base map is in EPSG:4326

We want to reproject these tiles to work on a base map with projection EPSG:27700

How would we do this?

I don't want to warp the tiles on the fly, would prefer to reprocess them for 27000. I am not very familiar with XYZ tile services as most of our stuff is done via a WMS service.


texts: Apply symbology to raster based on maximum raster value - arcpy

texts:
How can I apply different symbologies to rasters in a python script based on the maximum raster value?

The script loads one raster, applies symbology and outputs an MXD and PDF of the map, then moves on to the next raster.

I want to add a step which gets the maximum raster value and then selects one from three symbologies to be applied based on the maximum raster value.

Here is an extract of the script:

 # Get max raster value

timeMaxResult = arcpy.GetRasterProperties_management(Raster_lyr, "MINIMUM")
timeMax = timeMaxResult.getOutput(0)

# Apply symbology
if timeMax <= "31000":
    arcpy.ApplySymbologyFromLayer_management(Raster_lyr, symbology1)
elif timeMax > "31000" and timeMax < "160000":
    arcpy.ApplySymbologyFromLayer_management(Raster_lyr, symbology2)
else:
    arcpy.ApplySymbologyFromLayer_management(Raster_lyr, symbology3)


texts: How to make PROJ work via Anaconda in Google Colab?

texts:
I broke Google Colab with GDAL and Anaconda, how do I repair it?

So hear my story:
In a project I need to use the RSGISLib Python library (link) developed for object-based analysis of georeferenced images.

I also should work on google colaboratory for the sake of ease of access for all project members.

RSGISLib is only available on conda repository...

I managed to install Conda on Colab and get it running, then my notebook correctly finds the RSGISLib functions and executes them. So far so good. RSGISLib relies on GDAL for creation of spatial datasets. One of necessary things to get when creating such dataset is spatial projection, provided in GDAL by the PROJ library.

Once my notebook reaches a RSGISLib command, which needs PROJ at some point the entire Kernel/Linux/runtime (not sure really) of the Colab spectacularly crashes and I have to factory reset the runtime. The logs look like this, that's why I know it's a PROJ problem:

Jun 9, 2020, 12:09:38 PM    WARNING WARNING:root:kernel 11aacb5e-4fc8-4d62-8178-841f9a287a02 restarted
Jun 9, 2020, 12:09:38 PM    INFO    KernelRestarter: restarting kernel (1/5), keep random ports
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_create_from_database: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_create_from_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_as_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_as_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: pj_obj_create: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_create_from_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_create_from_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_as_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: pj_obj_create: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: createGeodeticReferenceFrame: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: proj_as_wkt: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: createGeodeticReferenceFrame: Open of /usr/condaFolder/share/proj failed
Jun 9, 2020, 12:09:37 PM    WARNING ERROR 1: PROJ: pj_obj_create: Open of /usr/condaFolder/share/proj failed

Nothing nice, as you can guess.
As I googled, I learned, that the reason is that PROJ when run from Anaconda requires an environment variable to be passed to it somehow for some reason. I thought "fine, I'll just:"

!conda activate base

Oh no, I wont:

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

    $ conda init <SHELL_NAME>

Currently supported shells are:
  - bash
  - fish
  - tcsh
  - xonsh
  - zsh
  - powershell

See 'conda init --help' for more information and options.

IMPORTANT: You may need to close and restart your shell after running 'conda init'.

!conda init bash doesn't work, it just doesn't change a thing. Funnily enough, when asked

!/usr/condaFolder/bin/conda info --envs

my on-line Jupyter gladly responds, that:

# conda environments:
#
base                  *  /usr/condaFolder

so technically I'm working in an active environment and there shouldn't be a problem... Yet there is. Currently I gave up on Colab and moved my notebook to a Conda environment on the local drive, but if you know how to get these things running please let me know.


texts: Custom python tool setting table field as a parameter?

texts:
I am somewhat new to Python, but I'm trying to create a custom python tool with 3 parameters, a workspace, a table and a Field from that table. I have the parameters set up like this:

def getParameterInfo(self):
    """Define parameter definitions"""
    parm0 = arcpy.Parameter(
    displayName = 'SSURGO geodatabase',
    name = 'geodatabase',
    datatype = "DEWorkspace",
    parameterType = "Required",
    direction = "Input")

    parm0.filter.list = ["Local Database"]

    parm1 = arcpy.Parameter(
    displayName = "Table",
    name = "Table",
    datatype = "DETable",
    parameterType = "Required",
    direction = "Input")

    parm2 = arcpy.Parameter(
    displayName = "Fields",
    name = "Fields",
    datatype = "Field",
    parameterType = "Required",
    direction = "Input",
    multiValue = True)

    parameters = [parm0, parm1, parm2]
    return parameters

I would like to create code in the updateParameters section that will allow the user to select a table and select a field from the fields within the table. Is there any way to do this?


texts: Extracting Sub-National admin boundaries or regions from GRUMPv1 or GPW v4

texts:
GRUMP and Gridded Population of the World (GPW), v4 provides population and demographic values for sub-national levels. However, I am unable to extract the polygons or boundaries from their datasets. They are a mix of admin3, admin2 and urban area extents so getting the boundaries from within these datasets is essential. More so, for African continent there is'nt any source, to my knowledge, covering entire continent at admin level 3.

How can I extract or get to the source polygons for which GRUMP or GPW data is calculated?


texts: Openlayers 4 - Load WMS image source require authentication

texts:
I try to load WMS image layer with openlayers 4.6 and angular 5, the code is:

const syr_layer = new ol_layer_Image({
  source: new ol_source_ImageWMS({
      url: 'serverurl', crossOrigin: 'anonymous', serverType: 'geoserver',
      params: { 'LAYERS': 'tst:syr'},
      projection: 'EPSG:4326'
    });
});

But it throws an error:

  GET (myserverurl) 401 (An Authentication object was not found in the
  SecurityContext)

How can I send authentication header with the GET request sent by openlayers?


texts: Creating Triangulated Irregular Network (TIN) in QGIS?

texts: I want to create a TIN  and view it in 3D in qgis. The tin am creating using interpolation plugin is a raster TIN but i want the normal TIN(vector TIN).


texts: Unable to get latitude and longitude of the US map from TopoJSON's us-atlas repository when using GeoPandas

texts:
To do some geospatial analysis I am using the USA counties map from topojson/us-atlas repository
After loading the link https://cdn.jsdelivr.net/npm/us-atlas@3/counties-10m.json in geopandas I get the geodataframe correctly.
import geopandas as gpd
counties = gpd.read_file('https://cdn.jsdelivr.net/npm/us-atlas@3/counties-10m.json')

Then I just do the following to get the longitude and latitude -
counties['lon'] = counties['geometry'].centroid.x
counties['lat'] = counties['geometry'].centroid.y

But I get the error -
IndexError: list index out of range

Only the centroid attribute seems to give the output :
counties['geometry'].centroid gives -
0       POINT (-113.75823 35.70499)
1        POINT (-90.40498 30.62627)
2       POINT (-114.13825 43.00235)
3        POINT (-96.40159 47.77421)
4        POINT (-97.24835 46.93312)

Any idea how to get the longitude and latitude of the counties / states from the us-atlas file?


texts: TMS in UTM grid with MapProxy & OpenLayers 4?

texts:
I am trying to cache a GeoServer layer with MapProxy and display it on an OpenLayers4 map. I created a new MapProxy grid in UTM32N and in the MapProxy preview the layer is shown nicely, but when I add it to my ol map, it won't display at all. Requests are sent but return empty tiles.

I have tried several grid definitions, two are shown below ("utm32n" and "germany") in the mapproxy config:

services:
  demo:
  tms:
    use_grid_names: true
    origin: 'sw'

layers:
  - name: mylayer
    title: mytitle
    sources: [mysource]
caches:
  mylayer:
    grids: [webmercator, germany]
    sources: [mysource]

sources:
  mylayer:
    type: wms
    req:
      url: http://mygeoserver/ws/wms
      layers: mylayer
      transparent: true

grids:
    utm32n:
        srs: EPSG:25832
        bbox: [239323.44497139292,4290144.074117256,760676.555028607,9503675.174689397]
        bbox_srs: EPSG:25832
        origin: sw
        min_res: 156367.7919628329
        num_levels: 18
    germany:
      srs: 'EPSG:25832'
      bbox: [6, 47.3, 15.1, 55]
      bbox_srs: 'EPSG:4326'
      origin: 'sw'

In my ol map, I have tried to define the layer as ol.source.TileImage:

var extent = ol.proj.transformExtent([6, 47.3, 15.1, 55], "EPSG:4326", epsg_code);
var size = ol.extent.getWidth(extent) / 256;
var no_zoomlevels = map.getView().getMaxZoom() - map.getView().getMinZoom() + 1;
var resolutions = new Array(no_zoomlevels);
for(var i=map.getView().getMinZoom(); i<=map.getView().getMaxZoom(); i++){
    resolutions[i] = map.getView().getResolutionForZoom(i);
}
var layer = new ol.layer.Tile({
    name: name,
    source: new ol.source.TileImage({
        projection: epsg_code,
        tileGrid: new ol.tilegrid.TileGrid({
            resolutions: resolutions,
            extent: extent,
            origin: [extent[0], extent[1]],
            minZoom: map.getView().getMinZoom(),
        }),
        tileUrlFunction: function(tileCoord, pixelRatio, projection){
            var z = tileCoord[0];
            var x = tileCoord[1];
            var y = - tileCoord[2];
            var url = 'http://mymapproxy/tms/1.0.0/mylayer/germany/' + z + '/' + x + '/' + y + '.png';
            return url;
        }
    }),
});

as well as ol.source.XYZ:

...
source: new ol.source.XYZ({
    url: 'http://mymapproxy/tms/1.0.0/' + name + '/germany/{z}/{x}/{-y}.png',
    projection: epsg_code,
    minZoom: map.getView().getMinZoom(),
    maxZoom: map.getView().getMaxZoom(),
}),
...

It seems like the ol grid and the mapproxy grid are not the same, but how do I match them?


texts: Handling geoExt popup when one or more layer fails to load

texts:
GeoExt info Popup is very easy to use but when one or more layer is missing or fails to load
Avoiding Pink Tyles. The popup is not filled with the values. While I can see that request is sent to the server and response is there but popup is blank. Is there any way to handle this condition without building custom component?

OpenLayers.Util.onImageLoadErrorColor = 1;
            OpenLayers.Util.onImageLoadErrorColor = "transparent";
            OpenLayers.Util.onImageLoadError = function(){
                 this.src = "/img/blank.gif";
            };

controls.push(new OpenLayers.Control.WMSGetFeatureInfo({
                autoActivate: true,
                infoFormat: "application/vnd.ogc.gml",
                maxFeatures: 3,
                eventListeners: {
                    "getfeatureinfo": function(e) {
                        var items = [];
                        Ext.each(e.features, function(feature) {
                            items.push({
                                xtype: "propertygrid",
                                title: feature.fid,
                                source: feature.attributes
                            });
                        });

                        new GeoExt.Popup({
                            title: "Info",
                            width: 200,
                            height: 200,
                            layout: "accordion",
                            map: app.mapPanel,
                            location: e.xy,
                            maximizable: true,
                            collapsible: true,
                            items: items
                        }).show();

                    }
                }
            }));


texts: GRASS GIS 7.8 trigometric function of vector attribtue

texts:
I need to calculate the sin and cos of a vector attribute column, but these functions do not appear to be available. In Python I've tried using v.db.update:
gscript.run_command('v.db.update', map='orientation', column='Z', query_column='COS(AZIMUTH)')

but this was unsuccessful. Does anyone have recommendations?


texts: Error 000539 running multiple con expressions in raster calculator in ModelBuilder?

texts:
I am using Raster Calcuator in ModelBuilder with ArcMap 10.5.

I used Get Raster Properties to extract Standard Deviation and Mean from a previously calculated raster.

I would like to use Raster Calc again to create a raster where raster = (mean + 3*(SD)) or  (mean -3*(SD)). I attempted to use a Con statement with a boolean or (|) operator. Con statement I am using is as follows:

Con(("difference_mean_normalized_.tif"  <  (float(%Mean%) - 3 * (float(%Standard Deviation%)))) | ("difference_mean_normalized_.tif"  > (float(%Mean%) + 3 * (float(%Standard Deviation%)))), 0, 1)

Where difference_mean_normalized_.tif is my output raster.

This fails and I am getting a 000539 Error:

Start Time: Wed Oct 17 07:12:01 2018
Con(("difference_mean_normalized_.tif"  <  (float(-4.6467203950133E-06) - 3 * (float(0.036299477806278)))) | ("difference_mean_normalized_.tif"  > (float(-4.6467203950133E-06) + 3 * (float(0.036299477806278)))), 0, 1)
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 5, in rcexec
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 263, in Con
    where_clause)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\sa\Functions.py", line 250, in Wrapper
    ["IfThenElse", in_conditional_raster, in_true_raster_or_constant, in_false_raster_or_constant])
RuntimeError: ERROR 999998: Unexpected Error.

Failed to execute (Raster Calculator (3)).

Some notes: I double-click or drag the variables in from the list in Raster Calc, and "float" is automatically added to the expression for each variable. I manually type in the output raster name. I have tried wrapping the raster name with the % symbol as well, although this is not shown in the screenshot. Also, my original bands were integer, and I had to use the float() operator to convert.

How do I use Raster Calc in ModelBuilder to create a raster +/- 3SD from mean?


texts: Assessing if two variables (measured at point locations) covary in space?

texts:
I have a number of (point) locations at which two variables are measured. The two variables proves correlated (when using Spearman's r; the direction of the correlation is irrelevant here).

How can we assess (graphically and/or statistically) if the two variables covary in space?

For instance, if the two variables tend to increase as one moves (say) to the north, or something along that line.


texts: Cropping maps (Raster source map and vector map as boundary)?

texts:
I have DEM named "HD_DEM" (asc file) and vector map for lakes that named "Lakes". I want to crop lakes from HD_DEM. By following instructions of this link https://grasswiki.osgeo.org/wiki/Cropping_maps ,I tried to crop.

v.to.rast in=Lakes@TEST out=Lakes_MASK3@TEST use=val value=1

and

r.mask vector=Lakes@TEST

the both gave me strange output. (Grey: Lakes, Yellows: output)


texts: Connecting to SQL Server with ArcGIS Desktop?

texts:
I am trying to connect to views in a SQL Server DB in Arcatalog / Arcmap.
I am using the "Database Connections\Add Database Connection" method in catalog.
I can connect fine but cannot see / preview data in catalog.
It's not an SDE database, but another company DB.

I have a suspicion it’s a driver issue (http://desktop.arcgis.com/en/system-requirements/latest/database-requirements-sqlserver.htm) but I get lost about what to use and how to get these files / install them (32 bit or 64 bit etc.)

SQL specs:
Microsoft SQL Server 2016 (RTM) - 13.0.1601.5 (X64)
Standard Edition (64-bit) on Windows Server 2012 R2 Standard 6.3  (Build 9600: ) (Hypervisor)

Arc specs:
ArcgGIS 10.5.0.6491

Windows 10


texts: Geodjango render_to_kml Multipolygons - Can be fast or faster as WMS?

texts:
I am learning now geodjango. I made basic web app with postgis database. With openlayers  iam displaying basic wms map and I want to add other layers from my postgis database. I use render_to_kml function and then I visualize this kml with openlayers. The problem is that I can only display point layers instead of multipolygons...I tried to make another kml template bud I misunderstand it and there is no good documentation... I paste here my files ( view, template, urls) and please give me advice if you can,,,where I have mistake ( I think it will be something with KML templates but i dont know.)

this is my views.py
this is my html template
this is my urls.py

This results into application with 3 layers (basic map and 2layers from my postgis database which are equivalent (both are point feature....Polygons feature is only copy of Point feature...Thats my problem))...
Anyone has a solution for me?


texts: How to prevent running part of a model when iterator is used?

texts:
My process is to run Buffer on two feature layers, intersect the buffers and then Clip intersected buffers with polygon object from Iterate Feature Selection iterator. Here's my model:

The model runs Buffers and Intersect firstly, then iterate and create Single Object, and lastly Clip features. But in the next iteration it again runs Buffers and Intersect, what is not needed.

Is it possible to prevent running Buffers and Intersect in each iteration?
In other words, how to create Intersected buffers only once and then reuse it in each iteration?

EDIT:
Following blah's suggestion I managed to create desirable model. I enclose it, because it's slightly different than web-help example:


texts: How to remove projection in terminal from .shp file with valid .prj available

texts:
I have .prj, .shp, .sbn, .sbx, .shx, .dbf files where the content of the .prj file is the following:

PROJCS["S-JTSK_Krovak_East_North",GEOGCS["GCS_S_JTSK",DATUM["D_S_JTSK",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Krovak"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Pseudo_Standard_Parallel_1",78.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Azimuth",30.28813975277778],PARAMETER["Longitude_Of_Center",24.83333333333333],PARAMETER["Latitude_Of_Center",49.5],PARAMETER["X_Scale",-1.0],PARAMETER["Y_Scale",1.0],PARAMETER["XY_Plane_Rotation",90.0],UNIT["Meter",1.0]]

I would like to generate a topojson file which is not projected so I will be able to define the projection in d3. How should I parametrize topojson in Terminal to remove projection?


texts:
I'm putting together a large map consisting of some big datasets, and unfortunately deleted some (around 25) polygons from a key layer when editing a smaller version of the map for a different request. The data for these polygons is still in the attribute table. I need these polygons back and plotted on the map. I don't want to remove the layer and add it back again as I spent a few hours inputting extra data to the attribute table, and I would have to do this again.

Is there a way I can replot these deleted polygons and keep the data I added in the attribute table?


texts: Fixing line breaks in cartoDB infowindow?

texts: I am having trouble with my cartoDB infowindow, the text breaks no matter if I resize the window. Please, take a look at the screen shot below. How can I resolve the issue? Where should I work, dataset or carto map editor?


texts: IDE (PyCharm) settings for developing with Processing in Linux (Ubuntu)

texts:
How to setup Python interpreter in IDE to develop using Processing module? I need this not for running from IDE, but to have tips about classes, methods, arguments etc.

IDE: PyCharm 2016.3 (but example in any IDE will be helpful).
OS: Ubuntu 16.04 (any Linux example will be OK too).

In my current PyCharm settings project interpreter located at: /usr/bin/python2.7 and imports from module QGIS works fine, but not from processing (red curved underlines):

I've discovered, that processing located at /usr/share/qgis/python/plugins/processing but don't know what to do with this info.

And second. If anybody knows how to setup python virtual environment with all modules, that needed for QGIS plug-in developing, please tell. It may be helpful too. Especially if I’ll develop standalone app in future.


texts: Importing spreadsheet CSV file with coordinates into QGIS

texts: I watched a simple video tutorial on how to accomplish it, but nothing shows up on my project when I add the layer. Even if I try zoom to layer, the screen just shows blue. I'm not sure if it has something to do with my Excel, computer, or program. I really need help as I am new to QGIS.


texts: Add a specific Domain to Multiple Feature CLasses at once (batch)

texts:
I have a domain named "Life Cycle" with coded- "active" or "inactive". I need to apply this domain to over 150+ Feature Classes (in a SDE if that matters) via ArcCatalog. I know how to do this manually by adding the domain to each feature class field; but it will take all a very long time to do this for 150+ Feature Classes and I know I will have more domains similar to this that I will have to add eventually to these Feature CLasses down the road.

I want to know if there is a tool/add-in or script someone can enlighten me on automating this process. Or allow me to add the domain to the specific Field "Life Cycle" for all desired Feature Classes at once.


texts: Creating population dot density map at the block group level in QGIS

texts:
I'm new to QGIS, always used ArcMap, and I'm struggling to find out how to best display population density at the block group level for a county (the most granular level of detail for census geography). I have the shapefile for the block groups, and I've joined the csv containing my population data per block group to the shapefile, but now I'm at a loss on:

how to best display the population density
the method for executing the task in QGIS

Most of the block groups have populations ranging from 25-100 people. Would it be best to have 1 dot equal 25 people per census block?  How do I create this type of map in QGIS? Any suggestions?


texts: How to rasterize a river vector file and merge it to another raster layer?

texts:
I got a river vector file.
I want to rasterize it so that the river is 2 pixels wide.
How do I achieve that in QGIS 3.2.3?

When I try to use vector to raster, I get this and no file is created:

Algorithm 'Rasterize (vector to raster)' starting…
Input parameters:
{ 'BURN' : 0, 'DATA_TYPE' : 5, 'EXTENT' : '413885.10813709395,678750.4296324956,6754403.012272906,6909638.055526119 [EPSG:32635]', 'FIELD' : None, 'HEIGHT' : 0, 'INIT' : 0, 'INPUT' : 'C:\\Users\\vvaris\\Desktop\\Uoma\\Uoma10.shp', 'INVERT' : False, 'NODATA' : 0, 'OPTIONS' : '', 'OUTPUT' : 'C:/Users/vvaris/Desktop/asd.tif', 'UNITS' : 0, 'WIDTH' : 0 }

GDAL command:

gdal_rasterize -l Uoma10 -burn 0.0 -ts 0.0 0.0 -a_nodata 0.0 -te 413885.10813709395 6754403.012272906 678750.4296324956 6909638.055526119 -ot Float32 -of GTiff C:\Users\vvaris\Desktop\Uoma\Uoma10.shp C:/Users/vvaris/Desktop/asd.tif

GDAL command output:

ERROR 1: Wrong value for -outsize parameter.

Usage:

gdal_rasterize [-b band]* [-i] [-at]

{[-burn value]* | [-a attribute_name] | [-3d]} [-add]

[-l layername]* [-where expression] [-sql select_statement]

[-dialect dialect] [-of format] [-a_srs srs_def]

[-co "NAME=VALUE"]* [-a_nodata value] [-init value]*

[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]

[-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

CInt16/CInt32/CFloat32/CFloat64}] [-q]

<src_datasource> <dst_filename>

Execution completed in 0.13 seconds

Results:

{'OUTPUT': <QgsProcessingOutputLayerDefinition {'sink':C:/Users/vvaris/Desktop/asd.tif, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
Algorithm 'Rasterize (vector to raster)' finished


texts: How to drap a texture over a surface with python and mayavi?

texts:
I want to add a picture over a surface that I made with mayavi and gdal.

I found a website (http://geoexamples.blogspot.ch/2014/02/3d-terrain-visualization-with-python.html) which explains how to do it but it doesn't work for me. I also saw that there were many unanswered questions about it on the forum.

Here is my program: Why isn't it working ?

import gdal
from mayavi import mlab
from tvtk.api import tvtk
import PIL
import numpy as np

# maido is the name of a mountain
# tipe is the name of a french school project

# 1) opening maido geotiff as an array
maido = gdal.Open('dem_maido_tipe.tif')
dem_maido = maido.ReadAsArray()

# 2) transformation of coordinates
columns = maido.RasterXSize
rows = maido.RasterYSize
gt = maido.GetGeoTransform()
ndv = maido.GetRasterBand(1).GetNoDataValue()

x = (columns * gt[1]) + gt[0]
y = (rows * gt[5]) + gt[3]

X = np.arange(gt[0], x, gt[1])
Y = np.arange(gt[3], y, gt[5])

# 3) creation of a simple grid without interpolation
X, Y = np.meshgrid(X, Y)

#Mayavi requires col, row ordering. GDAL reads in row, col (i.e y, x) order
dem_maido = np.rollaxis(dem_maido,0,2)
X = np.rollaxis(X,0,2)
Y = np.rollaxis(Y,0,2)

# 4) deleting the "no data" values
dem_maido = dem_maido.astype(np.float32)
dem_maido[dem_maido == ndv] = np.nan #if it's NaN, mayavi will interpolate

# delete the last column
dem_maido = np.delete(dem_maido, len(dem_maido)-1, axis = 0)
X = np.delete(X, len(X)-1, axis = 0)
Y = np.delete(Y, len(Y)-1, axis = 0)

# delete the last row
dem_maido = np.delete(dem_maido, len(dem_maido[0])-1, axis = 1)
X = np.delete(X, len(X[0])-1, axis = 1)
Y = np.delete(Y, len(Y[0])-1, axis = 1)

# 5) ploting the raster and draping a picture of maido over the surface
img_filename = "maido.jpg"
img = PIL.Image.open(img_filename)
bmp = tvtk.JPEGReader(file_name=img_filename)

t = tvtk.Texture()
t.interpolate=0
t.set_input(0, bmp.get_output())

mlab.figure(size=(1000, 1000))

surf = mlab.surf(X, Y, dem_maido, color=(1,1,1), warp_scale="auto")

surf.actor.actor.mapper.scalar_visibility = False
surf.actor.enable_texture = True
surf.actor.tcoord_generator_mode = 'plane'

mlab.show()

On the left you have what I get and on the right you have the picture I want to use.

Here is the GeoTiff file and the picture that I used (it is uploaded on my own google account don't worry): https://drive.google.com/drive/folders/0B7P95aWmH4DUbzJCaVZQM3Y0Y1U?usp=sharing


texts: Find distance between multiple discrete sets of points QGIS3.8

texts:
No similar question seems to relate to my issue but I desperately need help with calculating distances.
I have a CSV with the residential postcode of 200 000 patients. For each patient, I also have the postcode of their respective doctor. I need to find out how far, on average, patients are travelling to see their doctor.

So I've created centroids for all my postcodes thinking that a point layer would be easier to work with, but now I have no idea what to do. Old threads speak of plugins like Points2one that no longer exist (to my knowledge).
I'm assuming I need to produce a line between all individual sets of points, then calculate the length of the lines and get the average.
Help!


texts: Update layer with user input in Earth Engine

texts:
My goal is to request a threshold value from user and then use it to recompute a binary mask for NDVI raster. But when I change a value in the ui.Textbox ndvi_textbox it throws an error: Invalid argument specified for ee.Number(): 0.2

Here's my code:

var ndvi = ee.ImageCollection("LANDSAT/LC8_L1T_32DAY_NDVI")
  .first();

var threshold = 0.1;

function set_threshold(thrsh) {
  var thresholded = ndvi.gt(ee.Number(thrsh));
  Map.addLayer(thresholded);
}

var ndvi_textbox = ui.Textbox({
  value: threshold,
  onChange: set_threshold,
  style: {position: 'top-center'}
});

ui.root.setLayout(ui.Panel.Layout.absolute());
ui.root.add(ndvi_textbox);

set_threshold(threshold);

Map.setCenter(6.969, 47.448, 9);


texts: Display coordinates for a specific projection in ArcMap?

texts:
Below are the available options I can find in the settings:

How can I pick a specific projection such as Web-Mercator (EPSG:3857)?


texts: MapServer - map rendered at the wrong angle

texts:
I've just started using MapServer to check the topology in my database, and whenever I render a map I have to flip the image vertically and then rotate it quite a bit to get it facing North. I've set the projections to WGS84 but it doesn't seem to have any effect. Any ideas what I'm doing wrong?

The way I'm calling MapServer

http://localhost/cgi-bin/mapserv?map=test.map&layers=all&map_imagetype=JPEG&mode=map&mapext=52.8055+-2.0924+52.8175+-2.0785&map_size=1000+1000

The map file can be found at http://pastie.org/pastes/1647658/text

Here's an image showing how the map gets rendered and also what angle it should be at..

Also, when I try to overlay the image onto Google Maps I can't get any of the roads to line up with the Google ones. I know that my data is correct as I've checked the lon/lat values in the database against Google and they match up perfectly.. it's only when I output via MapServer that it doesn't match up.

Any thoughts would be appreciated!


texts: Change palette color in Google Earth Engine

texts:
I've modified this code: https://code.earthengine.google.com/?scriptPath=Examples:User+Interface/Split+Panel but I can't change the maps color pallet

This is the color pallet that I want to use:

var colorizedVis = {
min: -1.0,
max: 1.0,
palette: [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
'66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
'012E01', '011D01', '011301'
],
};

My code can be found here:

https://code.earthengine.google.com/180bc2a78ec92b433ea7c369bb2e249c

or here:

Demonstrates before/after imagery comparison with a variety of dates.

Configure the imagery

  var images = {

  'L5_01_98' : ee.Image('LANDSAT/LT05/C01/T1_8DAY_EVI/19990101'),

  '5_04_08' : ee.Image('LANDSAT/LT05/C01/T1_8DAY_EVI/20080414'),
  '5_01_11' : ee.Image('LANDSAT/LT05/C01/T1_8DAY_EVI/20110117'),

};

  var colorizedVis = {
    min: 0.0,
    max: 1.0,
    palette: [
      'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
      '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
      '012E01', '011D01', '011301'
    ],
  };

Set up the maps and control widgets

Create the left map, and have it display layer 0.

  var leftMap = ui.Map();
  leftMap.setControlVisibility(false);
  var leftSelector = addLayerSelector(leftMap, 0, 'top-left');

Create the right map, and have it display layer 1.

  var rightMap = ui.Map();
  rightMap.setControlVisibility(false);
  var rightSelector = addLayerSelector(rightMap, 1, 'top-right');

Adds a layer selection widget to the given map, to allow users to change
which image is displayed in the associated map.

  function addLayerSelector(mapToChange, defaultValue, position) {
  var label = ui.Label('Choose an image to visualize');

This function changes the given map to show the selected image.

  function updateMap(selection) {
  mapToChange.layers().set(0, ui.Map.Layer(images[selection]));
  }

Configure a selection dropdown to allow the user to choose between images, and set the map to update when a user makes a selection.

  var select = ui.Select({items: Object.keys(images), onChange: updateMap});
  select.setValue(Object.keys(images)[defaultValue], true);

  var controlPanel =
  ui.Panel({widgets: [label, select], style: {position: position}});

  mapToChange.add(controlPanel);
  }

Tie everything together

Create a SplitPanel to hold the adjacent, linked maps.

  var splitPanel = ui.SplitPanel({
    firstPanel: leftMap,
    secondPanel: rightMap,
    wipe: true,
    style: {stretch: 'both'}
  });

Set the SplitPanel as the only thing in the UI root.

  ui.root.widgets().reset([splitPanel]);
  var linker = ui.Map.Linker([leftMap, rightMap]);
  leftMap.setCenter(-65.49522105694524,-24.717770902233013, 13);


texts: How to change spelling of scalebar label using arcpy?

texts:
I have come to realise that hundreds of my MXD files have the scalebar units label spelled as "Meters" instead of "Metres".  I suspect this is something I am going to have to do manually (I don't know ArcObjects).

Looking at the MAPSURROUND_ELEMENT class documentation, it doesn't appear that this is something I can change using arcpy.

Is there a way to change the spelling of the units label in my scalebar using arcpy?


texts:
I'm trying to use arcpy.ExtractValuesToPoints() to get raster values to a points feature class for several raster files (8760). The goal is to get a time-slice view of how the values at certain locations change over time.

I'm fairly well versed in Python as well as the built-in datetime module. The trouble that I'm having is saving the date-time values to a table in ArcGIS.

(I've gotten a similar approach to work using arcpy.da.InsertCursor() but this approach is very process-intensive and I've concluded that it isn't a viable approach for what I'm trying to do. Plus I'd rather not re-invent the wheel and use the ArcGIS function: ExtractValuesToPoints() instead.)

To give a little background, I'm using dates that require datetime precision to the hour so I've elected to use a File Geodatabase to accommodate that feature (since shapefiles can only handle the date portion).

I'm trying to solve for a single case as I can create a solution for iterating over the 8760 files that I'm working with.

So far this is what I have:

import arcpy
import os

GRID_Table = "E:\HRRR.gdb\GRID_Table"               # GRID_Table is an ArcGIS table
SPP_LMP_Points = "E:\HRRR.gdb\SPP\SPP_LMP_Points"   # that contains datetime and
out_path = "C:\scratch.gdb"                         # raster_path pairs.

row1 = []                                           # Get date_time and raster_path
for row in arcpy.SearchCursor(GRID_Table):          # from every entry in GRID_Table
    row1.append((row.DateTime, row.RasterPath))     #

date_time = row1[0][0]                              # Set date_time to the first
raster_path = row1[0][1]                            # entry

out_path_shp = os.path.join(out_path, os.path.basename(raster_path))  # Generate
                                                                      # output path.

arcpy.sa.ExtractValuesToPoints(SPP_LMP_Points, raster_path, out_path_shp)

arcpy.AddField_management(out_path_shp, "DateTime", "DATE")

arcpy.CalculateField_management(out_path_shp, "DateTime",
                                "!DateTime! = {0}".format(date_time),
                                "PYTHON_9.3")

This all works fine and doesn't report any exceptions or errors. The trouble that I'm having is that the call to arcpy.CalculateField_management() doesn't populate the DateTime field. Instead I get a NULL value in the table.

Any thoughts?


texts:
What I'm trying to do is basically plot a rectangle on a runway that goes parallel with the runway.

When I start out everything works as expected:

But what I notice is that as I go further down the runway it varies significantly from the centerline as we see in the end of my rectangle:

This is not what I would expect, as I would expect the centerline to continue to match the runway centerline, or at least be closer.

My calculation for the points to plot is rather simplistic, as it is just the runway heading +/- the magnetic variation and it leaves me with this result. So my question is, is there anything else in my calculation I am missing in order to get this to be more precise?


texts: Displaying Image for built-in command added to Python Addin toolbar?

texts:
I'm trying to create a Custom Python Addin for ArcMap, which contains a custom tool, and also needs to have the standard 'Save As' button.

After referring to https://blogs.esri.com/esri/arcgis/2012/08/27/python-add-ins-how-to-add-built-in-commands-to-your-custom-toolbar-or-menu/,

I've added a Tag <Button refID="esriArcMapUI.MxFileMenuItem" subtype="4" />

This shows the button in the toolbar, but only as text.

How can I get it to show up as an Image? I know that the 'Save As' command doesn't have an Image in the File Menu, so I can supply a new image for it.


texts: Using kriging model to predict values from image pixel data

texts:
I have 1000 images that look like this:

I have the image data 64x128 saved for each image along with each image's compliance value (i.e. 107.33, 207.45, etc.)
Is it possible to create a kriging model that can predict the compliance values of new images?


texts: GXP Hidden layer switcher - composer.html

texts:
I am working with GeoServer, GeoExplorer, GXP. I am creating a WPS client and need to generate a vector layer that serves to make operations on vector data, such as buffer, split polygon.  I know how to create the vector layer but I want that is not on the menu layers GeoExplorer, I want this hidden.

I try modify this script wpsdemo.js from (http://suite.opengeo.org/docs/latest/processing/wpsclient/index.html ) where:

In wpsdemo.js replace:

....
target.on('ready', function() {
            // Get a reference to the vector layer from app.js
            this.layer = target.getLayerRecordFromMap({
                name: 'sketch',
                source: 'ol'
            }).getLayer();
            // Some defaults
....

for code that does the following:

create dynamically  OpenlayerLayer.Vector (with name: ''sketch', source: 'ol', displayInLayerSwitcher: false)  <---- I know how to do it
add this layer to the map   <---- I don't know (my map is defined in composer.html)
GIS operations do but with the hidden layer in the layers menu switch

I try the following settings (in the composer.html) using the property displayInLayerSwitcher: false but does not work :

{
   // A vector layer to display our geometries and processing results
   source: "ol",
   name: "Res",
   type: "OpenLayers.Layer.Vector",
   fixed: true,
   visibility: false,
   selected: false,
   displayInLayerSwitcher: false
}

In what way I can hide the layer in the layer switcher of GeoExplorer?


texts: Pre-Processing for multi-sensor Time series analyses based on Landsat Level 2 imagery

texts:
I downloaded TM, ETM+ and OLI Surface Reflectance Higher-Level Data Products (atmospheric corrected) based on L1TP Landsat images (radiometrically calibrated, orthorectified). I need to include images from those three different sensors to get a complete time series for my study area.

Based on the description given by the USGS and a very helpful paper i read regarding pre-processing (https://www.researchgate.net/publication/312202874_A_survival_guide_to_Landsat_preprocessing) i conclude that most of the usual pre-processing steps are unnecessary in this case by simply using the Level 2/Surface reflectance product, except for:

1) Cloud masking
2) Radiometric Normalization to account for radiometric differences between the three sensors and their different atmospheric correction methods (LEDAPS and LaSRC)

Now my question:

Will a "Relative Radiometric Normalization" (f.ex. PIF or the relnorm() function in R) help me to get rid of the aforementioned differences and make a time series analyses based on these different sensors possible?


texts:
I have an ImageCollection of Sentinel2 as follows.

var Sentinel2 = ee.ImageCollection('COPERNICUS/S2')
                    .filterDate('2018-01-01','2018-12-31')
                    .filterBounds(AOI);

From above ImageCollection (Sentinel2) of 100 (suppose) images, I want to randomly sample 10 and output this as ImageCollection of 10 images.

How can I do it?


texts: FutureWarning: Sorting because non-concatenation axis is not aligned. Converting a Google Earth Engine Dataset into a Pandas/Geopandas Dataframe

texts:
I am trying to get a pandas or geopandas dataframe returned from a Google Earth Engine FeatureCollection for specific point data. Basically like a spatial intersect where all values from the feature collection are returned for each point.
I found this code example below but it returns with a warning that the nonconcat axis is not aligned. My end goal is to perform Machine Learning so I need to make sure that the point data is perfectly aligned with the results returned.
import ee
import pandas as pd

ee.Initialize()
districts = points_features_collections

#selecting images with a data/time stamp and averaging results
chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").select('precipitation').filterDate('2018-01-01', '2018-12-31').mean()
evi = ee.ImageCollection('MODIS/006/MOD13A2').select('EVI').filterDate('2018-01-01', '2018-12-31').mean()

composite = ee.Image.cat([chirps, evi]) #concat multiple images into one

#stores multiple regions as one in a feature collection
dist_stats = composite.reduceRegions(districts, 'mean', 7000)

dist_stats = dist_stats.select(['precipitation', 'EVI'], ['tot_rainfall', 'med_evi'], retainGeometry=False).getInfo()

dtstats_df = pd.DataFrame()
for dist in dist_stats['features']:
    df = pd.DataFrame([dist['properties']],columns=dist['properties'].keys())
    dtstats_df = pd.concat([dtstats_df, df], axis=0)

dtstats_df = dtstats_df[['tot_rainfall', 'med_evi']]

dtstats_df


texts:
The Google hybrid layer from the OpenLayers Plugin is no longer working for any of our maps. It's gone completely blank, and uninstalling and reinstalling the plugin hasn't resolved the issue.

Has anyone else experienced this and knows how to address the problem?

Sorry for the vagueness of the question, or if this is a duplication, but I'm still something of a novice with QGIS!


texts: How to calculate the shared area between a raster and a shapefile in ArcGIS 10?

texts:
I have a raster file with 5 classifications developed from habitat modelling algorithms. I'd like to know how much of 2 classes of the raster is covered by the network of protected areas (Shapefile).

Could any one advise me how to do so?


texts: QGIS 2.16 GRASS geoalgorithms don't work?

texts:
I have problem again almost like version QGIS 2.14. I did not solve the problem.
QGIS 2.14 GRASS geoalgorithms don't work?

I install C:\MinGW and now QGIS grass geoalgorithms don't work?
There is 300 geoalgorithms and I have only 148, what is the problem?
Is there a solution??
PC Windows 10 x64.

I need other grass geoalgorithms like v.net and others


texts: Open QGIS prococessing algorithm dialog with Python

texts:
I would like to programmatically show a processing algorithm dialog from Python. I gived a look to processing toolbox and I found that this is possible importing and instantiating the class AlgorithmDialog gui/processingToolbox.py

from processing.gui.AlgorithmDialog import AlgorithmDialog
from qgis.core import QgsApplication

alg = QgsApplication.processingRegistry().algorithmById('qgis:extractbyattribute')
dlg = AlgorithmDialog(alg, False, iface.mainWindow())
dlg.show()
dlg.exec_()

This opens algorithm dialog and let me perform processing computations but once closed the dialog window QGIS become unstable and crash without any message interacting with the user interface.

What am I doing wrong? Is there a method to do this in other way?


texts:
Are there any advantages to the use of GRASS Tools (Creating a new mapset and using GRASS modules) instead of using GRASS algorithms available in the processing toolbar with QGIS Data?

I think there might be a speed advantage because that way the layers are already in GRASS format and can be directly treated. I'd like to know if my guess is right. Did somebody compare?

I think there might also be more options (flags) when using GRASS tools instead of GRASS algorithms in the Processing Toolbox. Are there?

I also think that there might be more algorithms available or more up to date algorithms in GRASS tools than in the Processing Toolbox. Are there?


texts: How to I create a .prj file with the following projection parameters?

texts:
I have received shape files that are missing .prj files. The documentation for the shape files describes the projection system as Lambert Conformal Conic and provides the following "Projection parameters".

Projection parameters :

Input
Projection geographic
units dd
parameters
output
projection Lambert
spheroid everest
units meters
parameters
23 09 00
28 48 00
90 00 00
26 00 00
2743185.699
914395.233
end

Can anyone tell me how I can translate these parameters into a usable .prj file?

As is probably clear from my question, I have a very weak (or no) grasp of of projections and spatial references.


texts: Processing Sentinel-3 SLSTR data using python via snappy API

texts:
I am trying to convert Sentinel-3 SLSTR S6 radiance data to a Geotiff file; However, I am finding converting the Sentinel-3 data which is in NCDF4 format to Geotiff format correctly very difficult. I am trying to use snappy API from Python. Since the lat-long information is provided in a separate file, my main problem is how to overlay geo-coordinates and radiance values from S6 and write it in a tiff file with the correct projection and resolution. Any help would be highly appreciated. The code that I am trying till now is below. If somebody could guide me in the process it would be very helpful

rad13 = p.getBand(‘S6_radiance_ao’)
w = rad13.getRasterWidth()
h = rad13.getRasterHeight()
rad13_data = np.zeros(w * h, np.float32) # numpy array of longitude
rad13.readPixels(0, 0, w, h, rad13_data)
lat = p.getBand(‘x_ao’)
long = p.getBand(‘y_ao’)
w1 = lat.getRasterWidth()
h1=long.getRasterHeight()
lat_data = np.zeros(w1 * h1, np.float32) # numpy array of latitude
lat.readPixels(0, 0, w, h, lat_data)
long_data = np.zeros(w1 * h1, np.float32)
long.readPixels(0, 0, w, h, long_data) # numpy array of longitude

###After I got the separate numpy arrays of lat-long and radiance values I have been trying to convert it into a geotiff file
xmin,ymin,xmax,ymax = [long_data.min(),lat_data.min(),long_data.max(),lat_data.max()
nrows,ncols = np.shape(rad13_data)
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName(‘GTiff’).Create(‘myraster.tif’,ncols, nrows, 1 ,gdal.GDT_Float32)
output_raster.SetGeoTransform(geotransform)
####I am getting a 0 value here so not able to proceed
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
###getting a 0 value here
output_raster.SetProjection( srs.ExportToWkt())
output_raster.GetRasterBand(1).WriteArray(rad_array)


texts: Need help regarding Hard disk space for Amazon EC 2 for ArcGIS server 10.1

texts:
It is regarding ArcGIS Server for Amazon EC2 instance.

We have ArcGIS Server 10.1 installed on AMI EC2 instance and now we are facing space problem. We have few map services around 12 to 15 map services and only 1.5 GB space remaining on C drive.

I spoke my IT admin and contacted support but not get satisfactory answer.

I have stored all data on D drive still wondering that few service definitions files are created on C drive.

The arcgisserver folder size is 2.35 GB still it is showing that all my C drive is approximately full (Please check the screen-shots and please explain me why it’s taking the too much space.)

According to research, I found following options:

Completely move (the arcgisfolder) to D drive
Increase the space of C drive ( I don’t want to restart my machine because its production machine and REST services has the
ip address so If I restart the machine the IP address will change and I need to update REST url with new IP address [Which is not possible : ( ]
Delete all services (& recreate it on D drive, which is not possible.)
My IT admin has checked the all settings of C drive and we are wondering that why it’s taking too much space on Amazon instance (We do not have any cache services)
Going through following links

https://forums.aws.amazon.com/message.jspa?messageID=288313

http://www.rdacorp.com/2011/03/how-to-increase-windows-2008-r2-ebs-volume-size/

https://stackoverflow.com/questions/5898308/how-to-safely-upgrade-an-amazon-ec2-instance-from-t1-micro-to-large

I remember that in previous version i.e. version 10.0 SP3 I did not face this type of problem.

Please check this and let me know suitable option so that we will do this is on Downtime.

Software's : ArcGIS Server 10.1 & ArcGIS Desktop 10.1 , O.S : win 2008 R2

Any help will be highly appreciable …thanks in advance : )

Please check the win_dart information of my machine

Image 1:


texts: Removing Shadows from Aerial Imagery in ArcGIS

texts: I am creating an automatic change detection tool in Python for ArcGIS 10.2 using aerial imagery.  The tool works but it also detects shadow change. I have been able to create a classified raster of just shadows and tried using raster calculator to eliminate the shadows using the Con statement and subtraction without success. Is there a way to remove shadows from an aerial image or make those areas NoDATA, using the tools in ArcGIS and Python.  I read about histogram matching and using a mask to eliminate shadows though I am not sure how to go about doing this.


texts: WGS84 Lon/Lat from pixel x/y in a georeference image

texts:
I'm quite a newbie in GIS.

I need to compute (without using any existing lib/tool) the longitude & latitude (WGS84) corresponding to the a pixel, inside a (georeferenced) map image.

I know:

the x,y coordinates or the pixel,
the width & height (pixels) of the map image,
the bounding box (lon/lat coordinates for the top-left, top-right, bottom-left & bottom-right vertices).

Intuitively I try to use longitude difference between top-left/top-right, in correlation with image width and pixel x to determine the corresponding longitude.

Please be open minded. I don't ask for a ready to use solution, but for documentation to understand the key concept required to do it.


texts: Selecting features within bounding box using QGIS?

texts: What is the best way within QGIS to select all of the features within a rectangular bounding box? Is it necessary to create a separate layer with the box, then to use the intersection tool within vector research -> select by location?


texts: Half disappearing line rectangle when defining base height

texts:
I am having an issue in ArcScene. I have a spline surface representing the topography of an area, and want to create a solid fill beneath this. I can't do this form the spline layer itself or a (rectangle) polygon of the spline area. I have now tried to converted the polygon to a line feature, however, upon defining the base height (from the spline layer), only one length of the rectangle appears.

Why is this?


texts: Custom function for label in QGIS Server GetPrint request

texts:
I'm trying use GetPrint request with QGIS Server and I find feature that you can assign value of function to label (with id) ex.

http://qgis.local/cgi-bin/qgis_mapserv.fcgi?MAP=/home/mario/wms.qgs&VERSION=1.3.0&SERVICE=WMS&SRS=EPSG:4326&REQUEST=GetPrint&LAYERS=wms&STYLES=&FORMAT=png&TEMPLATE=Template1&field1=[%25 now() %25]

From QGIS 2.8 there is possibility to make your own function so I try with built-in custom function func(value1)

http://qgis.local/cgi-bin/qgis_mapserv.fcgi?MAP=/home/mario/wms.qgs&VERSION=1.3.0&SERVICE=WMS&SRS=EPSG:4326&REQUEST=GetPrint&LAYERS=wms&STYLES=&FORMAT=png&TEMPLATE=Template1&field1=[%25 func('a') %25]

Is there way to make it work?


texts: Adjust Color Theme of CartoDB Basemap "Dark Matter"?

texts:
I wanted to know - but haven't found the solution elsewhere - whether it is possible to adjust the color theme of the CartoDB basemap "Dark Matter".

I'd only like to change the colors of the water and country areas as on some computer screens one unfortunately cannot see the differences between water and land. I would like to change the black color of the land to dark grey and the water color to a slightly brighter shade of grey in order to accentuate the contrast a bit more.

Is this possible? Or do I have to stick to the color theme that has been chosen by your cartographers?


texts: GeoServer 2.7.1 fails to start after installing 11-Aug-2015 scripting extension

texts:
This problem does not appear to be version specific for me. I have tried this several times with different versions and different system configurations. As the title states, when I add the script extension to GeoServer it stops working. Any suggestions? I am running this on Debian 8.

Here is an excerpt from data/logs/geoserver.log:

  2015-08-12 11:52:31,066 WARN [support.DisposableBeanAdapter] - Invocation of destroy method failed on bean with name 'geoServerLoader': java.lang.NullPointerException

  2015-08-12 11:52:31,067 ERROR [context.ContextLoader] - Context initialization failed
  org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'wcsURLMapping' defined in URL [jar:file:/var/lib/tomcat8/webapps/geoserver-test/WEB-INF/lib/gs-wcs-2.7.1.jar!/applicationContext.xml]: Cannot resolve reference to bean 'catalog' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'localWorkspaceCatalog' defined in URL [jar:file:/var/lib/tomcat8/webapps/geoserver-test/WEB-INF/lib/gs-main-2.7.1.jar!/applicationContext.xml]: Cannot resolve reference to bean 'advertisedCatalog' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'advertisedCatalog' defined in URL [jar:file:/var/lib/tomcat8/webapps/geoserver-test/WEB-INF/lib/gs-main-2.7.1.jar!/applicationContext.xml]: Cannot resolve reference to bean 'secureCatalog' while setting constructor argument; nested exception is org.springframework.beans.factory.BeanCreationException: Error creating bean with name 'secureCatalog' defined in URL [jar:file:/var/lib/tomcat8/webapps/geoserver-test/WEB-INF/lib/gs-main-2.7.1.jar!/applicationContext.xml]: Instantiation of bean failed; nested exception is org.springframework.beans.BeanInstantiationException: Could not instantiate bean class [org.geoserver.security.SecureCatalogImpl]: Constructor threw exception; nested exception is java.lang.NullPointerException
      at org.springframework.beans.factory.support.BeanDefinitionValueResolver.resolveReference(BeanDefinitionValueResolver.java:329)
   ...


texts: Multi-Band (RGBA) into single band raster

texts:
I have the raster which contains 4 bands (RGBA), which ilustrates the pupulation (traffic) density. I want to convert it to the simple gridded ascii, however I'm getting in the trouble when I try the export it to the different formats. I suppose it comes from multi-band property of this map, therefore:

Is there any way I can convert from RGBA to single value representing the data I need ?

E.g. when I want to convert using RGB2PCT:

  ERROR 1: TIFFOpen:traffic_przyciety2: Permission denied  Traceback
  (most recent call last):  File "C:\PROGRA~2\Quantum GIS
  Lisboa\bin\rgb2pct.py", line 136, in
  tif_ds.GetRasterBand(1).SetRasterColorTable( ct )  AttributeError:
  'NoneType' object has no attribute 'GetRasterBand'

BTW: AttributeError: 'NoneType' error is dropped in other functions as well when I want to do something on this map

Info about the map:

  Band 1 Block=903x1 Type=Byte, ColorInterp=Undefined   Min=0.000
  Max=255.000    Minimum=0.000, Maximum=255.000, Mean=254.667,
  StdDev=9.206 Band 2 Block=903x1 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=124.000    Minimum=0.000, Maximum=124.000, Mean=119.886,
  StdDev=7.547 Band 3 Block=903x1 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=0.000    Minimum=0.000, Maximum=0.000, Mean=0.000,
  StdDev=0.000 Band 4 Block=903x1 Type=Byte, ColorInterp=Undefined
  Min=0.000 Max=255.000    Minimum=0.000, Maximum=255.000, Mean=217.937,
  StdDev=67.026'

Alternatively, do you know if there's any standard interpretation of the geodata maps ?
E.g. if the resulting pixel value may be obtained by specific operation on the bands ?

Thanks for the answers

BR
Krzysztof


texts: PyQgis : Creates table coordinates of selected feature and add it to layout

texts:
I am working on a plugin that allows users to look for a feature in the database, select it, and then open the Layout designer with map, arrow legend and some labels.
I am looking to add a table that contains the selected feature coordinates if someone could help.
I already tried the code given here : Adding QgsLayoutItemTextTable to a print layout
Once the script executed, the Layout designer opened and Qgis crashes and shut down
NB : When I delete the last line, the layout is created. But when I open it manualy. Qgis crashes.
project = QgsProject.instance()
layout = QgsPrintLayout(project)
layout.initializeDefaults()

table = QgsLayoutItemTextTable(layout)
layout.addMultiFrame(table)

# Add columns
cols = [QgsLayoutTableColumn(),QgsLayoutTableColumn(),QgsLayoutTableColumn()]
cols[0].setHeading("heading 0")
cols[1].setHeading("heading 1")
cols[2].setHeading("heading 2")
table.setColumns(cols)

#Add only 1 row
table.setContents([['hello', 'world', '!']])

#Base class for frame items, which form a layout multiframe item.
frame = QgsLayoutFrame(layout, table)
frame.attemptResize(QgsLayoutSize(30, 30), True)
table.addFrame(frame)
# Open the layout designerfor preview
iface.openLayoutDesigner(layout)

After I deleted the #Add columns section from the code below :
the frame is added and the layout designer opened
But the frame is shrinked and the row isn't displayed.


texts: How to select among overlapping linear features by create date?

texts:
I have a table which has overlapping linear data; for example:

A----------C and B---------C.

Of these two records I would only like to return the record with the latest create date; in other words, the record created last.  How do I set up up my Search Cursor to produce the desired result?  Really what I'm attempting to perform is dynamic segmentation depending on the create date.  My end result can be A---B then B---C or if the create date of the first record (which is from A-----C) is newer than the create date for B-----C, I should just return record A----C.


texts: Adding association for location with city or zip code when using PostGIS to add markers using geography?

texts:
On my front end, I am using google maps JS api.  I am adding markers, the google maps event gives me the lat/lng for where I add the marker.

I'm adding that marker into my postgres database:

INSERT INTO markers(created_by, title, description, lat, lng, geography)\
values($1, $2, $3, $4::decimal, $5::decimal, ST_SetSRID(ST_MakePoint($5::decimal, $4::decimal), $6))

Now, I have the marker added into my database and from the front end I can show the user which markers are within a mile or whatever of his current location.

However, I would like to start showing results by zip code or the city associated with the user's current location.

For example, let's say I am adding this latLng to my database: 37.555880, -122.325580.  In my table I'm just going to add these points and the geography.  But, this is in the city of San Mateo, CA, and the zip code is 94402....

How can I add this information to my database given just the lat and long?

If it isn't possible, then I will look into using the reverse geocoding API from Google going forward for future points, so I can get this information from the front end itself and insert it.


texts: Working with orthophotos on an SSD drive vs HDD drive: is the performance gain significant?

texts:
I'm shopping for a server that will allow my coworkers to use QGIS through a terminal server. We have 2 Tb of orthophotos in tif, twf, ecw and shapefiles format. They tell me that actually when they work with orthophotos it can be sluggish. The way they describe it, the loading time of the orthophotos seems to be the culprit. They work with the orthophotos stored on an external HDD through an USB 2.0 connection.
I consider buying extra big SSD drives to store the orthophotos on them to assure the best performance, however an IT consultant tells me he thinks it would be relevant and the performance gain compared to using HDDs for those file's storage would not be significant (we will however put the OS and the main files on SSD drives). Has anyone experienced working with orthoimages on SSD drives, and is the performance gain significant compared to working with orthoimages on a HDD ?


texts:
I am trying to create a grid in QGIS and keep getting a error saying "Horizontal Spacing is too Large for covered area". I'm trying to create a 63.93m x 63.93m square grid.
I am using the projection UTM Zone 10N or EPSG:26910. Both the project CRS and Layer are in this Projection.
I feel my error is somewhere in the projection. Even if I try to simply measure a line the units are very small. I have attached some images for reference. Where am I going wrong?


texts: Create Filter Parameter in Python Script Tool from Unique Set of Values in Field

texts:
I currently have a python script tool that reads in a feature class and (among other processes) subsets the feature class based on a value in one field, 'REGION'. Instead of having all subsets created, I'm now seeking to be able to run the tool for certain regions and am trying to figure out how to have a parameter in the tool give a drop down list of the unique values within the REGION field. I've gotten up to the point where I can select the region (which I don't even need as a parameter since it will never change and always be based off of REGION), but I'm stumped on how to have a parameter of a drop down of a list of unique values in REGION.

Here is a the code I'm trying to use. I've been looking at http://desktop.arcgis.com/en/arcmap/10.5/analyze/creating-tools/customizing-script-tool-behavior.htm but don't quite understand how to have the set of REGION values appear in the python tool GUI.

import arcpy
from arcpy.sa import *

featureClass = arcpy.GetParameterAsText(0)

class ToolValidator(object):
        def _int_(self):
            self.params = arcpy.GetParameterInfo()

        def initializeParameters(self):
            return

        def updateParameters(self):
            rows = arcpy.da.SearchCursor(featureClass, REGION)
            self.params[0].filter.list = sorted(list(set(r.getValue('REGION') for r in rows)))
            del rows
            return

        def updateMessages(self):
            return

def new_feature_class():
# 1. For each feature class selected create a new feature class for every unique region within it.
    for j in featureClasses.split(';'):
        def unique_values(table, field):
            with arcpy.da.SearchCursor(table, field) as cursor:
                return sorted({row[0] for row in cursor})

        uniqueRegionList = unique_values(j, 'REGION')
        print (uniqueRegionList)

        for i in uniqueRegionList:
            in_features = j
            out_feature_class = str(i) + str(timePeriod) + in_features[-5:]
            where_clause = " \"REGION\" = '%s'" % (i)

            arcpy.Select_analysis(in_features, out_feature_class, where_clause)
            print(out_feature_class)
new_feature_class()

Edited to add Validation Code:

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

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

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

  def updateParameters(self):
    rows = arcpy.da.SearchCursor(featureClass, REGION)
    self.params[1].filter.list = sorted(list(set(r.getValue('REGION') for r in rows)))
    del rows
    return

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


texts: How to create polygons in ArcMap 10.0 using lat/long data from SQL Server 2008 database?

texts:
Forgive me if this is very basic, as I'm new to GIS work (though I've been a programmer for 10 years) and so I'm not really sure how to approach this, but here goes...

Can I take a table with xmin/xmin ymin/ymax coordinates and use it to create a bunch of polygons in ArcMap 10.0? I have a connection set up to SQL Server 2008 wherein there is a table that has several pieces of information regarding various maps that our team keeps track of.

Here's what I'm working with:

Said table contains the lat/long data in the form of xmin/xmax ymin/ymax, in addition to several other attributes - gmap_id, author, year, scale, etc. My boss asked me if I could create a connection to our SQL Server, use that information to create polygons on a map of the US, and set up a little attribute window with the other information (author, gmap_id, scale, etc) for each polygon.

I've got the first part covered, I think, thanks to this question. I set up the Data Query and pulled the columns that I initially wanted from SQLServer. That gave me a table like so:

I had a screen cap but the image won't post due to my rep < 10, so the table looks like:

gmapId, scale, usgsIdF, xmin, xmax, ymin, ymax
18015, 250000, 13139, -121.891, -123.071, 43, 45.5
18016, 125000, 13140, -117, -118.936, 34.25, 35.624
...

And I now have the data listed as a Source under the Table of Contents window.

How can I now translate that information into a bunch of squares on the map? And then can I use that to set up an attribute table?

I have a lot of python experience, so if that is the way to go I can do it. I am also good with SQL Server 2008, though I have not used it's spatial capabilities yet. So if there is a solution to be had involving TSQL I can also go that route.

I'm just scratching around trying to figure a good way to start.


texts: C# / .NET library for converting GeoJson to KML

texts:
I'm pulling in data from ArcGIS using ArcGIS.PCL which provides the geometries as GeoJson. I'd like to convert this GeoJson to KML. Does anyone know of any .NET libraries that provides this functionality? I've had a good search but nothing has popped up on my radar.

I guess a couple of libraries may be required if there's no catch all for this problem. One that can convert GeoJson to WKT then another that can convert WKT to KML.


texts: Calculating relative share of an area of circle that has slope >15% in Germany

texts:
I need to calculate the share of buildable land for the 15 major german cities within a 50 kilometer radius of the city center.  Part 1 of the definition buildable is that the area must have a slope <15% (Part 2 would be being no water surface or wetland, but this might be covered elsewhere).
As a student I am completely new to GIS Software and have already tried to mess around with a DEM and the Raster Calculator in ArcMap and QGIS after watching/reading some tutorials but encountered several problems. As my guess/hope is that this should be quite an easy task for an experienced user, I thought I would rather just ask if someone would be very nice and could give me some instructions how to do that first before talking about specific problems in detail(file formats, z-values, ArcGiS Online map..) which may just be there due to a wrong approach in the first place.  Software available would be ArcMap 10.3.1 and QGIS 2.12.3 and the DEM Model used should be 90m resolution or better.  Please recommend a publicly available DEM to use, as I don't have a working one as of now.  In case projections, reference systems, stuff I know not deeply about yet, might lead to different distortions, the results do not need to be extremely precise.
Perhaps there is a tutorial where this task is somehow covered?


texts: Set bar colors on a Google Earth Engine chart

texts:
I am trying to define a set of colours by month on a ui.Chart.image.series, however, after adding the .style line, no changes are made. I would like to have a different bar color for each month, the palette should repeats over different years.

My code so far looks like this:

var dataset = ee.ImageCollection('TRMM/3B43V7')
                  .filter(ee.Filter.date('1998-01-01', '2019-06-30'));
var precipitation = dataset.select('precipitation');
var precipitationVis = {
  palette: ['blue', 'purple', 'cyan', 'green', 'yellow', 'red'],
};
Map.setCenter(6.746, 46.529, 3);
Map.addLayer(precipitation, precipitationVis, 'Precipitation');

// Define a region of interest as a buffer around a point.
var geom = ee.Geometry.Point(-74.08384, 9.42503).buffer(500);

// Create and print the chart.
var mean_prec_chart = ui.Chart.image.series(precipitation, geom, ee.Reducer.mean(), 12);
mean_prec_chart.setChartType('ColumnChart');
mean_prec_chart.style({color:['blue','red','green','brown']});
print(mean_prec_chart);


texts:
I am using the MergeShapes plugin with QGIS 1.8.0 on a Windows 7 computer.  The issue I am having is when I perform a merge from a directory based on one geometry type, and then switch the geometry type to process a second time (ie: process all polygons in c:\Temp and then process lines in the same directory).  This is a limitation when you need to preform multiple merges based on geometry.  Any help would be appreciated.  Below is the error from MergeShapes.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/ryangarnett/.qgis//python/plugins\mergeshapes\mergeshapesdialog.py", line 141, in accept
    self.progressFiles.setRange( 0, self.inputFiles.count() )
AttributeError: 'NoneType' object has no attribute 'count'

Python version:
2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

QGIS version:
1.8.0-Lisboa Lisboa, 6416f38

Python path: ['C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/ryangarnett/.qgis//python', 'C:/Users/ryangarnett/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\PROGRA~2\Quantum GIS Lisboa\bin\python27.zip', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\DLLs', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\plat-win', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\bin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\Users\ryangarnett\.qgis\python\plugins\qgis2kml', 'C:\Users\ryangarnett\.qgis\python\plugins\qgis2kml', '/usr/share/qgis/python', '~/.qgis/python', '.', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\python\plugins\fTools\tools', '~/.qgis/python', '.']


texts: How does Ellipsoidal Height and Geoidal Height impact water flow?

texts:
I see this as a common theme

Ellipsoidal height doesn't dictate water flow.

Take a shot at explaining this. I think I understand this but it sounds like it would be of use to have a question about it.

Who has the best answer as to why it doesn't?

Rather than just repeating it and linking to ESRI and other docs about it.


texts: Adding parent folder name to feature name in ArcGIS ModelBuilder?

texts:
I have created many shapefiles with the same name and every shapefile is in its own folder with a specific name.

For example:

folder1\polygon.shp
folder2\polygon.shp
folder3\polygon.shp

And I want the features to be like: folder1_polygon.shp, folder2_polygon.shp, etc...

How to add a parent folder name to a feature name using the ModelBuilder?


texts: Accessing CRS info with pyqgis in a standalone application

texts:
I'm currently having trouble accessing and defining CRS information in a standalone application. I've read the API documentation and have tried it in the python console in QGIS 2.2.0. The python console works fine. I'm able to get the information I need. But the layers in my stand alone application do not have have valid CRS's. What is the reason for this?

Code:

import qgis
from qgis.core import *
from qgis.utils import *
from qgis.gui import *

import pdb

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

print "Testing Json crs:"
layer = QgsVectorLayer("E:/GIS_Data/GeoJsonRendering_Test/test.json" ,  'musym', "ogr")
QgsMapLayerRegistry.instance().addMapLayer(layer)
print layer.isValid()
print layer.crs().isValid()
print layer.crs().authid()

print "Testing QgsCoordinateReferenceSystem:"

t = QgsCoordinateTransform()
projSrc = QgsCoordinateReferenceSystem(3857, QgsCoordinateReferenceSystem.EpsgCrsId)
projDest = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId)

print projSrc.isValid()
print projDest.isValid()

print "Testing memory layer CRS:"
vTemp = QgsVectorLayer("Point?crs=EPSG:3857", "temporary_points", "memory")
print vTemp.isValid()
print vTemp.crs().isValid()

Output:

Testing Json crs:
True
False

Testing QgsCoordinateReferenceSystem:
False
False
Testing memory layer CRS:
True
False

So I have two valid layers without valid projections.

Here is the output for that code in the python console in QGIS:

Testing Json crs:
True
True
EPSG:3857
Testing QgsCoordinateReferenceSystem:
True
True
Testing memory layer CRS:
True
True

So my code only works in the QGIS python console. I'm assuming that I am not defining settings that need to be defined, but I"m not quite sure what I'm missing.


texts: ArcMap 10.4.1, How to show linked pictures in layout data frame for a map book (using Basic License)

texts:
I'm trying to create a map book for a tree inventory in ArcMap 10.4.1. I have a shapefile (no problem changing to a gdb) with 350 points, and the goal is to get a layout consisting of 3 data frames: two for maps showing each tree location with different layers underneath and the third one to show the picture of the tree itself. The first 2 are no problem, the issue is with the 3rd one.

The shp also has several fields that will show info in the layout via "Data driven page attribute", and what i'm having problems to achieve is a workaround to do this exact same thing but for images. Something in the attribute table that links the feature to its correspondent image, and for every page of the map book you'll end up with 2 data frames shoing the point location, and another one showing the picture linked to the point.

I know you can use attachments to do this, adding pictures and then setting data driven pages so each picture shows in the layout with the related point, but tied to work with a Basic License and this features are restricted to Standard and Advanced only.

Also i have tried working a gdb instead of a shp, creating a new field that lets attach raster images, but i have not found the way to show them in the data frame.

Is there any way to achieve this without having to use the attachment tools?


texts: Combining DEM file with building (height) vector files using ArcGIS Desktop?

texts:
I am using ArcMap 10.2. I need a DSM file to use for my thesis to run visibility analysis on it, but my case study (Tehran) do not have LiDAR data.

So I am thinking to combine a DEM file with building vector files which includes height.

Does anyone have a solution?


texts: Merging shapefiles with different geometries (point and polygon) to get single layer with multiple geometries?

texts: I want to merge two shape files in which one is of point geometry type and other is polygon type in order to get a single layer containing both points and polygon in a single map.


texts: Clipping multiple rasters by matching polygons from SPDF in R

texts:
I have about 50.000 small glacier rasters (rectangular shape) that need to be clipped to the glacier's extent. I have the glacier polygons stored in a large SpatialPolygonsDataFrame.

I want to use parallel computing or else it would take forever, I think. For previous tasks I've successfully used the mcmapplyfunction, but I am open for other approaches.

My (admittedly rudimentary) code so far is:

filenames <- list.files("/.../RGI60-13_reproj/", pattern="*.tif", full.names=F)
filelocations <- list.files("/.../RGI60-13_reproj/", pattern="*.tif", full.names=T)
glaciers <- readOGR("/.../13_rgi60_CentralAsia.shp",verbose=TRUE)

fun_clip <- function(filelocations, filenames, glaciers){
  r <- raster(filelocations)
  r <- crop(r,glaciers) # here I need to clarify the corresponding shp in the SPDF
  writeRaster(r, paste0("/.../RGI60-13_crop/",filenames))
}

mcmapply(fun_proj, filelocations, filenames, mc.cores = 50)

How can I give the crop-function the right iterative arguments? filelocationis of the same length as glaciers, so in a for-loopI would use something like
r <- crop(r,glaciers[i]), but how do I pass the iteration in my kind of function?
What would be the way to introduce the i, so to speak?


texts: Conversion of design file to shapefile using FME

texts:
I want to convert a MicroStation design file to shapefile using FME. I have two type of feature in a level 7. One has style 1 and another one style 3. I want them in different shapefiles.

How can I get FME to differentiate between them?


texts: Geoprocessing service fails if output parameter is empty?

texts:
I have made a network isolation model that lets the user pick a point on the network to simulate a pipe 'failure'/section that needs to be isolated and finds all valves that need to be turned and other affected features (hydrants, buildings, ect.). The model will output feature layers of those features found in the trace from the 4 feature classes, but only if a feature is found from the trace operations. This works fine in ArcMap when ran from the catalog, but not when ran from the server.

The problem is that not all outputs will contain a record for each FC. For instance, a trace may only return a pipe section and 2-3 valves, with no building connection points or hydrants returned. This works fine in ArcMap because i use a try/except to assign the results to the output parameters of the script, so if one of the FC has no selected records, the output is not assigned with anything. When ran on the server, this does not seem to work. When i run a process with no hydrants resulting from the trace, the model fails with Error Code 20010 - Invalid return value.

Does the server not adhere to the try/except function or is there any way to return an empty feature layer if no results are returned for that FC?


texts: Leaflet: Markers hidden when panning map

texts:
I've implemented a Leaflet map with a couple of different custom markers using DivIcons.

It works fine in desktop browsers, but when I pan the map on a mobile device the markers that aren't currently visible in the viewport seems to be hidden until I stop panning the map.

Is there any option for disabling this behavior? I would ideally like to keep the markers rendered all the time.

I should mention that I'm using the MarkerCluster plugin as well.


texts:
I am using Anaconda 2.1.0 and Basemap installation is successful

conda install basemap

However, I keep getting the following error when I try to project lat,lon to a planar:

Assertion failed: (0), function query, file AbstractSTRtree.cpp, line 286.

My code looks like this:

import shapely.geometry as sg
from pymongo import MongoClient
from mpl_toolkits.basemap import Basemap

p = Basemap(projection='sinu',lon_0=0,resolution='c') #fails here

projected_xy = p([24.4,45.1])

The code fails even before attempting to project.


texts: How to create a "Datum-less" projection in QGIS?

texts:
I want to use QGIS to do some georeferencing and image registration. On ArcMap10 I can define an "Unknown Projection", basically a default metric grid, and use it as a grid for georeferncing.

In QGIS however, when I tried to georeference an image with GDAL, it gave me the following error:

  Failed to compute GCP transform: Transform is not solvable

I've looked around and figured that the problem lies with the projections.

I've imported a 4 line CSV with 4 points demarcating the image extent in meters (BL 0,0, UR 450,450), but the layer properties tell me it's in WGS 84. because I do not know the actual coordinates of the image, I have to use just a basic internally coherent 5*5 metric grid, that I can georeference with corresponding marks in the image. This can only work with a datumless metric projection, which I do not know how to define (and save it as a custom defined projection)

I'm new to QGIS, so the problem might be remarkably easy to solve, but after a few hours of searching, I still haven't found it.

appreciate any help,
H.


texts: Sub expression not considered when using aggregate functions (QGIS 3.8)

texts:
I am attempting to use aggregate() or relation_aggregate() in order to populate a label on an atlas powered print composer using QGIS 3.8, but am having issues with both.

My goal is to show the number of children that satisfy a certain query on the an atlas that is powered by the parent.

I first created a relation in the project properties, saved, closed, and reopened project.

I added a label to my print composer, and have tried the following syntax:

Using aggregate

[%aggregate('child_layer', 'count', "child_layer_field" is 'OUI', attribute($currentfeature, 'foreign_key') = attribute(@parent, 'primary_key'))%]

Using relation_aggregate:

[% relation_aggregate('relation_key', 'count', "child_layer_field" is 'OUI')%]

In both cases, the label displays a number, but it is simply the total number of children related to the parent (not taking into consideration the expression "child_layer_field" is 'OUI').

Am I misunderstanding how this works or am I doing something wrong?
Should this work with 2.18 as well?


texts: How to customize date format in ArcGIS Operation Dashboard's list widget?

texts:
I'm trying to add to my dashboard a list widget to show a list of recently updated records. I'm able to easily add fields by putting the field name between braces (ex {LastEditDate}), like in arcGIS Javascript API's PopupTemplate.

What I want to do now is to customize the date format. In this case only the time is relevant, but in other cases I might want only the date.
I tried to use a syntax similar to javascript API ({LastEditDate:DateFormat(selector: 'date', fullYear: true)}), without success.
Is there a way to customize a field's date format in the list widget?


texts: return only one value per point from ST_Value() - currently returning 2 values for some points

texts:
I'm testing a point sample query with a subset of 5 points on a single band elevation raster. Using my usual approach (which I've previously not noticed this behavior), I get six records returned:

SELECT t1.pt_id, CEILING(ST_Value(rast, geom)) as elev_m
FROM (
    SELECT pt_id, geom
    FROM test.pt
    LIMIT 5
) t1
LEFT JOIN test.raster t2 ON ST_Intersects(t1.geom, rast);

Output:

It's clear which record is wrong. I understand that the offending point intersects with two pixels (of possibly two different tiles) and therefore, is returning two values.

I temporarily fixed this particular problem, by eliminating the duplicate pt_id and the value of 0 with:

SELECT t1.pt_id, CEILING(t2.elev_m)
FROM (
    SELECT pt_id, geom
    FROM test.pt
    LIMIT 5
) t1
JOIN LATERAL (
    SELECT t1.pt_id, MAX(ST_Value(rast, t1.geom)) as elev_m
    FROM test.raster
    WHERE ST_Intersects(t1.geom, rast)
    GROUP BY t1.pt_id
) t2 ON TRUE;

But this solution will not be reusable necessarily, because other rasters may have outliers that are high, not low like in this example. Also, I'm not concerned about a few outliers in my dataset. I'm more concerned with returning only raster value per point.

So, how can I return one value per point, preferably the value of the pixel that it intersects 'most' with? I need a general solution that can be re-used consistently across different datasets. Also, efficiency is key, as this will be scaled up to millions of points once I get this worked out.

I found this question, but its answers remain unaccepted and did not help me for this particular problem.

Again, my question is: How can I create a reusable query that will ensure that only one raster value is returned per point?


texts: PostGIS; won't update geometry when importing CSV file with psql

texts:
I installed PostgreSQL 12 and PostGIS on a new server.
I want to use this environment as source for different projects and services.
When I try to import CSV files in the database using psql, I get an error when trying to execute the following:
UPDATE bet_srk_2020 SET geom =
ST_SetSRID(ST_MakePoint(AISTrajectPositionLongitude,
AISTrajectPositionLatitude) ,4326);

Result:
FATAL:  failed to load summary "/usr/lib/postgresql/12/lib/bitcode/postgis-3.index.bc": Invalid summary version 7, 1, 2, 3 or 4 expected
server closed the connection unexpectedly
    This probably means the server terminated abnormally
    before or while processing the request.
The connection to the server was lost. Attempting reset: Succeeded.

I tried to find information about this error but can't really find solutions.
I am able to add data to the databases if I use the QGIS Toolbox feature:
Export to PostgreSQL
Is this simply a syntax error in my request to the server? Or is there a problem with the installation of PostGIS?
EDIT: 2020 10 30:
Based on CL. comments I investigated the problem in my installation:
I use a Ubuntu setup with the apt.postgresql.org repository.
For some reason Postgresql-12 was installed from the postgresql.org repository and postGIS 3.0 was installed from the standard Ubuntu repository. After installing postGIS from the same repo as postgreSQL the problem was solved and server operation is back to normal.


texts: Can WMS be used to get a local copy of raster data?

texts:
I have tried reading up on WMS on wiki and other sources but haven't managed to get my head around a couple of points.

1) Using WMS is it possible to get a local copy of the raster data (say a tiff image tile). The reason, I ask is, I need to use his raster data in a secondary process that needs a local copy of the file.

2) I am primarily interested in allowing users the option of downloading orthoimagery and elevation from the USGS repository. I've tried using the Seamless viewer but I am thinking more along the lines of user specifying a coordinate (or four corner coords) and getting the orthoimage and elevation data corresponding to that. I thought WMS (Geoserver being my first option) is the way to go for this, but not certain. So wanted to get your advice.

Thanks in advance.

Sincerely,
ST.


texts: Choosing correct Coordinate Reference System when doing image georeferencing in QGIS?

texts:
I am trying to follow the tutorial in http://www.qgistutorials.com/en/docs/georeferencing_basics.html to do the image georeferencing. But every time I got strange results.

My image is self created picture, with several known points (I obtain the coordinate of these points from handheld GPS receiver). I'd like to match the image to real world coordinate.

Following the tutorial, I got a geotiff file loaded into QGIS. The problem is that its position is totally wrong. Wrong position is around (0,0) in (lat, lng).

When I open the CSR property of the geotiff file, the coordinate reference system shows Selected CRS, USER100000 (as I'm expecting EPSG3857).

Could anybody help with this?


texts: netcdf query ocean temperature at different depths

texts:
I am looking for the best way to query a netCDF file - I would like to eventually expose this to an API to be consumed by a web service.

I would use SQL which I have experience in but I think its probably going to be resource heavy. I have looked into geoserver but I dont see any documentation on query this type of information only serving WMS tiles.

I have a grid of temperatures at different depths, I would like to for a given lat / lng find the nearest source (looking at just ground distance not vertical) of water at a certain temperature.  For example lat 30.083 / lng -14.5 find the nearest location of water at 20c.

My data has an array of data points for depths from 5 to 5000m the query will need to return the lat / lng and depth.


texts: What will GIS be like in 5 years (and in 50)?

texts: How will GIS technology be used 5 years from now ... and in 50 years?


texts:
I have added bar charts to my map and this works fine. But I have no idea how to add values to the bar chart. For example one bar shows the value of 2000 and another shows the value of 2100. Is it possible to show this values in the bar chart? I tried labeling but this is not what I want.

Any ideas? A Plugin or something else?

qgis 2.01 and 1.8

Thanks.


texts: How to get selected layers in QGIS3 from python

texts:
How can I get the selected layers in QGIS 3 with pyqgis?

I tried

selectedLayers = iface.legendInterface().selectedLayers()

as suggested in How to get programmatically selected layers from the QGIS legend?, but it seems that there is an API breaking change

  QgsLegendInterface was removed. It was replaced by layer tree API
  (QgsLayerTreeNode class and others). Methods that deal with custom
  actions in main window's layer tree context menu were moved to
  QgisInterface:

I then tried:

selectedLayers = QgsProject.instance().layerTreeRoot().selectedLayers()

but I get an error about:

AttributeError: 'QgsLayerTree' object has no attribute 'selectedLayers'

What's the correct way to get the selected layers in QGIS3?


texts: Choosing opengis or geotools imports dealing with geotools library?

texts:
In some cases, I have to use some classes of GeoTools library.

Most of class names are duplicated in org.geotools and org.opengis package. for example, StyleVisitor, Symbolizer, PointSymbolizer, etc.

I want to know which one of them I should use in my project in usual coding.

import#1:

import org.geotools.styling.Graphic;
import org.geotools.styling.LineSymbolizer;
import org.geotools.styling.Mark;
import org.geotools.styling.PointSymbolizer;

import#2:

import org.opengis.style.Graphic;
import org.opengis.style.LineSymbolizer;
import org.opengis.style.Mark;
import org.opengis.style.PointSymbolizer;

My code:

public class MyStylePainter extends StylePainter {

    private MyStylePainter(int width, int height) {
        super(width, height);

    }

    public void setWidthAndHeight(PointSymbolizer pointS) {
        double d = 15;
        double size = size(pointS.getGraphic());
        double realSize = d * size;
        int width = (int) Math.round(realSize);
        this.setWidth(width);
        this.setHeight(width);
    }

    protected static double size(Graphic graphic) {
      double size = 1;
      if (graphic.getSize() != null) {
          size = Filters.asDouble(graphic.getSize());
      }
      if (size <= 0) {
          size = .2d;
      }
      return size;
    }
}

Do I use import#1 or import#2 ?


texts: Merging multi-polygons gives me error in shapely

texts:
I am trying to merge a bunch of either polygons or multipolygons by the merge_multi_polygons function in shapely, I am using my below code inside Spark

        elems=row['_c0']
        area_multipolygons=[]
        for y in elems:
            json_polygon=[]
            searchObj = re.finditer( r'[(]{1}([^()]*)[)]{1}', y)
            for i in searchObj:
                if i:
                    st_polygon= i.group(1)
                    list_polygon=[float(i) for i in st_polygon.split()]
                    json_loop=list(zip(list_polygon[::2], list_polygon[1::2]))
                    json_polygon.append(json_loop)

            forest_area_multipolygon=[json_polygon]
            area_multipolygons.append(forest_area_multipolygon)

        merged_parttion_polygon=merge_multi_polygons(*area_multipolygons)

I get below error on some of my data

 File "C:\Users\Mahdi\AppData\Local\Programs\Python\Python37\lib\site-packages\pyspark\util.py", line 55, in wrapper
    return f(*args, **kwargs)
  File "main.py", line 147, in <lambda>
    udf_calculate_intersection=udf(lambda x: calculate_intersection(x), FloatType())
  File "main.py", line 138, in calculate_intersection
    merged_parttion_polygon=merge_multi_polygons(*area_multipolygons)
  File "C:\Users\Mahdi\eclipse-workspace\EliizaTest\mahdi\polygon_utils.py", line 75, in merge_multi_polygons
    res = [to_polygon(union)]
  File "C:\Users\Mahdi\eclipse-workspace\EliizaTest\mahdi\polygon_utils.py", line 47, in to_polygon
    e = coords(polygon_shape.exterior)
AttributeError: 'GeometryCollection' object has no attribute 'exterior'


texts: Fishnetting many polygons using each feature extent in ArcMap

texts: I am attempting to place a 10X10 grid over 400 30mX30m polygon squares. I have been successful using the Create Fishnet tool in ArcGIS to create a polyline grid feature within a specified polygon but only if the polygon is in its own shapefile. I used the Split Layer by Attribute tool that I found on this site to create 400 separate shapefiles but have not found a way to create a fishnet over each polygon all at once. I am dreading performing this function 400 times.


texts: differences in distance between shapely and geopy

texts:
I am getting totally different distances from shapely and geopy. On one hand - it makes total sense, shapely measures on a plane and geopy measures on a spheroid. But the differences are so dramatic, I am wondering if something is truly off here.

Using:

Shapely 1.3.2
geopy 1.9.1
Fiona 1.l.5

my prj (for shapely/fiona) IS Web-Merc EPSG: 3857 (w/o flattening):

     >>print input_crs
     >>{u'a': 6378137, u'lon_0': 0, u'k': 1, u'y_0': 0, u'b': 6378137, u'proj': u'merc', u'x_0': 0, u'units': u'm', u'no_defs': True}
     >>print to_string(input_crs)
     >>+a=6378137 +b=6378137 +k=1 +lon_0 +no_defs +proj=merc +units=m +x_0 +y_0

Test Pts (*Arizona/N.Mex general area):

city: POINT (-12352678.6153664 3384049.040913342)
capital: POINT (-12318614.85118365 3417715.06365641)

My code, using shapely for .buffer and Vincenty for the ellipsoid measure (geopy)

def evaluate_point(city_pt, capital_pts):

    #tests whether the capitals are within a certain distance of a
    #city point

    # set the ellipsoid for measurement
    distance.VincentyDistance.ELLIPSOID = 'WGS-84'
    d = distance.distance

    # create a wkt obj of city pt
    city_wkt = dumps(city_pt)

    dist = 300000

    # create a shapely obj of all the capitals
    for pt in capital_pts:
        # create a shapely obj of all the captial_pts
        capital_pt = shape(pt['geometry'])
         # create a wkt obj of the admin1_pt
        capital_to_measure = dumps(capital_pt)

        ## BUFFER USING SHAPELY - DIST = 300000 meters due to prj
        capital_buffer = capital_pt.buffer(dist)
        # city_pt is a shapely obj
        if city_pt.within(capital_buffer):
            distance_between_pts = d(capital_to_measure, city_wkt).meters
            print"within the buffer" + " " + str(distance_between_pts)
            print "city: " + city_wkt
            print "capital: " + capital_to_measure

        else:
            print "outside the buffer"

FOR this data I get this:

city: POINT (-12352678.6153664 3384049.040913342)
capital: POINT (-12318614.85118365 3417715.06365641)

   within the buffer # shapely buffer set to (300000 meters)
   14815370.7894 # Vincenty ellipsoid measured distance reported in meters


texts: Configurable Map Viewer (CMV) - Navigation Toolbar Widget - Zoom box issue in IE10

texts: In playing around with adding the Navigation Toolbar Widget to the CMV I noticed the zoom box doesn't show when dragging a box to zoom in/out in IE10 (the red box).  It works as expected in Firefox, Chrome and Safari.  I don't see an obvious issue.  I'm curious if others have seen/resolved this or have any thoughts on where I may look to resolve the issue?