texts:
I'm currently working on a project where I have to use the map Corine land cover.
https://land.copernicus.eu/pan-european/corine-land-cover/clc2018
(You can download the map in the .tif format)
I'm not familiar with gdal and rasters and I'm using Python to extract the data out of the raster. Can you provide me the most efficient way to extract the data?

Gdal info:

Coordinate System is:
PROJCRS["ETRS_1989_LAEA",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["Lambert Azimuthal Equal Area",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",52,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",10,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",4321000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",3210000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",3035]]
Data axis to CRS axis mapping: 1,2
Origin = (900000.000000000000000,5500000.000000000000000)
Pixel Size = (100.000000000000000,-100.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  900000.000, 5500000.000) ( 56d30'18.51"W, 56d29' 4.75"N)
Lower Left  (  900000.000,  900000.000) ( 23d49'33.58"W, 24d17' 3.04"N)
Upper Right ( 7400000.000, 5500000.000) ( 72d54'22.09"E, 58d57' 9.90"N)
Lower Right ( 7400000.000,  900000.000) ( 40d39'45.75"E, 25d32'40.96"N)
Center      ( 4150000.000, 3200000.000) (  7d30'57.52"E, 51d53' 2.21"N)
Band 1 Block=65000x1 Type=Int16, ColorInterp=Gray
  Min=111.000 Max=999.000   Computed Min/Max=111.000,999.000
  Minimum=111.000, Maximum=999.000, Mean=326.518, StdDev=118.029
  NoData Value=-32768
  Metadata:
    DESCRIPTION=clc18
    RepresentationType=THEMATIC
    STATISTICS_MAXIMUM=999
    STATISTICS_MEAN=326.51842078382
    STATISTICS_MINIMUM=111
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=118.02878635921
    STATISTICS_VALID_PERCENT=24.58

Python Code

import gdal
import numpy
from affine import Affine

lons=[15.174866]
lats=[43.169129]

fn="C:/path-to-the-map/map.tif"

ds=gdal.Open(fn)

transform=ds.GetGeoTransform()
xOrigin=transform[0]
yOrigin=transform[3]
pixelWidth=transform[1]
pixelHeight=transform[5]

aff=Affine.from_gdal(xOrigin,pixelWidth,0.0,yOrigin,0.0,pixelHeight)

x_coords,y_coords=aff*(numpy.array(lons),numpy.array(lats))

band=ds.GetRasterBand(1).ReadAsArray()

x=int(x_coords[0]/pixelWidth)
y=int(y_coords[0]/pixelHeight)

value=band[x][y]
print(value)

I get some value from the raster but the value is not correct.
My guess is that I am not converting the coordinates in the right way. I need to convert the coordinates 43.169129 lat, 15.174866 lon to coordinates used in the map to extract data at that exact spot.


texts: Calculate the class limits from the kernel density depending on the percentage of the source dataset

texts:
My goal is to make the class boundary of a kernel density dependent on a percentage of the input points.
Specifically I calculate the kernel density from 3000 points in a 50m x 50 m grid with a search radius of 16000 m in an area of about 357386 km². Now I want to divide the density result raster into two classes.  One class should contain 75% of the original points and the other class with the highest density values should contain 25% of the original points. So I want to set the threshold value for a density to a percentage of the original data.
Unfortunately, I have not yet found a way to make the class boundary dependent on the percentage of the original data points. Is this even possible?
I would be very happy about support. I work with ArcGIS 10.6.


texts: Openlayers 3 circles from geojson not working with layer tile source different from TileWMS

texts:
I am working with openlayers 3 with a layer from mapbox or openstreetmap and I need to draw some circles on it.

I am able to print the circles with a view with a projection EPSG:4326, but then I have no map. Changing the projection with a transform the map is displayed, but the points are all together. The only way I am made it work is using TileWMS as a source, but I am not able to use it in production environment

Here is my code:

version 1: working with TileWMS

  var source = new ol.source.GeoJSON({
    url: 'geojson url'
  });
  var pointsLayer = new ol.layer.Vector({
        source: source,
        style: new ol.style.Style({
          image: new ol.style.Circle({
            radius: 15,
            fill: new ol.style.Fill({color: 'rgba(170,33,33, 0.8)'}),
            stroke: new ol.style.Stroke({color: 'rgba(170,33,33, 0.3)', width: 8})
          })
        })
      });

  var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
              title: "Global Imagery",
              source: new ol.source.TileWMS({
                url: 'http://maps.opengeo.org/geowebcache/service/wms',
                params: {LAYERS: 'bluemarble', VERSION: '1.1.1'}
              })
        }),
        pointsSource
    ],
    view: new ol.View({
      projection: 'EPSG:4326',
      center: [-82.3, -10.65],
      zoom: 3
    })
});

This is the result

Using mapbox or osm, failing:

  var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
                    source: new ol.source.XYZ({
                        url: 'http://api.tiles.mapbox.com/v4/XXXXX.kcp9cpdn/{z}/{x}/{y}.png?access_token=token'
                    })
                }),
        /*
          new ol.layer.Tile({
          source: new ol.source.OSM()
        }),*/
        pointsSource
    ],
    view: new ol.View({
      projection: 'EPSG:4326',
      center: [-82.3, -10.65],
      zoom: 3
    })
});

This is the result

And finally, changing the view the map is displayed but the circles

  var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
                    source: new ol.source.XYZ({
                        url: 'http://api.tiles.mapbox.com/v4/XXXXX.kcp9cpdn/{z}/{x}/{y}.png?access_token=token'
                    })
                }),
        /*
          new ol.layer.Tile({
          source: new ol.source.OSM()
        }),*/
        pointsSource
    ],
    view: new ol.View({
      center: ol.proj.transform([-82.3, -10.65], 'EPSG:4326', 'EPSG:3857'),
      zoom: 3
    })
});

And the result

Is there a way to make this work?


texts: How can I create a new function inside a existent class in OpenLayers?

texts:
I want to add some functions inside a existent class in OpenLayers. For example, I want to add function A() and B() and a variable C inside OpenLayers.Layer.Vector.

var X = function() {
    A() = funtion { /* something */ },

    B() = funtion { /* something */ },

    C = 'something'

}

So when I create a new layer inside a map object I could call those function by:

map.layer[1].A();

or could change:

map.layer[1].C = 'something else';

I tried many ways but didn't worked well. And most of the time, when I tried to change the C var inside one layer it changed the other too.


texts: Extracting data from OSM with ogr2ogr --- blank output

texts:
I have an issue while extracting data from osm via ogr2ogr.

What I tried:

When I run from command line:

ogr2ogr -f "ESRI Shapefile" example.shp source.pbf -lco ENCODING=UTF-8 -progress -sql "SELECT name, tourism from points where tourism is not null"

I am getting shapefile with all tourism tags, but I wanted to select only specific tags to avoid doing selection afterwards.

So I run (I added zoo and museum to OSMCONFIG.ini):

ogr2ogr -f "ESRI Shapefile" example.shp source.pbf -lco ENCODING=UTF-8 -progress -sql "SELECT name, tourism from points where tourism in (zoo, museum)"

and my output is blank, In the shapefile I will get name and tourism collumns, but they are empty.

Question:

How one can select tags that one only need?


texts: Scale 16 bit to 8 bit within range using GDAL

texts:
I am looking to scale a 16bit unsigned tiff image to an 8bit image using GDAL (OSGeo4W). However I want to constrain the pixel values so that they stay within the range of the 8bit imagery. What I mean is, I want to ensure the visual proportion in the 16bit image is represented in the 8bit image (0-255). Using gdal_translate to convert from 16bit to 8bit cuts the pixel values and does not produce a RGB image. I am unsure of what function and setting to call from within gdal_translate, or if there is a better option.

GDALINFO

C:\>gdalinfo C:\Projects\Stormwater\ForPCI\images\1537TORO_6cm_04761_RGBI.tif
Driver: GTiff/GeoTIFF
Files: C:\Projects\Stormwater\ForPCI\images\1537TORO_6cm_04761_RGBI.tif
Size is 17310, 11310
Coordinate System is `'
Metadata:
  TIFFTAG_DATETIME=2015:05:03 12:27:13
  TIFFTAG_IMAGEDESCRIPTION=UltraCam-Lvl03
--------------
CAM_ID: UC-SXp-1-50215465 [3]
IMG_NO: 4761
RECORD_GUID: 8fc3d7a9-d2e5-40e5-babf-939f1f803dcf
IMG_GUID: C850D0A4-39CA-4BC8-ACD7-81A363D78C6E
FILE_GUID: 86C4F8E4-A2FD-43FB-B3B1-4C4ED8DF0E42
LICENSE_ID: 820342819
SOFTWARE: UltraCam Aerial Radiometry Core 12.6.1408.2501
PIXEL_SIZE_WIDTH: 6 [micron]
PIXEL_SIZE_HEIGHT: 6 [micron]
APERTURE: F_8
EXPOSURE_TIME: 0.002000000000 [s]
HIGH_ISO_MODE_CAPTURE: off
HIGH_ISO_MODE_PROCESSING: off

----- Inner Orientation -----
PRINCIPAL_DISTANCE: 100.500000000000 [mm]
PRINCIPAL_POINT_X: -0.000000000000 [mm]
PRINCIPAL_POINT_Y: 0.180000000000 [mm]
SENSOR_AREA_WIDTH: 103.860000000000 [mm]
SENSOR_AREA_HEIGHT: 67.860000000000 [mm]
-----------------------------

----- Exposure Annotation Data -----
MID_EXPOSURE_CORRECTION: 0.000762000000 [s]
FMS_SENSOR_CODE: UCXp
FMS_CAMERA_PORT: 1
FMS_EXPOSURE_NUMBER: 4761
FMS_PROJECT: 1537TORO
FMS_AREA: 1537TORO_6cm_6030
FMS_LINE_NUMBER: 47
FMS_SEGMENT_NUMBER: 1
FMS_WAY_POINT_NUMBER: 0
GPS_DATE: 150412
GPS_TIME: 162517
GPS_LATITUDE: N43.745319 [degree]
GPS_LONGITUDE: W079.358953 [degree]
GPS_ALTITUDE: 1159 [m]
GPS_POSITION_SOLUTION: GPS
GPS_TRACK_OVER_GROUND: 73 [degree]
GPS_ABOVE_GROUND_LEVEL: 1008.7 [m]
GPS_GROUND_SPEED: 80.4 [mps]
------------------------------------

----- Level-3 Parameters -----
IMG_TYPE: High resolution multi channel RGBI
ROTATION: 0 [degree]
------------------------------------
TIFFTAG_RESOLUTIONUNIT=1 (unitless)
TIFFTAG_SOFTWARE=UltraCam Aerial Radiometry Core 12.6.1408.2501
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,11310.0)
Upper Right (17310.0,    0.0)
Lower Right (17310.0,11310.0)
Center      ( 8655.0, 5655.0)
Band 1 Block=17310x1 Type=UInt16, ColorInterp=Red
Band 2 Block=17310x1 Type=UInt16, ColorInterp=Green
Band 3 Block=17310x1 Type=UInt16, ColorInterp=Blue
Band 4 Block=17310x1 Type=UInt16, ColorInterp=Undefined


texts: Select a distinct attribute with maximum shape area

texts:
I am trying to manipulate a large dataset containing a lot of sliver polygons resulting from an overlay operations. (The eliminate tool is very time consuming to generalize them.)

I identified the sliver polygons using:

SELECT * FROM data1 WHERE
  [CODE] IN ( SELECT CODE
              FROM data1
              GROUP BY CODE
              HAVING Count ( CODE ) >1)

Now I want to select the slivers with maximum area from the ones isolated through above query.

I need to select each CODE value (from the multiple instances) where the associated shape area is the largest. This will help me eliminate the sliver polygons without having to merge them separately.

Any help is extremely appreciated.


texts: How can transform UTM Geographical coordinates from Zone 19 to Zone 18?

texts:
I want to transform UTM Geographical coordinates from Zone 19  to Zone 18 in southern hemisphere. How can I do it using QGIS?

5622274 -   249667  Zone 19S  Datum 1984

5623190 - 249646  Zone 19S  Datum 1984

I can do it in ArcGIS using Project (data managment) and saving the shape zone 19 as zone 18. In order to see the new coordinates, in atributte table I make 2 colummns , X and Y (float), then I use geometry to calculate x coordinate of point (meters) and the same for y on the second column.


texts: Can ArcGIS Server 10.1 GetFeatureInfo return the geometry?

texts:
Using GeoServer I can make a WMS GetFeatureInfo request using the application/vnd.ogc.gml info_format and it returns the full GML of the feature, including the geometry.  This is really useful as I can then use this GML in the browser to highlight features.

Example GeoServer request:

  http:///geoserver/wms?LAYERS=&QUERY_LAYERS=&STYLES=&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&BBOX=&FEATURE_COUNT=10&HEIGHT=500&WIDTH=500&FORMAT=image/png&INFO_FORMAT=application/vnd.ogc.gml&SRS=EPSG:2157&X=500&Y=500

However when I run this against an ArcGIS Server 10.1 WMS the result doesn't include the geometry.  I know there are ways of customising an ArcGIS Server 10.1 GetFeatureInfo response (http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//005300000228000000), but I don't see any way of getting it to return the geometry.

I've tried the same with all the ArcGIS formats included in the GetCapabilities, but none of them return the geometry coordinates, they just give the geometry type.  E.g. the GeoJSON gives "GEOMETRY": "Polygon"

Is it possible to get ArcGIS Server 10.1 GetFeatureInfo to return the geometry coordinates?


texts: Changing brightness of WMS layer in Leaflet?

texts:
I need to change brightness in images of a Leaflet wms layer. There is several plugins to do it like Leaflet.TileLayer.GL, Leaflet.TileLayer.Colorizr, Leaflet.TileLayer.PixelFilter, but they only works with TileLayer, not WMS layers.

Can I do this using css?


texts: Contour lines from AutoCAD to QGIS: how to automatically assign the right value to the Z

texts:
I have a plate of contour lines made on AutoCAD (contour lines do not have a Z). The Z is indicated only for some curves but in a text format.
If you import the table in QGIS the situation is this (the points correspond to where the Z is indicated):

My question is: knowing that the contour lines shown are every 20 meters is there a way to automatically assign the right value to each curve having only this little information?
I thought about merging attributes between the layer of the contourlines and the text layer, but this would solve the situation for only a few curves.


texts: Creating DEM Contour files from shapefile line data in QGIS

texts:
There are many question occurring on this forum regarding this issue, however the issue I am facing has not been answered in Generating DEM from contour Shapefile using QGIS?, which is that I have contour line data. Not point data.

Following the steps in Generating DEM from contour Shapefile using QGIS?. I have used SAGA Tool in the Processing Toolbox. Selected the Triangulation Tool from the Raster Creation Tools.

At this point It asks for the 'points' required for rasterising. But I want to input 'line data' does anyone know how to do this?

below shows the 'line data expressed in QGIS'

the purpose behind this is due to the fact I need to create a 'height profile' using the Profile tool to create a cross section area from contours as seen in the example below.

I am using QGIS V 3.4


texts: Google Maps v3 with Openlayers - First Zoom Level

texts:
I use Google Maps v3 in my application and I noticed a strange behavior, When I refresh the web page, the map zoomed out to display the whole map then get back to the zoom level I specify in my code. This happen fast but every time I refresh the page. How can I fix it?

EDIT

Part of mu code:

var geographic = new OpenLayers.Projection("EPSG:4326");
var mercator = new OpenLayers.Projection("EPSG:3857");
maxExtent: new OpenLayers.Bounds(-20037508, -20037508, 20037508, 20037508.34).transform(geographic,mercator);

map = new OpenLayers.Map('map',{
        projection: new OpenLayers.Projection('EPSG:3857'),
        displayProjection: new OpenLayers.Projection("EPSG:4326")
 });

 gmap = new OpenLayers.Layer.Google("Google",{ numZoomLevels: 20});
 map.addLayers([gmap]);


texts: Cannot see widgets and templates in ArcGIS Online?

texts:
I want to use and extend some of widgets and themes of web app builder in arcgis.com or ArcGIS Online, but I can't see them when I click on web app builder button for developing app.

How can I solve this problem?

Is this problem due to my account or my browser or something like this?


texts: Homebrew (OS X), gdal, and gdal2tiles.py/python

texts:
I had a previous installation of gdal (v. 1.11.5_3) that included in its binaries python files like gdal2tiles.py. I recently upgraded gdal with to 2.2.4 using Homebrew, but the result seems to have removed the python functions.

Is there an additional step that I should take?


texts: Creating dynamic, interactive map with QGIS?

texts:
How can we create an interactive map with QGIS?

I would like to keep aerial photographs as a background layer. When the mouse moves over elements of a shp layer, I'd like them to be highlighted and that a windows shows the information contained in the attribute table of my layers. I would also like that the export is done in a good resolution to be able to zoom in and zoom out dynamically.

  Comment peut-on créer une carte interactive avec Qgis ?

  je souhaiterais conserver en fond des photos aériennes. Lorsque la
  souris passe sur des éléments d'une couche, je souhaiterais qu'ils
  apparaissent en surbrillance et qu'une fenêtre indique les
  informations contenues dans la table attributaire de mes couches shp.
  Que l'export permette se fasse dans une bonne résolution pour pouvoir
  zoomer et dé-zoomer de manière dynamique.


texts:
This is only part of a bigger code block that I'm building.

The goal of these lines is to import a new layer, and select one feature from the feature class(TOWNSSURVEY_POLYM). The code imports the .shp as a layer feature successfully, however, it throws an error after running (regardless of whether or not the I included the bit about selectLayerByAttribute_management() ).

Here's the error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 13, in <module>
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\mapping.py", line 48, in AddLayer
    assert isinstance(data_frame, DataFrame)
AssertionError

If I'm not mistaken, it seems to me as though the problem has to do with the my dataframe not being accepted as a string. Here's my code block.

I am still new to ArcPy.

from arcpy import mapping

## setting variables to parameters

map = mapping.MapDocument("CURRENT")
df = mapping.ListDataFrames(map)[0].name = "Data Frame"
data = 'C:\Users\GIS\FilePath\TOWNSSURVEY_POLYM'

## adding the new layer as a '.lyr' file, as ArcMap won't allow you to directly import '.shp' files

lyr = arcpy.MakeFeatureLayer_management(data)

#mapping.AddLayer(df, lyr, "AUTO_ARRANGE")

arcpy.SelectLayerByAttribute_management(mapping.AddLayer(df, lyr, "AUTO_ARRANGE"), NEW_SELECTION, ' [TOWN] = "ASHBY" ' )

arcpy.RefreshActiveView()
arcpy.RefreshTOC()


texts: When creating graduated color legends, other than manually formatting value ranges under the symbology legend label, is there any way of automatically showing the full numeric value range in the legend rather than the less than or equal to upper value range?


texts: gdal2tiles for both raster and vectors with a directory of input files (shps and .tifs)

texts:
I am looking to tile a large dataset that was split into raster areas based on political borders (ie. 1 raster file for each area). I have a few hundred .tiffs and shapefiles for the area and am trying to find the easiest way to tile them rapidly without needing to merge them first (would be an unruly large .tiff and .shp) likely why they were originally split. I have to do this for both datasets the large series of .tif's and the .shps.

I have tried QTiles but it runs incredibly slowly on a high end machine for tile generation and am hoping to use something like gdal2tile_parrelel but it doesnt work through a directory. Anyone had a similiar struggle and find a solution?

machine is i7-7700k, 64gb ram and Raid 0 drives.


texts: Displaying NLCD landcover data in QGIS with landcover class type?

texts:
I am trying to display land cover layer (.tif) into QGIS. But some reason the layer appears in grayscale when I load in QGIS (Figure 1). I am expecting to see the somewhat similar display as in ArcMap (Figure 2). In QGIS, when I looked into the layer's properties, the Singelband gray is a default Band rendering for this layer. I wonder how can I display the layer as in AcrMap or is there steps in QGIS to process the layer to look like in Figure 2.

I'm using
QGIS 2.14.3-Essen
OS: Windows 8

 Figure 1

 Figure 2


texts: Querying on Household Income ranges in a dataset

texts:
I've recently been assigned a project to create an ArcGIS database to hold responses from a survey which will be used to perform analysis on the data. One column requires me to enter a dollar-range for household income. The possible ranges are:

25,000 //for > $25,000
25,000-49,000
50,000-99,000
100,000-149,000
150,000-199,000
200,000 //for >$200,000

When trying to select by attributes, for example

Select * From Respondent Where HHI > '150,000-199,000'

everything is returned except for the 100,140,000 and 150,000-199,000 ranges is returned.

This makes sense because it's not less than or equal to. I'm more familiar with regular SQL, not esri's version, however I'm thinking it's not matching the full range, and just performing the > operation on the first character in the range (or hanging on the -), and since the ASCII value is less than that of other ranges, 25,000-49,000, for instance, they are assumed to be greater than the initial query 150,000-199,000.

Does anyone know a way to properly select on these? This has had me stumped for a while now.

Edit: Here's an example image showing how the data is entered into the specific column.


texts: gdal_translate or gdalwarp: Preserve noData value when using DISCARD_LSB option

texts:
I want to reduce the file size of a 32bit floating point 1-Band GeoTIFF file that currently uses -99 for noData pixels. The relevant part from gdalinfo -mm:

Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
    Computed Min/Max=-5.101,3.317
  NoData Value=-99

I want to use the DISCARD_LSB option to throw away the least significant bits:

gdal_translate -co COMPRESS=DEFLATE -co DISCARD_LSB=16 input.tiff output.tiff

But then the gdalinfo -mm gives:

Band 1 Block=36000x1 Type=Float32, ColorInterp=Gray
    Computed Min/Max=-99.250,3.320
  NoData Value=-99

In other words: The original noData pixels got treated as normal data in the conversion process (Minimum data value of -99.250).
How can I avoid this?
The same behaviour occurrs when I use gdalwarp instead of gdal_translate. The version is 2.4.0.

An uncompressed sample image can be found here.


texts: Attributing height data from points to Polylines for Tin creation?

texts:
I have to make a TIN of water surface height from about 170 points which have been collected from RTK along a river. The points however are few and far between and rarely on the opposite banks to allow creation of triangles to create a decent TIN. This therefore produces a TIN with gaps. I should mention that i have created a polygon of the water area.

What i have been asked to do is create lines from one point (which has x,y,z data) to the other side of the river and connect it to the polygon (this does not contain x, y, z data). I then need to attribute the z value from the point to that line. This will then give me a water surface height. I have to do this to all the points and i cant figure out what to do apart from going through each value in excel and adding it manually.

So far, i have created a shapefile with polylines, and this is where i am stuck.

I don't understand how to use Python and my knowledge of ArcMap is relatively elementary.

Can you explain step-by-step how it can be done?


texts: Is Esri high resolution (0.3m) imagery service still available?

texts:
Neither the World Imagery basemap in ArcMap nor the service at ArcGIS Online provide high resolution (0.3m) images anymore.

Has anybody been able to access the 0.3m resolution imagery recently?

Here is the longer question statement:

According to the metadata here, Esri High Resolution Imagery
should have 0.3m resolution imagery in the continental United States. However, neither that or the new alternative, World Imagery (Clarity) look like 0.3m imagery. Below are example images from Massachusetts, USA, near coordinates (-70.655505, 41.515357). The first two are from World imagery, and the other are from Clarity. Note the scale bar in the bottom left corner. The resolution does not seem to be 0.3m.

Here is a comparison of a 1m resolution imagery from another source (top) to the World Imagery for the same area (bottom). In the top image one can clearly see the raster cells of size of 1m. Not the same for the Esri image. It is blurred and does not seem to be 0.3m resolution in contrast to what the image info states.

Can anybody successfully see a 0.3m resolution imagery as a basemap or through any of the services?

Two ways to reproduce this issue:

You can follow the Esri High Resolution Imagery link to reproduce this issue online.
Or in ArcMap go to menu Add Basemaps -> Imagery.

Here is a list of posts pointing to the same problem:

Esri basemap no high resolution
High resolution basemaps not loading
High resolution imagery service not available


texts: QuantumGIS EPSG:3857 doesn't match with Leaflet

texts:
Leaflet documentation said that they have EPSG:3857 as default projection.

Using QGIS with OpenLayers plugin, I point to a location to find it's coordinate. I used EPSG:3857 in QGIS as well.

The problem is, the location doesn't match. I know that the location in Leaflet using "lat-long" format, while location in QGIS is using "long-lat" format.

Leaflet documentation said that it used "degree", while (I don't quite sure) I think QGIS use decimal.

Does anybody has experience with this? What is the solution?


texts:
I updated my QGIS 2.18.11 version with latest QGIS 2.18.23.

I have a problem, I can open my old projects, but when I try to zoom in or zoom out with the mouse wheel, everything disappears, and QGIS  tries to refresh the screen, which never refreshes. When I try to close QGIS,  not responding and freeze screen.


texts: Triggering the retrieval of a feature on a feature layer by something other than click using OpenLayers 2

texts:
I have added a feature map to my OpenLayers 2 map, by using the following code. When I click a spot on the map, it correctly retrieves the feature name for that place according to the feature layer.

Code:

       featureControl = new OpenLayers.Control.WMSGetFeatureInfo({
           url: featureLayer.url,
           title: '',
           layers: [ featureLayer ],
           infoFormat: 'text/xml',
           queryVisible: true
       });
       featureControl.events.register('getfeatureinfo', this, function(event) {

           var xml = new window.DOMParser().parseFromString(event.text, 'application/xml'),
               path = featureLayer.$featureInfoXPath,
               featureName;
           try {
               featureName = xml.evaluate(path, xml, null, XPathResult.ANY_TYPE, null).iterateNext().value;

                console.log('got the feature:', featureName)
           } catch ( error ) {
               console.log('error', error);
           }
       });
       map.addControl(featureControl);
       featureControl.activate();

However, I would like to trigger this retrieval of the feature name on other events as well, such as when I type in an address in an input field and I select an address in an auto complete. So I am looking for a mechanism to trigger the GetFeatureInfo event by something other than "click".

Is this possible, and if so, how?


texts: How to create Online IDW(Inverse Distance Weighted) Interpolation?

texts:
I have some point data  which store as JSON.As follow:
                [["Lat":"12.12","Long":"23.23","ZValue":"123"],...].

I tried to interpolate them as contour line and displayed them on my website ? I do know how to interpolate data in ArcInfo,so I tried to publish a geoprocessing  task and using ArcGIS JavaScript API to generate them but it failed to published to ArcServer.

However,I search for a long time if there exits some other solution but found nothing about how to implement this.

Thanks for any suggestion!


texts:
I'm trying to update a given field in this case the last one of the table by concatenating all values in the row excepting null values or '' values.

So far when execute the code until line: print (row[ListnamesLeng-1]) i get the concatenation but has been impossible to update the field. I get the message StopIteration: iteration not started. i've tried put cursor.updateRow(row) inside and outside the loop but i keep receiving the same message.

Sorry about my english im still learning! ;)

import arcpy
arcpy.env.overwriteOutput = True

arcpy.PivotTable_management(in_table="C:\Users\daniel.cortes\....", fields="idRegistro", pivot_field="LLAVE_NESTLE_TIENDAS_VENDEDOR", value_field="LLAVE_NESTLE_TIENDAS_VENDEDOR", out_table="C:\Users\daniel.cortes\.....")

FieldListnames = [f.name for f in arcpy.ListFields("C:\Users\daniel.cortes\..")]
ListnamesLeng = len(FieldListnames)
print (ListnamesLeng)
fielditerator = range(2,ListnamesLeng-1)
print (fielditerator)

arcpy.AddField_management("C:\Users\daniel.cortes\..","LLAVE_NESTLE_TIENDAS_VENDEDOR_F", "TEXT",field_length=300)

with arcpy.da.UpdateCursor("C:\Users\daniel.cortes\...",FieldListnames) as cursor:
    for row in cursor:
        ConcList = []
        for i in fielditerator:
            if row[i] != '':
                ConcList.append(row[i])
        row[ListnamesLeng-1]=str(ConcList)
        **print (row[ListnamesLeng-1])**
    cursor.updateRow(row)


texts:
I have the following requirements to manage raster data

1) Convert raster-maps from BSB format to PNG or GeoTIFF
2) Cut the borders of the map
3) Merge individual maps to form a contiguous coverage map
4) Tile them for caching and loading
5) Serve them up on a WMS service to clients.

I am using OpenLayers->GeoServer->PostGIS stack to manage all my vector data. I have been thinking of loading, storing & retrieve raster data the same way (in PostGIS).

So my questions are

1) What is the benefit of storing rasters in DB over doing this in the file system?
2) What are the steps to achieve this and what tools/ commands I would use?
3) Would I leverage the DB cache or would I still need a tile-cache ?
4) What are the bugs/ gotchas in PostGIS raster handling that I need to watch out for?

Thanks,
Ramesh


texts: Can SQLite and SpatialLite files be linked

texts:
Coming from a ms-access background, it is common to split up larger applications into several database files.  An ms-access database can then access data that is in another database via a linked table.  This allows for better information management in many cases.

I am wondering if I can do something similar with SpatiaLite/SQLite, for example use Spatialite for the GIS part of projects and a separate SQLite database for the non-spatial components.  One of the reasons for doing this would be to hide the Spatialite system tables from the user.  The other might be that it simply makes sense to separate datasets into different databases (to maintain a single point of truth, etc).  For example, my spatialite database may have a list of plant locations recorded during a field survey but my SQLite database may have a data dictionary with all the plant names for my country (~40 000 species).  If linking tables is not something that SQLite does, what is the alternative?

Andrew M


texts: Planet CLI - problem with authorization in planet init

texts:
I try to use Planet CLI or pclip plugin for download a Planet data. But I've got a problem with authorization in the client.

When I try a command planet init or pclip planetkey in the command line and enter email and password.

I get an error:

  Unexpected response: invalid parameters, login process has changed

How do I resolve this problem?


texts: Resampling with Python console (QGIS): number of parameters

texts:
I am trying to resample two raster layers with Python console in QGIS with this code:

outputs_0=Processing.runAlgorithm("saga:resampling", 'G:/DOCTORAT/2014/Incendis/Dades/LST/LST_TERRA/Tiff/2012/UTM30N/2012153_day.tif', True, 0, 0, 0, "-81084.51, 1138915.49, 3972094.6, 4861094.6", 1000.0, 'G:/DOCTORAT/2014/Incendis/Dades/LST/LST_TERRA/Tiff/2012/FINALS/2012153_day_FINAL.tif')

I got this:

  Error: Wrong number of parameters

However, it seems to contain all needed parameters:

INPUT <ParameterRaster>
KEEP_TYPE <ParameterBoolean>
TARGET <ParameterSelection>
SCALE_UP_METHOD <ParameterSelection>
SCALE_DOWN_METHOD <ParameterSelection>
OUTPUT_EXTENT <ParameterExtent>
USER_SIZE <ParameterNumber>
USER_GRID <OutputRaster>

Someone could please help me?

Thanks!


texts: Using ArcPy or SQL methods to truncate data from ArcSDE feature classes and tables?

texts:
There are a handful of ArcPy methods for truncating data in a table. It's also possible to do a direct SQL command to truncate a table (e.g. TRUNCATE TABLE DBO.VALVES) via a Python library (e.g. cx_Oracle or pyodbc).

What are advantages or disadvantages to ArcPy methods over SQL alternatives?

Do these considerations change when data is versioned or not?

If it's relevant, right now we're working with Python 2.7 and database 10.2.1; we're transitioning to Python 3 and database 10.6+ sometime in the next year.

Note: For the purposes of this question, the ability of Delete Rows to work on a subset is not relevant — we're interested in completely emptying the table to repopulate with data, and debating what method(s) would be best to standardize on.


texts: Converting continuous to discrete data using ArcGIS Desktop?

texts:
Does anyone know how to convert continuous data (in this example a raster file) into discrete data in ArcGIS Desktop?

I need to get the values from each pixel (emissivity from satellite image) into a table, so I can join them with another layer.


texts: GRASS GIS not importing any Shp File - UnicodeEncodeError

texts:
I tried to work on a project to get isochrones with GRASS Tools. I keep getting error messages in QGIS so I tried using GRASS as a standalone but I get the following error message when I try to import any layer from any project I have on my Harddisk:

UnicodeEncodeError: 'ascii' codec can't encode character u'\xfc' in position 25: ordinal not in range(128)

The only thing that is changing it the position number in this error message. What is going on? My shapes are in UTF-8.


texts: Why is my spatial join count 0 in ArcGIS Desktop?

texts: In arcMap, I have a list of locations with lat/long. From a previous post, I learned about the Spatial Join tool to determine the amount of joins that I have. I tried to use the tool with my states layer and it worked perfectly. However, when I tried to do it on the county layer, it is always giving me with 0 joins. What is the cause of this? The county layer that I am using is from http://www.gadm.org/. The one that I used for the states came from arcgis online. There was one for the countries, but it would give me an error when trying to do a spatial join on the data.


texts: Why does an OpenLayers map.getNumLayers() return 0 even though I've added a layer?

texts:
When I try the following code, OpenLayers returns a zero result. What have I done wrong, here?

var map;

function loadMap () {

    // Based on example code: view-source:
    //     http://openlayers.org/dev/examples/vector-formats.html

    map = new OpenLayers.Map(el);
    var googleStreetMap = new OpenLayers.Layer.Google('Google Street Map',
                                                      { numZoomLevels: 20 });
    map.addLayers(googleStreetMap);
    var result = map.getNumLayers(); // Result == 0 :sad-panda:
}

In the code above, I thought that the map.getNumLayers() should return 1 - considering i've asked to add a google street map.

Why is it returning zero, please?

UPDATE: something else to note: there is also no baseLayer :( I thought this would be the base layer???


texts: I've been blindly using Corpscon to convert from NGVD29 to NAVD88.  Have any of you kind folks found issues with your outputs?


texts: Reading a raster map cells values in C#?

texts:
Suppose I have a raster map file (TIFF, TGA, JNG... etc) loaded with geo-info manifested by the cells values.
I need to read the values of the raster map cells "programatically" in a C# program.

P.S: I don't want to use "heavy" tools or libraries (like ArcGis and ArcObjects) for such a simple task.


texts:
I'm trying to install gdal on my new computer in the anaconda environment using conda but having issues. All I've done so far is install Anaconda for python 3.6 and install gdal as described here + a few other libraries using both conda and pip.

Trouble is when I try to run from "osgeo import gdal" the iPython console in spyder I get this error, which I'm not smart enough to interpret on my own:

        Traceback (most recent call last):

        File "", line 1, in
           from osgeo import gdal

        File "/anaconda/lib/python3.6/site-packages/osgeo/init.py", line 21, in
           _gdal = swig_import_helper()

        File "/anaconda/lib/python3.6/site-packages/osgeo/init.py", line 17, in swig_import_helper
           _mod = imp.load_module('_gdal', fp, pathname, description)

        File "/anaconda/lib/python3.6/imp.py", line 242, in load_module
           return load_dynamic(name, filename, file)

        File "/anaconda/lib/python3.6/imp.py", line 342, in load_dynamic
           return _load(spec)

        ImportError: dlopen(/anaconda/lib/python3.6/site->>>>packages/osgeo/_gdal.cpython-36m-darwin.so, 2): Library not loaded: >>>>@rpath/libicui18n.58.dylib
         Referenced from: /anaconda/lib/libgdal.20.dylib
         Reason: image not found

Is there a problem with the conda install of gdal I'm using, or am I missing a piece?

Is there a different/easier way to start using gdal than this?


texts: Permission denied when creating shapefile

texts:
I'm having problems creating new shapefiles - in QGIS 2.14.20 any new fields I create and any associated data are not saved.  I've now got the same problem in QGIS 3.0, but there I get a message saying layer creation failed.  There is an error message but it just seems to say 'warning  :'   But under OGR I get:

2018-05-10T10:44:30    WARNING Creation of OGR data source test.shp failed: Failed to create file test.shp: Permission denied

I am running on an iMac High Sierra 10.13.4.  I also sometimes can't save projects.  I am the only user of this machine and have admin privileges, and am trying to save the shapefiles and projects locally.


texts: Python Code to return column data in multiple feature classes and gdb's

texts:
>>> import arcpy
... in_table = r"K:/MAPPING/PA/LEIDY/REPL/0188_18/MDL/GIS/MASTER/0188_18.gdb/DP_Centerline"
... with arcpy.da.SearchCursor(in_table, ["last_edited_date"]) as scursor:
...     for srow in scursor:
...         print srow[0]
... del srow, scursor
...

This works for DP_Centerline "last_edited_date" column and I am trying to PRINT "last_edited_date" across multiple feature classes within multiple geodatabases.


texts: How to use overview map with atlas plugin in QGIS 2

texts: I'm using atlas plugin to batch generate maps. I insert an overview map but it doesn't get updated over iteration and generation of maps.Is it possible? If it's not, is there any alternative?


texts:
I have a lot of widgets and my dashboard in Flex Viewer 2 has gotten cluttered. How can I add menus to the dashboard rather than buttons to organize the widgets?
Any help will be appreciated.


texts: What are arc's guarantees on data precision, XYTolerance (i) when displaying (ii) in arcpy?

texts:
Interesting behaviour:

Two line segments in a dataset have almost, but not quite, identical endpoints.  The difference between them is less than the XYTolerance for the data set.

Reading the coordinates into arcpy, they come out as different to one another.

But if I zoom in and in and in, far beyond the XYTolerance, ArcMap still displays them as identical - no gap appears between the lines.

Is this to be expected?  Does arc make any guarantees on precision of data drawn on screen, or data passed to arcpy?  What will happen if the points are separated by the XYTolerance exactly?

Update:  I tried to test the boundary case, but the plot thickens.  The following script creates a feature class that displays correctly in ArcMap despite having features smaller than the cluster tolerance.  Still I have other data that doesn't display correctly under the same circumstances.

import arcpy

out_polyline_feature_class = arcpy.GetParameterAsText(0)

# t is tolerance
t=0.001

ic = arcpy.InsertCursor(out_polyline_feature_class)

edges = [ [ (0,0), (0,t/2) ],
          [ (0,0), (t,0) ],
          [ (0,0), (-t,-t) ],
          [ (0,0), (-2*t,0) ]
        ]

line = arcpy.Array()
pnt = arcpy.Point()
for edge in edges:
    line.removeAll()
    for point in edge:
        pnt.X, pnt.Y = point
        line.add(pnt)

    row = ic.newRow()

    row.shape = line
    ic.insertRow(row)
    del row

del ic


texts: Georeferencing Himawari-8 in GDAL (or other)

texts:
Does anyone know the appropriate project information to georeference (and hence reproject) Himawari-8 AHI data in the IR?
The images are 5500x5500 pixels.

Currently I try these:

ulx=-5570248.832537
uly=5570248.832537
lrx=5567248.429179
lry=-5567248.429179
gdal_translate -a_srs  "+proj=geos +a=6378169 +b=6356583.8 +lon_0=140.7 +h=35785831" -a_ullr $ulx $uly $lrx $lry HDF5:"TEST"://img tmp.tif

gdalwarp -srcnodata -32767 -s_srs '+proj=geos +lon_0=140.7 +h=35785831 +x_0=0.0' -t_srs '+proj=latlong +datum=WGS84' -tr 0.03 0.03 -te 55.7 -80 225.7 80 -order 3 tmp.tif FINAL.tif

This, however, results in an image that's not quite right.
Australia and Indonesia almost match a coastline boundary overlay but as one moves further from the subsatellite point there exists a larger and larger offset. See this comparison image:

If someone could point me to the correct code to use for reprojection that would be great. I suspect it's the a_ullr values that I have wrong.


texts: Calculating coordinates given a reference point and distance

texts:
I am using QGIS 3.6 and trying to calculate coordinates (lat&long) from a point given the distance. I have the lat&long for reference point and have the distance from the reference point.

Points A, B, and C are reference locations of a road

X, Y, and Z are building within a distance from A toward point B
How can I get the lat and long for X, Y, Z?


texts: Select MapImageLayer along with FeatureLayer

texts: I have FeatureLayer on top of MapImageLayer so both layers are intersecting each other on some points. I have to show popup on map click which will contains all features available on click point. It means popup will have multiple pages. Right now when I click on map I get only feature layer's value and popup also contains only feature layer's value. Please help me how can I achieve this?


texts: arcpy geometry difference method fails randomly, not reproducible

texts:
Has anyone used the arcpy geometry difference operator and seen unexplained errors like what I describe below?

This is a weird one. I have a script that calculates wave fetch distance in the ocean using a polygon file with land features (coastal areas and islands) and a point file with points for which fetch is to be calculated. It does this by creating a multipart geometry object with numerous individual line segments in a star-shaped pattern radiating around the point (these are bearing lines at user-specified intervals). I am "erasing" this pattern of lines using the geometry.difference method described here: http://bit.ly/1vBMSjP (scroll down to "difference" under the "Methods" section). The object that is used to erase the lines is a multi-part polygon geometry created from the user-defined land features. Both geometries are in the same projection, which is the same as the coordinate system defined in arcpy.env.outputCoordinateSystem. After completing the geometry difference operation, the script moves on to the next point and generates a new set of line geometries. All in_memory objects and variables are erased after each iteration to free up memory.

Randomly, the script throws this error:

Traceback (most recent call last):
      File "c:\data\scripts\py\FetchGeomV3.py", line 182, in <module>
    lineErased = currentLine.difference(landGeom)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\arcobje
    cts.py", line 793, in difference
    return convertArcObjectToPythonObject(self._arc_object.Difference(*gp_fixarg
    s((other,))))
    ValueError: <geoprocessing describe geometry object object at 0x02447640>

The really weird part is that if I run it again, using the exact same input files and parameters, I cannot reproduce the error. Sometimes it runs all the way through without ever throwing an error, sometimes every 100 points or so will generate the error--on different points than on the previous run.

Original code is here: http://bit.ly/1rtMNtc


texts:
I have a script in QGIS which will not add to the User Scripts when I choose Add script from file.

It does, however, drop into the 'scripts' folder. It still does not display though.

Any ideas?

Happens on both QGIS 2.18.17 and QGIS 3.0, but does not always happen on different machines.


texts: How do you create either a geotiff or a .tfw in QGIS using python?

texts:
Probably I just don't understand the API documentation, I don't use object oriented programming to its full capability.

I do see this has to be done from a composer object in order to automatically apply the on-the-fly projection transformations since the input layers have different CRSs.  (Right?)

I found the c++ code that generates the worldfile, but, how do you generate the worldfile from the python console?

More generally, how do you export from a composer to geotiff (or, alternately, write a .tfw worldfile) using python?


texts: Arranging GUI elements in QGIS 3 processing model user interface

texts:
Using QGIS 3 graphical modeler I created a model with a variety of input parameters (> 10). When starting the model, a nice user interface is generated. But the arrangement of the GUI elements is totally confused and apparently unpredictable during creating the model. It does neither seem to depend on the arrangement of input parameters in the graphical modeler nor on the chronological order in which the input parameters are added to the model.

Is there a way that let's me arrange the GUI elements in the order I desire?


texts: How to create a convex hull with Leaflet points?

texts:
Using leaflet, given a set of points, I need to represent the area (polygon) containing all points.
It's a similar approach of building a bounding box but the result instead of being a rectangle it's a polygon with the external points as vertices.

Perhaps this problem can be fixed on the server side using arcgis?


texts: Rasterization tools for vector datasets with styles

texts:
I have hundreds of Personal Geodatabases and I want to create one georeferenced TIFF for each one. Each PGDB has the same feature classes ('CotourLines', 'Hidrography', etc) for what I have some layers with the styles. The style (or symbol) for each feature could depend on many fields. And the final TIFF files must have a fixed pixel size.

I was searching for tools I can use. From ArcToolBox there are some tools like 'Feature to Raster' but doesn't fit to my requirements.

Maybe should I make an automation?

Having an MXD containing the layers with the styles, loading the feature classes of each PGDB, and using 'Export Map' to create the TIFF file?

This is the only solution I have in mind.

I think there should be an straightforward solution I could use, maybe with an ArcGis extension. I could export the PGDB to ShapeFiles (or a PostGIS database) and have the styles in SLD files.


texts:
I loaded the following multipart polygon in QGIS which consists of three separate polygons. The highlighted (yellow) polygon consists of the bigger upper part and the smaller lower part.

My aim is to render the same multipart polygon in an Android app. Therefore, I extract the points of each polygon. Then I render the coordinates (latitude, longitude) in the app.

However, the Android rendering shows up with an error as shown in the following screenshot:

It looks like there is a false connection between two points. I cannot see this line in QGIS - neither in the polygon nor in the points layer.

Is there any way I can detect this error in QGIS and fix it?


texts: TauDEM: How to extract the longest flow path from an outlet to the basin divide as a shapefile

texts:
Using TauDEM I've created a number of derived raster products from an initial 30m SRTM DEM with a statewide (South Carolina) extent. Using the D8 modelling methods, I've created both the D8 Flow Direction and the D8 Flow Accumulation secondary inputs. Next, using the Gage Watershed method, I've rendered the catchment above a pour point. However, I'm struggling to create an output of the longest watercourse in the basin.
I have executed the Grid Network method, which creates several outputs, including a "grid of the longest flow length upstream of each point" (PLEN), as well as Strahler Order raster. I have also executed the Stream Reach and Watershed method (i.e. StreamNet), which also returns several outputs, including a vector of the stream links that also includes a Strahler Order attribute.
While I know I can get a value for the distance to the divide using any cell from the PLEN raster, what I want is just the longest flow from my outlet up to the divide as a discrete result--raster or vector doesn't matter. I'm not seeing a way to derive this result using just the TauDEM methods, which makes me think I'm overlooking or misunderstanding the language describing the proper method.
Is there an existing method in TauDEM that renders the longest flow in a basin as a single, discrete feature/value?
(FWIW, I tried using the LongestFlowPath method from WhiteboxTools just for this one step in the workflow, but it fails to run to completion. At this point I'd prefer to stick with TauDEM for the solution.)


texts: Choosing between Dynamic layers and Feature layers

texts:
I'm building a client using the ArcGIS WPF libraries, that's displaying my entities on a world map. I'm using an ArcGISDynamicMapServiceLayer to display the geographic and topographic base map, and I'm now wondering what's the best way to show my custom entity layer.

At first, I used a FeatureLayer object, linked to the layer's published URL, to display the entities on the map. This worked great in that it showed the entity symbology and allowed me to interact with it (catch right-clicks, select/deselect), but it doesn't show the labels associated with the entity, just the symbology.

So I tried using a second ArcGISDynamicMapServiceLayer to display the entities, which gives me data labels, but without the interactivity.

So now, since these are both client requirements, I find myself using BOTH, and synchronizing visibility between the two instances of the layer, once as a FeatureLayer, and once under a dynamic layer. This is not a good practice, and bugs rear their head constantly.

So my question is - what's the best way to achieve those two goals? Is there a simple way to add my own data labels on the map for each entity? There could be hundreds of them on 5-6 of these entity layers, and this feels like replicating functionality that ArcGIS already offers. On the other hand, I can simulate the right-click-on-entity behavior using a DynamicMap, catching the MouseClick event, and running a Query action on the map to find the entities clicked, but this, again, is reinventing parts of the wheel, and leaves some unsolved issues.

So what would you do? How do I get those two requirements together without ugly hacks?


texts:
I have been trying to install Cartaro the drupal GeoCMS for a while without success. I have properly installed Geoserver (running fine), Postgresql and postgis 1.5 enabled but when i try to install Cartaro it does not detect PostGIS version and neither does it allow me to change anything manually.

As this is a new software and the community and support is minimal still i do not know what is the right approach to solve this issue.

Any advice would be helpful.


texts: I have two separate polygons A and B in ArcMap. Which tool I have to use in order to find the mean line between them?


texts:
I have a simple CSV of {lat,lon,text} and I want to export it as a WGS84 Shapefile.

How can I achieve this?


texts: IR/R could not detect when using landsat 8 in ERDAS 2014 Is there a alternative way around?

texts:
I'm using ERDAS 2014 and I upgraded from ERDAS 2013.

I am doing a lot of band indices for various kinds of such as clay, ferrous, Mineral composite and so on.. from landsat 8.

Using the indices in ERDAS 2014, the problem is that it does not detect IR/R (infrared divided by red) if you change the sensor to Landsat 8 MS but if you look at other landsat such as 4 TM and 5 TM.

I can see that the formula for IR/R is IR/ Red which is Band 5 / Band 4. This band is for Landsat 8.

OR would I have to take a BAND 5 and BAND 4 separately and run the formula in ERDAS ?

Anybody have that same problem ?

Edit:

Recall this is a landsat 8 and there should be a list for IR/R but I do not find one.

Your thoughts ?


texts: ogr2ogr imports GML data (not from OSM) to Postgis but without the other_tags "hstore"-column

texts:
The problem

I'm encountering some issue trying to use ogr2ogr to import some GML input data (not OSM) into Postgres. I'm trying to use COLUMN_TYPES=other_tags=hstore to store all the fields not referenced in the -select option in a hstore-type-column.

The data gets imported, but there is no column of type hstore in the resulting table containing the rest of the fields...

The command I try to use

ogr2ogr -f PostgreSQL "PG:dbname=rem host=192.168.XX.YY user=XXXXX password=XXXXXXX" test.xml -append -nln public.t_test_ogr2ogr_v26 -select gml_id,identifier,beginnt -lco COLUMN_TYPES=other_tags=hstore -nlt geometry

The current result

'ogc_fid';'wkb_geometry';'gml_id';'identifier';'beginnt'
1;'0102000000030...';'DERPTL4a0000zZFt';'urn:adv:oid:DERPTL4a0000zZFt';'2015-11-25T14:28:13Z'
2;'0102000000020...';'DERPTL4a0000AIzO';'urn:adv:oid:DERPTL4a0000AIzO';'2015-11-25T14:28:13Z'
3;'0102000000020...';'DERPTL4a0000AIzP';'urn:adv:oid:DERPTL4a0000AIzP';'2015-11-25T14:28:13Z'

The expected result

The same as above with an additional column of type hstore.

'ogc_fid';'wkb_geometry';'gml_id';'identifier';'beginnt';'tags'
1;'0102000000030...';'DERPTL4a0000zZFt';'urn:adv:oid:DERPTL4a0000zZFt';'2015-11-25T14:28:13Z';'"anlass"=>"3000", "advStandard"=>"DLM50", ...'
2;'0102000000020...';'DERPTL4a0000AIzO';'urn:adv:oid:DERPTL4a0000AIzO';'2015-11-25T14:28:13Z';'"anlass"=>"3000", "advStandard"=>"DLM50", ...'
3;'0102000000020...';'DERPTL4a0000AIzP';'urn:adv:oid:DERPTL4a0000AIzP';'2015-11-25T14:28:13Z';'"anlass"=>"3000", "advStandard"=>"DLM50", ...'

Test Input GML-File (55 KB)

https://haste.zertrin.org/raw/tafomudifi

(if needed, the XSD-files can be found here: http://repository.gdi-de.org/schemas/adv/nas/6.0ea/ )

Corresponding gfs file created automatically by ogr2ogr

https://haste.zertrin.org/raw/diwowideji


texts: How to fix AttributeError: 'QgsVectorLayer' object has no attribute 'provider'?

texts:
When I try to start this algorithm I get an error:

Traceback (most recent call last):
  File "C:/Users/trdina/.qgis//python/plugins\sextante\core\GeoAlgorithm.py", line 145, in execute
    self.processAlgorithm(progress)
  File "C:/Users/trdina/.qgis//python/plugins\sextante\algs\ftools\SumLines.py", line 85, in processAlgorithm
    spatialIndex = utils.createSpatialIndex(lineLayer)
  File "C:/Users/trdina/.qgis//python/plugins\sextante\algs\ftools\FToolsUtils.py", line 32, in createSpatialIndex
    provider = layer.provider()
AttributeError: 'QgsVectorLayer' object has no attribute 'provider'

I have no clue woh to solve this error. Any help is welcome.
Best regards


texts: Rotating point object so that it is perpendicular to line on other layer in ArcGIS Pro?

texts:
I need to Rotate point object so that it is perpendicular to line on other layer.

need to line up valve symbol(points) perpendicular with the underlying watermains (polyline)

on ArcGis Pro


texts:
Working with QGIS, I saved my Rule Based Label as a .qml file. Opened a new Project Loaded the Shapefile and then Loaded the Style .qml file but it lost its rule based styling, x and y displacements and its Priority. Where am I going wrong?came in without the rule settings

This is the Attribute file:

Rule1

Rule2


texts: Layer from difference gets clipped when exporting

texts:
I have a shapefile of the US (5m resolution from the US Census Bureau), but I need the entire globe except for the US. To do this, I opened the shapefile (which I will call usa_layer) in QGIS, and then added a rectangle for the entire globe:

Vector > Research Tools > Polygon from layer extent

I then added a buffer to make sure that the shape would include the poles:

Vector > Geoprocessing Tools > Fixed distance buffer
Input layer: extent_layer
Distance: 10
Segments: 5

Next, I took the difference of usa_layer and  extent_layer:

Vector > Geoprocessing Tools > Difference
Input layer: extent_layer
Difference layer: usa_layer

The resulting layer looks like this:

This is what I was going for. However, when I right-click this layer and Save As... a KML, that layer looks like this:

Obviously, this is not the same as the difference layer. When I show this layer on my Google map, I see this:

This is even worse. It seems as if the edges that were not present on the KML layer (blue) but were on the difference layer (purple) are the only ones shown. Comfortingly, on the right of this sort of capsule shape, one can see that Alaska is indeed cut out of it as intended.

Why is this weirdness happening, and how can I fix it?


texts: If i clip a raster in Arcmap (10.5) and the export it so as to import it into autocad 2018(vanilla using georefimg - not map3d or civil).  it gives me a error message and will not attach the image.   it also will not open in adobe lightroom or infraview.  If I open it in QGIS and save it from there it works... Any thoughts on why?  would love to skip this step and keep all of my processing in arcmap.


texts: Filter collection with millis data on MODIS

texts:
I'm trying to split a modis collection in two others, using like a splitter a data time that's defined in millis. I don't want to do it manually.

var maskMOD13Q1 = function(image) {
  var qa = image.select('SummaryQA');
  var mask0 = qa.bitwiseAnd(1 << 0).eq(0);
  var mask1 = qa.bitwiseAnd(1 << 1).eq(0);

  var mask = mask0.and(mask1)

  return image.updateMask(mask);
}

var point = geometry

var col = ee.ImageCollection('MODIS/006/MOD13Q1')
.filterBounds(point)
.map(maskMOD13Q1)

print(col)

Map.setOptions('satellite')
var chart1 = ui.Chart.image.series(col.select('NDVI'), geometry, ee.Reducer.first(), 10)
  .setOptions({title: 'titleGraph',
               lineWidth: 0,
               pointSize: 2})

//print(chart1)

var chart1 = ui.Chart.image.series(col.select('NDVI'), geometry2, ee.Reducer.first(), 10)
  .setOptions({title: 'titleGraph',
               lineWidth: 0,
               pointSize: 2})

//print(chart1)

var chart1 = ui.Chart.image.series(col.select('NDVI'), geometry3, ee.Reducer.first(), 10)
  .setOptions({title: 'titleGraph',
               lineWidth: 0,
               pointSize: 2})
//print(chart1)

var im1 = ee.Image(ba.select('BurnDate').toList(1000).get(156))
print(im1)
Map.addLayer(im1,{min:0,max:1},'im1',true)

var im2 = ee.Image(ba.select('BurnDate').toList(1000).get(224))
print(im2)
Map.addLayer(im2,{min:0,max:1},'im2',true)

Map.addLayer(ba.reduce(ee.Reducer.count()),{min:0,max:1},'count',true)

print(ba)

var colbd1 = ee.ImageCollection('MODIS/006/MOD13Q1')
.filter.('system:time_start','1570043461095894')
.filterBounds(point)
.map(maskMOD13Q1,'try')

print(colbd1)

but I get this error:

  Received too many arguments to function filter(). Expected at most 1
  but got 2.


texts: Using the MOD expression on an SQL Feature Class is not possible?

texts:
I use the MOD expression quite frequently in my definition queries in ESRI to extract fields that are divisible by a specific vale

for example: MOD( [contour], 50) = 0 OR contour = MOD(0, 50)

However I have just moved my feature class from an ESRI file geodatabase onto an ESRI SQL Server and by doing this is appears to make the MOD expression redundant.

The error that comes up when i enter in the normal MOD expression is as follows

  "There was an error with the expression. Underlying DBMS error
  [[Microsoft][SQL Server Native Client 11.0][SQL Server]'MOD' is not a
  recognized built-in function name."

Has anyone experienced this and if so, were you able to find another expression that does the same thing as the MOD expression?


texts: Why ExcelWriter is not working from pands library in QGIS Python console?

texts:
When I am trying to write excel file using python Pandas library in  QGIS(version 3.12 & 3.10), it is through up error saying "ModuleNotFoundError: No module named 'openpyxl'"
How to resolve this issue in QGIS?
When I run my code in QGIS python console, it is showing the error as shown below.
Here my code is
>>> import pandas as pd
>>> writer  = pd.ExcelWriter("output.xlsx")
Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "C:\PROGRA~1\QGIS3~1.12\apps\Python37\lib\site-packages\pandas\io\excel\_openpyxl.py", line 18, in __init__
    from openpyxl.workbook import Workbook
  File "C:/PROGRA~1/QGIS3~1.12/apps/qgis/./python\qgis\utils.py", line 743, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'openpyxl'


texts: QGIS 2.2 Identify tool highlight settings

texts:
Does anyone know how to change the highlighting settings of Identify tool in QGIS 2.2?
I mean, when i click on a feature with the Identify tool, the feature gets highlighted in red. Can i change it to, lets say, blue?
And where are the highlighting settings settings stored?

W. regards,

Edit- i add a screenshot of the Options|Map Tools window. No Identify tool highlighting there...


texts: Inserting python variables into gdal_calc expression?

texts:
Is it possible to insert a python variable into gdal_calc expression?

For example:

mean = [stats[2]] ##mean value
gdal_calc.py -A Input.tif --outfile=Output.tif --calc="1*(A>=mean)" --NoDataValue=0 --type='Byte'


texts:
I'm trying to create a buffer using OpenLayers with JSTS library.

However, even after transforming geometry from EPSG:4326 to EPSG:3857, the buffer method doesn't create distance in meters as expected.

Take this code snippet, which runs on an OpenLayers Map with view with projection EPSG:3857:

const geoInput = {
    type: 'Point',
    coordinates: [
        14.432940,
        50.077371
    ]
};
const geoJsonFormat = new GeoJSON();
const features = geoJsonFormat.readFeatures(geoInput, {
    dataProjection: 'EPSG:4326',
    featureProjection: 'EPSG:4326'
});
const parser = new jsts.io.OL3Parser();
parser.inject(Point, LineString, LinearRing, Polygon, MultiPoint, MultiLineString, MultiPolygon, Circle);
const originalFeatures = [];
const bufferedFeatures = [];

for(let i = 0; i < features.length; i++) {
    const featGeom = features[i].getGeometry().clone();
    featGeom.transform('EPSG:4326', 'EPSG:3857');
    const bufferedJsts = parser.read(featGeom).buffer(251300); // in km ?
    const buffered = parser.write(bufferedJsts);
    originalFeatures.push(new Feature(featGeom));
    bufferedFeatures.push(new Feature(buffered));
}
const source = new VectorSource();
source.addFeatures(originalFeatures);
source.addFeatures(bufferedFeatures);
const vectorLayer = new VectorLayer({source});

This code snippet tries to put a circle from Prague with a distance of 251 km.
The circle obviously doesn't have radius of 251 km, as that's the distance between Prague and Vienna (according to Google Maps) and the buffer doesn't match that.

Even the official example in OpenLayers page doesn't seem like producing correct buffer - https://openlayers.org/en/latest/examples/jsts.html

Am I missing something out, or is it a problem report for JSTS library?


texts: Calculation of mean and standard deviation of raster data in R?

texts:
I have a set of raster images (GeoTIFF, Landsat 1984-2018) which I cropped with my AOI using a shapefile.
Now I want to stack them and calculate the mean and standard deviation of every pixel of the stacked images (in space and time). How can I stack them and calculate the mean and standard deviation? For the outputfile I need one image for the mean and one for the standard deviation.
library(raster)
library(rgdal)
#load data
setwd("C:/Users/cathe/Documents/GEOTIFF_offsetcorrected") dbase = "C:/Users/cathe/Documents/GEOTIFF_offsetcorrected"
#polygon with crop-extend
shape_data <- readOGR("C:/Users/cathe/OneDrive/Documents/ArcGIS/Projects/Cologne/WGS_1984_UTM_Zone_32N/Cologne.shp", stringsAsFactors=FALSE)
#load tif files
ALL_FILES <- list.files(path = dbase, pattern = ".tif$|.TIF$")
#output
outfiles = file.path("C:/Users/cathe/Documents/Output", paste0(basename(tools::file_path_sans_ext(ALL_FILES)), ".tif")) #crop
for (i in seq_along(ALL_FILES)) { r = crop(stack(ALL_FILES[i]), shape_data) writeRaster(r, filename=outfiles[i], bylayer=FALSE, format="GTiff", datatype= "INT1U", options="COMPRESS=ZIP", overwrite=TRUE) }
mean <- mean(r)
plot(mean)
writeRaster(x = mean, filename = "mean1.tif", driver = "GeoTiff", overwrite=TRUE)
I get a result, but i am not sure if it is correct?


texts:
I have preexisting caches of maps that I would like to turn into a WMS. I am trying to implement this through Geoserver/Geowebcache. I tweaked the web.xml file to set the Geowebcache_Cache_Dir to the directory with the preexisting cache. I then put the appropriate information into geowebcache.xml and Geoserver can see the cache under Tile Layers, from here I can open it in OpenLayers, however, I can't ping the cache through a WMS call from arcgis or google earth. I checked the box for direct integration with Geoserver WMS but I still cannot seem to get external softwares to see the WMS. I checked fiddler and it seems to be hitting the tiles in openlayers through /geoserver/gwc/service/wms?blahblahblah. If I try to hit that WMS endpoint from anything else I get the ever so helpful error 400.

Does anyone know the proper method for building a WMS from a cache using this software?


texts:
I am using QGIS 2.12 Lyon on a Mac running OS X 10.10.5 Yosemite.

I load a topographic map raster layer which happens to be in MrSID format (36105B5C.sid is the file). Every time, this makes QGIS crash.

I have followed instructions here, but most of it is for Linux or Windows. I have the GDAL 1.11 from KyngChaos and also the MrSID plugin, and have installed Lizardtech's Unified SDK 9.5 as instructed. (Except I installed Unified SDK 10, I think, 'cause 9.5 was not available.)

Still I cannot load the .sid file because it crashes QGIS.


texts:
I want to install QGIS on an iMac OSX 10.7.5. According to the KyngChaos archive, I have to use QGIS 2.12.0-1, with GDAL_complete-1.11EB, matplotlib and Numpy. This is what I am trying to do, unsuccessfully.

I did installation in the following order: I installed GDAL 1.11EB. This put GDAL, GEOS, SQLite, PROJ frameworks in Frameworks folder in Library. Then installed Numpy (I have tried the one in the GDAL1.11EB package and 1.9.2-1). Then installed  matplotlib (I have tried both 1.3 and 1.4.3). Had no trouble with any of these and both the Frameworks and Python 2.7 directories were populated correctly.

However, after this, when I tried to finish installing QGIS, QGIS 2.12.0-1 fails to install unless I install SQLite framework 3.8, for which it needs PROJ framework 4.9. QGIS 2.12.0-1 installs only if I install these on top of the earlier GDAL 1.11EB. But then it does not launch. Simply 'quits unexpectedly'.

What I understand from this process is that, there is a mismatch between what QGIS 2.12.0-1 is expecting to see regarding the SQLite and PROJ frameworks, and what is contained in GDAL-1.11EB. Simply installing the frameworks it wants to see is obviously not enough to get it to work.

How do I resolve this?

I hope this clarifies my problem. It really is only about how to install the basic QGIS 2.12.0-1 on a Mac OS 10.7.5 successfully.


texts: Are layer files compatible between ArcGIS 10 and 9.3?

texts: If I save layer files in ArcGIS 10, can they be opened and modified in ArcGIS 9.3?


texts: ArcGIS Server Object Extensions - Visual Studio Templates

texts:
The  Visual Studio Templates for ArcGIS Server Object Extensions show up correctly in VS 2008, but not in VS2010.  Is there an easy way to resolve this?  Outside of copy and pasting code from VS2008 to VS2010.

Thanks,
Seth


texts: Making a custom US forest topo map with R

texts:
I am planning a hike, and of course by luck the trail I want to do is spread over 4 US forest topo sheets. So I am trying to make a single geotiff using the raster package in R.

# each file is about 15MB.

f1 <- 'http://data.fs.usda.gov/geodata/rastergateway/downloadMap.php?mapID=472212115&mapType=tif'
f2 <- 'http://data.fs.usda.gov/geodata/rastergateway/downloadMap.php?mapID=472212122&mapType=tif'
f3 <- 'http://data.fs.usda.gov/geodata/rastergateway/downloadMap.php?mapID=473012122&mapType=tif'
f4 <- 'http://data.fs.usda.gov/geodata/rastergateway/downloadMap.php?mapID=473012115&mapType=tif'
download.file(f1,destfile = '~/chikamin.tif',method='wget')
download.file(f2,destfile = '~/sno.tif',method='wget')
download.file(f3,destfile = '~/sno_l.tif',method='wget')
download.file(f4,destfile = '~/bsm.tif',method='wget')
require(raster)
require(spatial)
r1=raster('~/chikamin.tif')
r2=raster('~/sno.tif')
r3=raster('~/sno_l.tif')
r4=raster('~/bsm.tif')
r1
tiled=extend(extend(extend(r1,r2),r3),r4)
tiled
plot(tiled)

but in my case the tiled map is only showing one of the original maps, and colors seems off. When I try the mosaic, it complains about the origins being different.

How can I accomplish this?


texts: How to load dxf2Shp Converter plugin from the in build python console in QGIS 2.6?

texts:
I'm trying to write a script to automate loading of dxf2Shp Converter using the build-in python console. The script should open the pop up of that plug-in. I've tried using .available_plugins/.loadPlugin() but it's not working for me.

Any suggestions?


texts: How to sum pixels of two overlapping rasters with different extents

texts:
I have two DEM rasters which overlap and have different extents. I would like to create a new raster which sums the pixels where there is an overlap.
I tried using gdal_calc first by doing A+B but seeing the rasters are different dimensions it errors.
I then tried using gdalwarp using the -r sum. This creates a raster but it has not preformed a sum.
I have successfully performed this in QGIS using the raster calculator, but as I will eventually be running it on multiple files I'm trying to get a command line solution to work.
The closest solution I've found on Stack Exchange is this question Using rasters with different dimensions in GDAL_calc.py
Is there another way?


texts: ArcGIS Server geoprocessing service examples for shapefile data export

texts:
I have a web application client built using ArcGIS Server Manager which is currently showing a single map composed by three data layers. The map data are provided as a map service by my ArcGIS Server 10.1 instance.

My task is to provide the map's users a way to download the full map locally as a shapefile (e.g: zipped) via her/his web browser.

I was hinted to read here and to build a geoprocessing service to accomplish the task. I decided to write a specific toolbox myself using Python, then publish it as a geoprocessing service on my ArcGIS Server instance and finally include it into my web application.

Can you point me towards any third-party (meaning: non-ESRI) code sample to start me up?

I can find no examples other than the ones in ESRI's support website!


texts: Reading through shapefile table with Geotools

texts:
Goal: Read through a shapefile table with Geotools.

I'm confused about FeatureSource, FeatureCollection, and FeatureIterator in Geotools.
In ArcGIS, I can open a table in an existing shapefile, and use a cursor to read through it, row by row. Here is the successful code using arcobjects jar & ESRI syntax:

IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactory();
IFeatureWorkspace featureWorkspace = (IFeatureWorkspace) workspaceFactory.openFromFile(fileDirectory, 0);
// field not found = -1
ITable table = new Table(featureWorkspace.openFeatureClass(fileName));
ICursor cursor = table.ITable_search(null, false);
IRow row = null;
int indexID = table.findField("ID");

while ((row = cursor.nextRow()) != null) {
    int idInt = (int) row.getValue(indexID);
}
// then do things with that data from the table

But in Geotools, I'm struggling to figure out, what is the equivalent? The FeatureIterator, maybe?

Summation: In order to read through a shapefile table and get data, what is the Geotools syntax?


texts: How to combine lines from different sources reperesenting the same road network but with slightly different geometry?

texts:
I have got two layers of polylines which represent the same road network but they are from different sources and have different geometry.
The two layers are mostly the same tracks with larger or smaller offsets and a segmentation that is very different. But each layer contains objects that do not appear in the other. Concerning the attributes, they are different and both contain interesting information that would be desirable to retain the maximum.

I'd like to combine these two layers to have a complete network, which includes the maximum of information contained in the source layers attributes.

I have not been able to get satisfactory results in ArcGIS10 nor in Safe's FME.


texts:
I have a layer in GeoServer called "pois". The layer is in EPSG:4326. Also, I have another layer but it is in: EPSG:31256

The problem that I'm having is that I'm trying to reach an all points around coordinates, but when I run query a get dozen of results but not the one point I need.

Here is how my query looks like I'm testing it trough postman:

http://localhost:8080/geoserver/wfs?request=GetFeature&version=2.0.0&srsName=EPSG:4326&typeName=omv_workspace:osm_pois&count=100&cql_filter=dwithin(geom,Point(14.81850 48.77621),500,meters)

My query executes but returns bunch of points but not the one I need the closest one to Point(14.81850 48.77621), that one is shown only after I remove count.

Also I passed that point from OpenLayers view in GeoServer but point is not returned.


texts: Why is satellite imagery data delivered as signed 16 bit integers?

texts:
Landsat 5 tm reflectance product images are delivered as signed 16 bit images, why is this preferred to unsigned 8 bit that the DN's are stored as level 1 products for the same satellite? Look at table 6-1 on page 19 of the following document:

https://landsat.usgs.gov/sites/default/files/documents/ledaps_product_guide.pdf.


texts: Issue on Getting Geometry of pgRouting Query in PostGIS

texts:
I am running this query against a Table victoria enabled with pgRouting and PostGIS

SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                SELECT gid AS id,
                source::integer,
                target::integer,
                length::double precision AS cost
                FROM victoria',
9411, 3986, false, false);

which in result return 269 rows in a table like

I also need to get geometry multi string of nodes so I tried to achieve them like

SELECT route.*, ST_AsText(geom) FROM victoria
JOIN(
SELECT seq, id1 AS node, id2 AS edge, cost FROM pgr_dijkstra('
                         SELECT gid AS id,
                         source::integer,
                         target::integer,
                         length::double precision AS cost
                         FROM victoria',
                         9411, 3986, false, false)) AS route
ON
victoria.gid = route.seq;

this returns 268! rows plus the exactly same result of geom for all queries no matter if I change the 9411, 3986 any other number the geometry is always same. can you please let me know how I can get the correct geometry

Update


texts: Normalizing or parsing feature info from GeoServer?

texts:
I'm building a website to visualize layers from services on an OpenLayers Map.

The services come from different geoservers some of them which I don't own.

All services are WMS and don't have WFS.

I don't know what attributes exist on the layers.

I can only make getfeatureinfo request with text/xml or text/html formats

With this in mind, I'm trying to create a method to show the features info on click but different geoservers send me different results is there a way to ask GeoServer for a default template html?

Or use the readFeatures() from OpenLayers to read the XML returned from a getfeatureinfo request?

I can solve this with C# or JavaScript as I'm able to do some server side coding


texts: MMQGIS Hub Distance Error - 'NoneType' object has no attribute 'boundingBox'

texts:
I'm trying to use the MMQGIS Distance to Nearest Hub tool but keep getting this error.

I'm on Windows 7, QGIS 2.18.2 and MMQGIS 2017.2.28. Any ideas?

An error has occurred while executing Python code:

AttributeError: 'NoneType' object has no attribute 'boundingBox'
Traceback (most recent call last):
File "C:/Users/EJO31/.qgis2/python/plugins\mmqgis\mmqgis_dialogs.py", line 1504, in run
   nameattributename, units, addlines, savename, evenly_distributed, 1)
File "C:/Users/EJO31/.qgis2/python/plugins\mmqgis\mmqgis_library.py", line 3816, in mmqgis_hub_distance
 wgs84_point = feature.geometry().boundingBox().center()
AttributeError: 'NoneType' object has no attribute 'boundingBox'

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.2 Las Palmas, 102ee87

Python Path:
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python
C:/Users/EJO31/.qgis2/python
C:/Users/EJO31/.qgis2/python/plugins
C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\nose-1.3.3-  py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\bin\python27.zip
C:\PROGRA~1\QGIS2~1.18\apps\Python27\DLLs
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\plat-win
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGIS2~1.18\bin
C:\PROGRA~1\QGIS2~1.18\apps\Python27
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/EJO31/.qgis2//python
C:\Users\EJO31\.qgis2\python\plugins\mmqgis/forms
C:/Users/EJO31/.qgis2/python/plugins\NNJoin
E:/Fixed Broadband Model/pcd_2_cab_2_exchange_data


texts: Force openlayers to not use browser cache for tiles refresh

texts:
I use OpenLayers.Layer.XYZ to display tiles from TileStache server without using the cache option of the server. However, I notice that the tiles stay in cache (probably the browser cache) until the whole page is refreshed with Ctrl-F5.

If I want to redraw the XYZ layer, it does not works because the map use the browser cache. Is there a way to force the map to not use that cache? So could refresh the layer by asking fresh tiles from the server.

var map = new OpenLayers.Map('map',
          { projection: new OpenLayers.Projection("EPSG:3857"),
            numZoomLevels: 20 });

var tiledLayer = new OpenLayers.Layer.XYZ('TMS',
                "{{ tmsURL }}1.0/layer/{{ shapefile.id }}/${z}/${x}/${y}.png?"+ new Date().getTime());

map.addLayer(tiledLayer);
tiledLayer.redraw(true); //does not make new GET requests to the server


texts: I have a point shapefile containing elevation values in the id field. The shapefile is overlaying a raster layer of the topographical slope. Is it possible to multiply the elevation value of each point by the respective slope value of the raster cell laying immediately below the point?


texts: Export ArcGIS Model to Python - Iterate Rasters, Collect Values, Calculate Values

texts:
I want to batch process a large number of raster files, calculate some values and properties and export these to a .csv or .xls. I managed to get a good working model that does this perfectly.

However, the next step is to export this to Python and with this I encounter some problems.
I am not an expert in Python code. I also want to try this to get a little bit acquainted with it. The most important aim of exporting it is that the model should be compatible to use in different working environments.

Exporting the code is difficult because I used some tools that are not exportable to Python: Iterate Rasters, Collect Values and Calculate Values. Therefore I need to make a work-around.

What is the best way to accomplish that?

Here is the code posted on a GIS SE question.

I hope this is clear enough. Otherwise just ask for clarification.


texts: MMQGIS plugin causes QGIS crash while geocoding

texts: When trying to geocode data with MMQGIS pulgin (Google Maps) the QGIS crashes (not responding). I even tried to divide data into smaller sets but it does not help.


texts: how to get concave surface which close to the outer edge point

texts: In qgis, how to generate the blue outer edge (the concave surface close to the outer edge point) in Fig. 2 from the point in Fig. 1. Any suggestion?


texts: No margins displayed in Layout view of ArcMap

texts:
I have an issue with Layout View in ArcMap.  I need the layout view to look like this (note the white margins around the map):

At the moment, it looks like this:

My version of ArcGIS Desktop is 10.7.1.


texts: How to calculate the basin length - which point to which point and how in ArcGIS 9.3?

texts: This is the area of my research. I have calculated area using area and perimeter using ArcGIS 9.3 for doing morphometric analysis of the area.but I don't know how to calculate the length of the basin that is which start and end point I have to take to calculate using ArcGIS 9.3.


texts: Overlapping legend: How to put geopandas legend next to the map?

texts:
I am using geopandas and contextily to plot building data (polygon) from OpenStreetMap. The problem is, the legend is overlapping the map. Is there any way to put the legend next to the map? or plot the map without legend, but export the legend separately from the map?

I hope the image helps.

here's my code to produce the map (I know the code is a little different, but it yields the same map):

from matplotlib import pyplot as plt
import geopandas as gpd
import contextily as ctx
%matplotlib inline

buildings = gpd.read_file("buildings.shp")
fig, ax = plt.subplots(1,1,figsize=(9,9))
buildings.to_crs(epsg=3857).plot(column='building', ax=ax, cmap="Paired", legend=True)
ctx.add_basemap(ax, url=ctx.providers.Esri.WorldTopoMap)
ax.axis('off')
plt.show()


texts: Why does arcpy.ImportToolbox() not work in ArcGIS 10 (SP5) when calling custom model for topology?

texts:
I have an arcpy/python script using ArcGIS10 and using a custom toolbox and model. The toolbox has an alias of defaultTopo and my model within the toolbox is called AddRuleToTopologyModel. This model was created because there is a known bug when calling the AddRuleToTopology_management tool for the Must Be Covered by Feature Class of (Area-Area) rule from python and the work around is to use this model. I'm simply updating my script from 9.3 to 10.0.

The code goes like this:

# Defining location of toolbox which is in the same location as the script
toolboxLocation = sys.argv[0][:sys.argv[0].rfind('\\')]

# Importing custom toolbox
arcpy.ImportToolbox(toolboxLocation+"\ApplyDefaultTopology.tbx")

# Calling custom model for polygon feature class - not all the variables are posted here
arcpy.AddRuleToTopologyModel_defaultTopo(featureDataset+"\\"+topologyName,'Must Be Covered By Feature Class Of (Area-Area)',featureClass,featureDataset+"\\"+boundaryFC)

Now if i run the script from Eclipse it does not pick up my featureDataset location I use as argv1 in the script throwing the error I have in place.

If I try and run the script from ArcCatalog I get the following error:

<type 'exceptions.AttributeError'>: 'module' object has no attribute 'exists'

Eclipse also tells me that arcpy.AddRuleToTopologyModel_defaultTopo() is an undefined variable when trying to only use arcpy.

The strange thing is if i use gp.arcgisscripting(9.3) with gp.AddToolbox() to replace arcpy.ImportToolbox() the script runs fine from Eclipse and ArcCatalog because arcpy still supports the 9.3 geoprocessing object. I've been asked to remove all references to the arcgisscripting module and make the new script pure arcpy to function on our new servers.

Thoughts?

Thanks in advance.


texts: Using pyshp to create line shapefiles with color

texts:
I am using PyShp to create line shapefile in Python.

w = shapefile.Writer('shapefiles/line')
w.field('name', 'C')
w.line([ [[3,2],[2,6]] ])
w.record('linestring1')
w.close()

Is there a way to specify a different color (and width) for each line object?


texts: Update extent in MultiPolygon postgis layer layer do not work

texts:
Wen trying to automatically update extent, nothing happens and extent still covers whole planet.

Name
...
Source
dbname='postgres' host=localhost port=5432 user='postgres' sslmode=disable key='"id"' srid=4326 type=MultiPolygon checkPrimaryKeyUnicity='1' table="public"."..." (geom) sql=
Storage
PostgreSQL database with PostGIS extension
Comment

Encoding
UTF-8
Geometry
Polygon (MultiPolygon)
CRS
EPSG:4326 - WGS 84 - Geographic
Extent
-180.0000000000000000,-90.0000000000000000 : 180.0000000000000000,90.0000000000000000


texts: How to use Button with Behavior in Silverlight Legend?

texts:
I'm trying to use a button with behavior in the MapLayerTemplate of an Esri Silverlight SDK 2.1 legend.

When I click the button and then drag on the map, the map pans instead of showing a selection rectangle.

The behavior code is pulled from this sample, while the legend is from this sample.

If I put the same xaml code within the layoutroot, it works as expected - I see a selection rectangle and the selection completes successfully.

<esri:Legend x:Name="MyLegend" Map="{Binding ElementName=MyMap}"
            LayerIDs="ProjectPolys,AssetPoints,AssetLines,BaseLayers,WaterInfrastructure,SewerInfrastructure"
            LayerItemsMode="Tree"
            ShowOnlyVisibleLayers="False"
            Refreshed="Legend_Refreshed">
    <esri:Legend.MapLayerTemplate>
        <DataTemplate>
            <StackPanel Orientation="Horizontal">
            <CheckBox Content="{Binding Label}"
                    IsChecked="{Binding IsEnabled, Mode=TwoWay}"
                    IsEnabled="{Binding IsInScaleRange}">
                </CheckBox>
                <Button Content="Select">
                    <i:Interaction.Triggers>
                        <i:EventTrigger EventName="Click">
                            <esriBehaviors:SpatialQueryAction
                                DrawMode="Rectangle"
                                LayerID="MyGraphicsLayer"
                                Url="http://hq-gis-01/ArcGIS/rest/services/WaterInfrastructure2/MapServer/15"
                                Symbol="{StaticResource GraphicsLayerFillSymbol}"
                                TargetName="MyMap" />
                        </i:EventTrigger>
                    </i:Interaction.Triggers>
                </Button>
                <Slider Maximum="1" Value="{Binding Layer.Opacity, Mode=TwoWay}" Width="50" />
            </StackPanel>
        </DataTemplate>
    </esri:Legend.MapLayerTemplate>
    <esri:Legend.LayerTemplate>
        <DataTemplate>
            <CheckBox Content="{Binding Label}"
            IsChecked="{Binding IsEnabled, Mode=TwoWay}"
            IsEnabled="{Binding IsInScaleRange}" >
            </CheckBox>
        </DataTemplate>
    </esri:Legend.LayerTemplate>
</esri:Legend>

Update:
Following Seth's suggestion I tried wiring the trigger in the code-behind.  I get the same results.  The code behind approach fails for the button in the datatemplate, but works for the one on a stackpanel sibling to the map.

private void Assign_Behavior_Click(object sender, RoutedEventArgs e)
{
    var button = sender as Button;
    ESRI.ArcGIS.Client.Actions.SpatialQueryAction sqa = new ESRI.ArcGIS.Client.Actions.SpatialQueryAction();
    sqa.DrawMode = DrawMode.Rectangle;
    sqa.LayerID = "MyGraphicsLayer";
    sqa.Url = @"http://hq-gis-01/ArcGIS/rest/services/WaterInfrastructure2/MapServer/15";
    sqa.Symbol = LayoutRoot.Resources["GraphicsLayerFillSymbol"] as Symbol;
    sqa.TargetName = "MyMap";

    var et = new System.Windows.Interactivity.EventTrigger("Click");
    et.Actions.Add(sqa);
    var triggerColl = System.Windows.Interactivity.Interaction.GetTriggers(button);
    triggerColl.Add(et);
}

I'm beginning to think this is not an Esri issue - I can't find an example anywhere showing how to attach trigger actions from buttons in a datatemplate.


texts:
I am trying to clip a raster using a geoJSON without writing to my disk.
The following code compiles and produces the correct result.
import tempfile
import fiona
import rasterio
from rasterio.io import MemoryFile
import rasterio.mask
from scipy import stats

def getSummStats(raster, shapefile):
    #Open Polygon
    with fiona.open(shapefile, "r") as shapefile:
        shapes = [feature["geometry"] for feature in shapefile]

    #Clip the raster
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta

    #Setup to save and write the clipped raster in memory
    with MemoryFile() as memfile:
        with memfile.open(driver = "GTiff", height = out_image.shape[1], width = out_image.shape[2], count=1,
                          dtype=rasterio.float64, transform = out_transform) as dataset:
            #write the clipped raster to memory
            dataset.write(out_image)
            #read the clipped raster from memory
            ds = dataset.read()
            return stats.describe(ds, axis = None)

with rasterio.open("B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2") as src:
    print(getSummStats(src, "C:/Users/Tom/Downloads/canada_provinces/canada_provinces/canada_provinces.geojson"))

However, I would like to be able to call this script from the command line and input the raster and geoJSON when the script is called. I made some changes to accomplish this.
import tempfile
import fiona
import rasterio
from rasterio.io import MemoryFile
import rasterio.mask
from scipy import stats
import os, sys

def getSummStats(raster, geoJSON_path):
    #Open Polygon
    with fiona.open(geoJSON_path, "r") as shapefile:
        shapes = [feature["geometry"] for feature in shapefile]

    #Clip the raster
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta

    #Setup to save and write the clipped raster in memory
    with MemoryFile() as memfile:
        with memfile.open(driver = "GTiff", height = out_image.shape[1], width = out_image.shape[2], count=1,
                          dtype=rasterio.float64, transform = out_transform) as dataset:
            #write the clipped raster to memory
            dataset.write(out_image)
            #read the clipped raster from memory
            ds = dataset.read()
            return stats.describe(ds, axis = None)

def main(raster_path, geoJSON_path):
    with rasterio.open(raster_path) as src:
        print(getSummStats(src, geoJSON_path))

if __name__ == '__main__':
    #
    # example run : $ Polygon_CMDVersion.py /<full-path>/<raster-name>.grib2 /<full-path>/<geoJSON-name>.geoJSON
    #
    if len( sys.argv ) < 2:
        print ("[ ERROR ] you must two args. 1) the full raster path and 2) the full geoJSON path")
        sys.exit( 1 )

    main(sys.argv[1], sys.argv[2])

However, this second method gives me the following error
File "C:\Users\Tom\Downloads\Polygon_CMDVersion.py" line 21, in getSummStats
out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
NameError: name 'src' is not defined

If this helps, here is what I wrote in the command line
Polygon_CMDVersion.py B:/dd.weather.gc.ca/model_gem_global/15km/grib2/lat_lon/00/000/CMC_glb_TMP_TGL_2_latlon.15x.15_2020052500_P000_Copy.grib2 C:/Users/Tom/Downloads/canada_provinces/canada_provinces/canada_provinces.geoJSON

How can I resolve this issue?


texts: I want to manually scale selected polygon feature in QGIS so I don't need to set any parameters, just click on the feature and scale it moving the mouse. In ArcGIS there is a Scale tool. I want the same thing for QGIS. Note that I don't need v.transform and Affine plugin. As I said I want to scale it manually.


texts: How to Calculate NDWI Value in Arcgis 10.1?

texts:
I want to delineate surface water feature of my Landsat raster image and also want to calculate the area of water body. I have learned through internet that NDWI can be used to delineate open water feature.

How can I calculate NDWI value using arcgis 10.1 ?


texts:
I am kind of a beginner. I have exported my webmap using QGIS2Web and while playing with it, I was able to change the settings of the items already present on the map (for example position, language, zoom levels, etc.) Although, I would like to add some functions to it.

I want to add the Leaflet.Dialog plugin (https://github.com/NBTSolutions/Leaflet.Dialog), but I'm failing somewhere.
In my export folder, I added both the .css and .js to the project and imported it with a link and script tag.

I then added this to the index.html :

var dialog = L.control.dialog(options)
          .setContent("
Hello! Welcome to your nice new dialog box!
")
          .addTo(map);

Why does this code not show anything?

Here is the code to link the plugin :

<link rel="stylesheet" href="css/Leaflet.Dialog.css">


texts: Creating a local coordinate system in QGIS

texts:
I'm trying to make a local coordinate system using the PROJ.4 coding. I thought I had it figured out but come up with an error when trying to check parameters. Can anyone tell me what is wrong with the following parameters? I'm using Cortuna 3.10.1. Thanks.
+proj=tmerc +datum=NAD83 +lat_0=40.71232748 +lon_0=-117.2831823 +k=1.000239646 +x_0=50000 +y_0=150000 +ellps=GRS80 +unit=us-ft
+no_defs


texts: OpenLayers2 Vector Event Popup with Textbox

texts:
The OpenLayers map here:
http://www.dassh.ac.uk/demonstrations/valmer/demo5.html

has an event to add a vector layer, plus creates a textbox on the page below.
Ideally I'd like the click event to add the marker, plus display a popup with a textbox/textarea where the user can fill in information on the location.  This would then be submitted to the textbox below.

Can anyone help?


texts:
I have a few indoor map construction algorithms and the ground truth map. I am trying to compare their construction accuracies w.r.t. the ground truth map.
Like any other scientific evaluations, mere visual comparisons are not enough. I need to quantize the errors so that I can claim which one is superior.
I am only interested in the corridors. So basically by construction accuracy, I mean how accurately the corridors (the corridor walls) are constructed. So I am indifferent about the room walls, etc.
I find this not easy, because this problem lie in between image processing and graph theory. On the one hand, I cannot simply compare the two images with image processing techniques, because I am not interested those trivial details mentioned above. On the other hand, the constructed map is not a collection of vertices and edges. So graph theory may not work here, either.
I have thought about overlapping the constructed on over the ground truth and see how well they are overlapped. A problem with this is what if my constructed map is of a slightly different scale with the original one? Then they will overlap very poorly, albeit the constructed map is actually ok.
So is there any good standard way of comparing two indoor maps?

An Example
Ground truth is the background floor map. The red curve is what the algorithm constructs. I wish to compare and quantify the construction errors. (The thickness of the red curve is only for visual purposes. It is actually a deterministic mathematically expressed curve)