texts: We have recently upgraded our ArcSDE 9.3 instance to ArcSDE 10.  I was hoping to upgrade postgresql 8.3 to 8.4 at the same time.  Although ArcSDE 10 supports 8.4, it turns out that the ArcSDE wizards do not.  Everything must be done via the command line.  I made the call to stick with 8.3 just to keep things easier.  Am I missing out?  Does anyone have any experiences with 8.3 vs 8.4 in ArcSDE, if so, what are your thoughts?


texts: How to combine multiple overlapping rasters and outputs a new raster that carries the maximum value

texts:
I have 7 classification rasters.  Each raster has probability values from 1 to 10. All of the rasters overlap each other at some areas.  I would want to create a new raster that carries the maximum probability value on the overlap of the 7 rasters, and at the same time keeps the classification of the original raster.  For example:

In one location:
raster 1 (max is 10) overlaps with raster 2 (max is 9) and raster 3 (max is 5) ----> expected output of the pixel should have the value of raster 1

In another location:
raster 1 (max is 2) overlaps with raster 2 (max is 9) and raster 3 (max is 4) ----> expected output of the pixel should have the value of raster 2

I tried "Mosaic To New Raster" tool with "maximum" but the output is simply a map showing the maximum value of the pixels, without me knowing from which raster it was coming from.

So it is possible to use ArcGIS to create this new raster (without the overlap because the overlap has already been assigned to the classification raster that carries the maximum value)?

I have been looking at existing questions, but I am not finding one that answers my question.  Or if there is any, can someone point me to it?

I am using ArcGIS 10.4.

From comments by @FelixIP, this could be solved using Spatial Analyst Tools --> Local --> Highest Position!


texts: Problem reprojecting annotation in ArcGIS 10

texts:
I have a large annotation feature class that I want to reproject from decimal degrees to a specific iteration of Albers.  Using the Project tool, it appears to work, all annotation features are present in the output, attributes look fine, but the features don't draw.  I can select them in the data frame, but can't see them.

What to do?


texts: How to rotate text label anti-clockwise in QGIS?

texts: i'm relatively new to qgis. I am using version 2.14.7. I would like to rotate a text label anti-clockwise in map composer but the Rotation settings in Item Properties only allow for clockwise rotation. What do I do?


texts: How to efficiently export specific files from .hdf to .img

texts: I am working with MODIS .hdf grid files and am looking for some help on exporting specific files from the .hdf. There are 200+ .hdf's and I need to extract the 0, 1, 2, and 11 (NDVI, EVI, quality, and reliability) layers, ideally to .img's and with the nomenclature either intact or editable. Does anyone have a quick way to pull this off, either in ArcGIS, Grass, or QGIS?


texts:
Is the below query correct to get a convex hull shape:

DATA "the_geom from (select ST_ConvexHull(ST_Collect(the_geom)) As the_geom FROM xyz) as foo using unique id using SRID=XXXX"


texts: Same field names as given in text file (in a script)

texts:
I have a raster (ras) and a text file (T, DA, etc). When I join both then field names are like following:

ras.vat:T
ras.vat:DA

I have to access these fields using SearchCursor for further processing.
Is there a way to avoid these long combined names and simply obtain "T" and "DA" after joining?


texts: Geoserver WFS GetFeature doesn't work - "SyntaxError: missing ; before statement"

texts:
I'm trying to use GeoServer WFS JSON for creating a vector layer (instead of using WMS) because this is a layer I redraw every 20 seconds.

This is the code I use :

function getJSON() {

    var geoJsonUrl = "http://localhost:8080/geoserver/test/ows?service=WFS&version=1.1.0&request=GetFeature&typeName=test:Coords&srsName=EPSG:4326&format=text/javascript&outputFormat=json&format_options=callback:loadGeoJson";

    $j.ajax({
        url: geoJsonUrl,
        dataType: 'jsonp'
    });
}

function loadGeoJson(data) {
    console.log('success');
}

The JSON from the request is :

{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "id": "DIView.fid-25487b92_142a9c99a7e_-5722",
            "geometry": {
                "type": "Point",
                "coordinates": [
                    34.8923,
                    32.1893
                ]
            },
            "geometry_name": "DeviceLocation",
            "properties": {
                "Id": 1,
                "DeviceSerial": "046009117212",
                "DateTime": "2013-12-03T12:20:48.080Z",
                "Spd": "4.72175   ",
                "Brn": "143.3     ",
                "Alt": "91.6      ",
                "DeviceName": "12",
                "DeviceReversedName": "21",
                "DeviceBattery": 3.71
            }
        },
        {
            "type": "Feature",
            "id": "DIView.fid-25487b92_142a9c99a7e_-5721",
            "geometry": {
                "type": "Point",
                "coordinates": [
                    34.8025,
                    31.2813
                ]
            },
            "geometry_name": "DeviceLocation",
            "properties": {
                "Id": 2,
                "DeviceSerial": "046009117224",
                "DateTime": "2013-12-03T12:41:12.313Z",
                "Spd": "0         ",
                "Brn": "136       ",
                "Alt": "349       ",
                "DeviceName": "24",
                "DeviceReversedName": "42",
                "DeviceBattery": 0.29
            }
        },
        {
            "type": "Feature",
            "id": "DIView.fid-25487b92_142a9c99a7e_-5720",
            "geometry": {
                "type": "Point",
                "coordinates": [
                    34.8026,
                    31.2812
                ]
            },
            "geometry_name": "DeviceLocation",
            "properties": {
                "Id": 3,
                "DeviceSerial": "047180050069",
                "DateTime": "2013-12-03T12:41:14.140Z",
                "Spd": "0.690988  ",
                "Brn": "161.8     ",
                "Alt": "340.2     ",
                "DeviceName": "4",
                "DeviceReversedName": "4",
                "DeviceBattery": 3.65
            }
        }
    ],
    "crs": {
        "type": "EPSG",
        "properties": {
            "code": "4326"
        }
    }
}

and I get error in Firebug :

Any idea ?!


texts: How to symbolize based on feature class and assign a unique value per point?

texts:
I am trying to produce a map in ArcGIS 10.1 that symbolizes by feature class but also assigns a unique number (as part of the symbol) to each location that will relate to a table embedded in the map that has more detailed descriptions for each location (i.e. location name, address, etc.).  My current map (below) has unique symbols for monitoring wells, endangered species, churches, schools, daycares, etc., but it doesn't leave any place for a unique value per location (starting at 1, then 2, 3, etc.).

Would it be possible to create a custom symbol with a unique color (e.g. blue in-filled circle for hotels) for each feature class and then have it automatically assign a consecutive, unique number in the center for each location?  An easy out would be to keep the symbology I have and just have number labels for each location.  I was just hoping to avoid having to place labels at each location.  If the solution requires heavy scripting, then I may not have the skills to do it.  Thanks.

Example symbol I would like:

Example table key I would like embedded in the map somewhere:


texts: Calculate radius from magnitude of earthquake on leaflet map?

texts:
I have a working leaflet map that is showing earthquake data from http://earthquake.usgs.gov/. I'm trying to figure out the "right" way to calculate the radius of the circle for an earthquake given its magnitude (and possibly depth).

The circle is being created with

L.circle(latlng, radius);

(where radius is in meters). The previous map developer did this:

var radius = 30000 * magnitude;

Which looks good on a map but doesn't seem right - he just made up a number to make it look decent. I googled and found two other earthquake maps using leaflet and they did this:

var radius = Math.sqrt(Math.pow(10, magnitude))/50000);

and the other did this:

var radius = Math.pow(2, magnitude);

but both of those draw very tiny circles on the map. As an example, if the magnitude is 5.0, I would get radii of 150000, 1.414, and 32, respectively.

So, my question... is there a standard/accepted formula for drawing an earthquake circle on the map using the magnitude (and maybe the depth)?

UPDATE: I found one more:

var radius = (Math.exp(magnitude/1.01-0.13))*1000;

which for a magnitude of 5.0 gives a radius of 124026 meters (and "looks" good on the map). Barring any insights from the community, I may just go with this one.


texts: Homogeneous vs Heterogeneous schema in Postgis

texts:
Introduction

I'm working on creating a spatial database in Postgis, to be viewed and edited trough QGIS to help managing some farmland areas.

One of my goals is the have a centralized place where all actions to develop in farm are listed. The actions all have (kind of) the same structure (attributes) and would be stored in a non spatial table. Each action will be related to a particular feature in the farm (a road, a sign pole, a gate, a fence...).

The problem

Since features will have to be represented in different geometry types, (at least) two options arise.

The homogeneous schema

Each geometry type is stored in a different table (one for points, one for lines and one for polygons). This makes opening and editing the tables in QGIS (and other GIS clients) much easier as qgis does not accept tables with multiple geometries. In the other hand I can't seem to find a way to consistently manage a unique key in the three tables to use as foreign key in the actions non-spatial table.

the heterogeneous schema

All features are stored in the same table, where the geometry column is of the type geometry(geometry, 3763). This will make very easy to manage a unique key to relate to the actions non-spatial table. On the other hand, viewing and editing the table trough QGIS will require the creation of VIEWS to separate the different geometry types, and INSTEAD TRIGGERS to redirect edits on views into the original table.

The question

Are this my only options? Is it possible to use the homogeneous schema for this case? If yes, how?

Please enlight my database schema creativity


texts: Explaining gap between Robinson WGS84 (ESRI:54030) in QGIS and ArcMap?

texts:
Does anybody know why there a gap appearing between two point layers create with ESRI and QGIS when I a reproject, from the same source file, to Robinson?

The source coordinate system and destination coordinate system share the same datum. It seems there's a shift depending on the technology used to reproject.

Here are the steps to reproduce:

I'm using a shapefile point layer in Lat/Long WGS84 Coordinate System and I  reproject it to Robinson (ESRI:54030) using the Tool "Project" in ArcMap

I do the same thing in QGIS. I take the same point layer in WGS84 Lat/Long and reproject it to Robinson (ESRI:54030) using the "Save as.."

Download the sample point data here

At the request of mkennedy, here is the cs2cs conversion of regular 10° points:

deg_E   deg_N   sph_X     sph_Y      ell_X       ell_Y
 10       0   943711.34    0        944768.52       0
 10      10   939370.28 1068322.39  940422.59   1069519.16
 10      20   926913.31 2136644.78  927951.66   2139038.32
 10      30   905962.87 3204967.3   906977.76   3208557.61
 10      40   869724.36 4271566.44  870698.65   4276351.58
 10      50   819047.09 5320935.04  819964.61   5326895.73
 10      60   753647.9  6328948.79  754492.16   6336038.68
 10      70   678150.95 7267177.19  678910.63   7275318.12
 10      80   586327.84 8093403.95  586984.66   8102470.44
 10      90   502243.15 8615503.11  502805.77   8625154.47

Scaling the results leads to the factors in https://en.wikipedia.org/wiki/Robinson_projection.


texts:
I have a multi-band GeoTIFF with >10 bands, for example 11. I would have expected these to be shown in an ordered numerical sequence from 1 to 11. However, ArcCatalog 10.3.1 wants to show the string sorted list of bands:

How can I fix this to make ArcCatalog sort bands by their band number?

One manual fix within ArcCatalog is to rename (e.g.) Band_1 to Band_01 by clicking each of these, repeating this until Band_09. But that's a lot of clicking and manual effort! And it adds a .aux.xml sidecar file (yuk!)

My other attempt to fix this is to use GDAL to set each band's Description metadata to Band_01, Band_02, ..., Band_11, which looks like this (via gdalinfo):

Band 1 Block=365x1 Type=Float64, ColorInterp=Gray
  Description = Band_01
Band 2 Block=365x1 Type=Float64, ColorInterp=Undefined
  Description = Band_02
...
Band 11 Block=365x1 Type=Float64, ColorInterp=Undefined
  Description = Band_11

ArcCatalog is now happy to sort every band, except for Band_1.

Furthermore, ArcCatalog allows me to rename any band except for Band_1, so I can't change this name or sorting.


texts: ModelBuilder submodel not showing all parameters

texts:
I created a master model with many submodels in it. I'm having problems with the last submodel "Table Join Shapefile". In this submodel there are three input parameters: "TDM_Join_INRIX_Merge","TDM_Network_Projection","INRIX_Network_Projection".The input of these parameters have already been existed in the master model. "TDM_Join_INRIX_Merge" is the output of submodel "TDM Data Preparation". "INRIX_Network_Projection" is the output of submodel "INRIX Data Preparation".
So I just need to connect the existing output to the parameters of submodel "Table Join Shapefile". However, after adding the submodel "Table Join Shapefile" into the master model, it only shows two input parameters instead of three. I have to make an independent variable from parameter of "TDM_Network_Projection" instead of connecting the existing "TDM_Network_Projection" (output of submodel "TDM Data Preparation") to it. I couldn't figure out why.


texts:
I am having a hard time to figure out what this information refers to. It is stated at the product description of Terra MODIS MOD11A1 product  (Wan et al.,2015). (from https://lpdaac.usgs.gov/products/mod11a1v006/)

(It is also mentioned at the dataset page of Earth Engine Data Catalog at
https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD11A1)

"Above 30 degrees latitude, some pixels may have multiple observations where the criteria for clear-sky are met. When this occurs, the pixel value is a result of the average of all qualifying observations."

From Version 5, temporal averaging was removed from the MOD11A1 products, and they have now LST data by instantaneous observations. I am saying this based on;

  For daily MOD11A1 product, in the MODIS LST Products User's Guide of
  Version 4 products (by Wan Z. (2006) ) there was this information  "
  As latitude increases beyond 30 degrees, the LST value at some grids
  in the MOD11A1 product may be the average value of LSTs retrieved from
  multiple MODIS observations in clear-sky conditions in day or night.
  If LST values from single clear-sky MODIS observations are required
  for all grids, the MOD11_L2 LST product and MODIS geolation product
  (MOD03) should be used instead."

  Which was changed into  "As latitude increases beyond 30 degrees,
  there may be multiple MODIS observations in clear-sky conditions. In
  V5 MOD11A1 product, LST values at all grids are from single clear-sky
  MODIS observations by selecting LSTs in MOD11_L2 files at smaller
  viewing zenith angles or the LSTs at larger zenith angles but their
  values being larger by at least 2K" in the Collection-5 MODIS Land
  Surface Temperature Products Users' Guide (Wan Z., (September 2006).

In this context, I would like to ask if anybody could help with the above information? What type of averaging is being referred to?


texts: How to draw Rectangle and Circle using free hand mouse in openlayers3?

texts:
How to draw Rectangle and Circle using free hand mouse in openlayers3 ?
I know we can draw point, lines, and Polygons. I saw this example:
http://openlayers.org/en/master/examples/draw-features.html
but I need to draw rectangles and circles too.


texts: GRASS: r.to.vect outputs empty layer when creating lines

texts:
I am trying to convert a binary raster to polyline using r.to.vect in QGIS 3.4.

It keeps giving me an empty output, however. This is my input raster:

And I am expecting an output (more or less) like this:

All of my attempts result in an empty polyline layer, however, with this in the logs:

Processing algorithm…
Algorithm 'r.to.vect' starting…
Input parameters:
{ '-b' : False, '-s' : False, '-t' : False, '-v' : True, '-z' : False, 'GRASS_OUTPUT_TYPE_PARAMETER' : 2, 'GRASS_REGION_CELLSIZE_PARAMETER' : 0, 'GRASS_REGION_PARAMETER' : None, 'GRASS_VECTOR_DSCO' : '', 'GRASS_VECTOR_EXPORT_NOCAT' : False, 'GRASS_VECTOR_LCO' : '', 'column' : 'value', 'input' : '/home/isaac/raster_to_line/dummy_img.tiff', 'output' : '/tmp/processing_f785fcf97bd74b71ae07fd43d69b03cc/984b737cf59b428380ad5d6b4a2ec9c3/output.gpkg', 'type' : 0 }

g.proj -c proj4="+proj=longlat +datum=WGS84 +no_defs"
r.external input="/home/isaac/raster_to_line/dummy_img.tiff" band=1 output="rast_5d4331d87fa689" --overwrite -o
g.region n=0.0 s=-5.0 e=5.0 w=0.0 res=1.0
r.to.vect input=rast_5d4331d87fa689 type="line" column="value" -v output=outpute51df7d2d5ef43198c5d492a6ac8df29 --overwrite
v.out.ogr type="line" input="outpute51df7d2d5ef43198c5d492a6ac8df29" output="/tmp/processing_f785fcf97bd74b71ae07fd43d69b03cc/984b737cf59b428380ad5d6b4a2ec9c3/output.gpkg" format="GPKG" --overwrite
Starting GRASS GIS...
Executing </tmp/processing_f785fcf97bd74b71ae07fd43d69b03cc/grassdata/grass_batch_job.sh> ...
Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default
Projection information updated
No projection information available
Over-riding projection check
Reading band 1 of 1...
Link to raster map <rast_5d4331d87fa689> created.
WARNING: Raster is not CELL, '-v' flag ignored, raster values will be written to the table.
Extracting lines...
20..40..60..80..100
120..Building topology for vector map <outpute51df7d2d5ef43198c5d492a6ac8df29@PERMANENT>...
Registering primitives...
r.to.vect complete.
WARNING: No lines found, but requested to be exported. Will skip this feature type.
Exporting 0 features...
WARNING: Output layer is empty, no features written
v.out.ogr complete. 0 features (Line String type) written to <outpute51df7d2d5ef43198c5d492a6ac8df29> (GPKG format).
Execution of </tmp/processing_f785fcf97bd74b71ae07fd43d69b03cc/grassdata/grass_batch_job.sh> finished.
Cleaning up default sqlite database ...
Cleaning up temporary files...
Execution completed in 0.61 seconds
Results:
{'output': <QgsProcessingOutputLayerDefinition {'sink':/tmp/processing_f785fcf97bd74b71ae07fd43d69b03cc/984b737cf59b428380ad5d6b4a2ec9c3/output.gpkg, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
Algorithm 'r.to.vect' finished

I think that this is the culprit:

WARNING: Raster is not CELL, '-v' flag ignored, raster values will be written to the table.

But what does this mean? I have indicated in the tool dialog that it should "Use raster values as categories", which is what the -v flag is supposed to do, but because the raster is "not CELL" it gets ignored. The docs don't seem to say anything about what makes it "CELL".

Can anybody shed some light on this?

This is a screenshot of my inputs:

Part 1:

Part 2:


texts: Does ogr2ogr Convert from PostGIS to File Geodatabase (FGDB_BULK_LOAD YES) work?

texts:
I am trying to export a lot of data from a PostGIS database to a File Geodatabase.

The command I used to use was

ogr2ogr -progress -append -f "FileGDB" -gt 66536 -sql "SELECT * FROM line" C:\Temp\data.gdb  PG:"host=localhost user=postgres password=xxx dbname=xxx" -nln Line -nlt LINESTRING -a_srs "EPSG:4326"

This took 34 minutes.

In the 1.9.2 release of GDAL which I am using OSGeo4W there is a new switch to bulk load data into a File Geodatabase. So I changed my code to be

ogr2ogr --config FGDB_BULK_LOAD YES -progress -append -f "FileGDB" -gt 66536 -sql "SELECT * FROM line" C:\Temp\data.gdb  PG:"host=localhost user=postgres password=xxx dbname=xxx" -nln Line -nlt LINESTRING -a_srs "EPSG:4326"

And this has no performance improvement at all.

I have run it several times with and without the other switches and still nothing seems to improve.

Has anyone got the Bulk Load working?


texts: Value maps (comboboxes) doesn't follow when exporting to .qlr (layer definition)

texts:
I added value maps (comboboxes) to some fields in QGIS 2.18.15 but when exporting to .qlr (layer definition) the comboboxes doesn't follow.
Any idea why?


texts: Generating Python Code from user clicks in QGIS

texts: Is there any way to auto-generate Python code in QGIS meanwhile I do clicks, so at the end I can copy/edit it the way I want?


texts: Query to Map GPS points with multi-linestring

texts:
I have a column of geometry of type multi-linestring in a database as in the following format

Now I want to add a new column of speed by mapping GPS points against geometry column (as described below that has been listed on separate CSV file).

Tell me how I query to merge the speed column (captured against GPS coordinate points) with existing geometries


texts:
I have some polygonal features that have minimum and maximum z values.

I'm trying to calculate the difference between a fixed extruded polygons and the sum of extruded related polygons.

Here's how I'm doing it:

Notes: the cursor is non recycling.

IExtrude extrude = new GeometryEnvironmentClass();

        while ((tempFeature = cursor.NextFeature()) != null)
        {
            var lower = Classifier.GetNamedLevelValue(tempFeature.get_Value(lowerIdx).ToString());
            var upper = Classifier.GetNamedLevelValue(tempFeature.get_Value(upperIdx).ToString());
            var geometry = extrude.ExtrudeFromTo(lower, upper, tempFeature.ShapeCopy);
            geometries.Add(geometry);
        }

        IGeometryFactory fact = new GeometryEnvironmentClass();
        IGeometry unioned = null;
        fact.CreateEmptyGeometryByType(esriGeometryType.esriGeometryMultiPatch, out unioned);

        for (int i = 0; i <= geometries.Count - 1; i++)
        {
            var topo = (ITopologicalOperator)unioned;
            var tempGeom = topo.Union(geometries[i]);
            unioned = (IGeometry)((IClone)tempGeom).Clone();
        }

        var extrudedParent = extrude.ExtrudeFromTo(0, 999, feature.ShapeCopy);
        var topoParent = (ITopologicalOperator)extrudedParent;
        var remaining = topoParent.Difference((IGeometry)unioned);

I get a COM Exception in topoParent.Difference. I'm not sure what is the problem. Can anyone help me? Any other ways to do it?

This is the COMException: - "Exceção de HRESULT: 0x80040229"}   System.Exception {System.Runtime.InteropServices.COMException}

I even whipped up a little unit test to check if it was my code or if the operation is illegal. Still crashes in the same spot.

    [TestMethod]
    public void VolumeTest()
    {
        object missing = Type.Missing;

        var p1 = PointHelper.ConstructPoint(0, 0);
        var p2 = PointHelper.ConstructPoint(0, 1);
        var p3 = PointHelper.ConstructPoint(1, 1);
        var p4 = PointHelper.ConstructPoint(1, 0);

        var p5 = PointHelper.ConstructPoint(0, 2);
        var p6 = PointHelper.ConstructPoint(2, 2);
        var p7 = PointHelper.ConstructPoint(2, 0);

        IPointCollection pointCol = new PolygonClass();
        pointCol.AddPoint(p1, ref missing, ref missing);
        pointCol.AddPoint(p2, ref missing, ref missing);
        pointCol.AddPoint(p3, ref missing, ref missing);
        pointCol.AddPoint(p4, ref missing, ref missing);
        pointCol.AddPoint(p1, ref missing, ref missing);

        IPointCollection pointCol2 = new PolygonClass();
        pointCol2.AddPoint(p1, ref missing, ref missing);
        pointCol2.AddPoint(p5, ref missing, ref missing);
        pointCol2.AddPoint(p6, ref missing, ref missing);
        pointCol2.AddPoint(p7, ref missing, ref missing);
        pointCol2.AddPoint(p1, ref missing, ref missing);

        IExtrude ext = new GeometryEnvironmentClass();
        IGeometry extruded = ext.ExtrudeFromTo(0, 1, (IPolygon)pointCol);
        IGeometry extruded2 = ext.ExtrudeFromTo(0, 2, (IPolygon)pointCol2);

        ITopologicalOperator topo = (ITopologicalOperator)extruded2;
        IGeometry result = topo.Difference(extruded);
    }

Thanks


texts:
OK so last time had a big issue with running WFS. Seems the issue was neither with my OpenLayers code nor GeoServer. The issue was with the servers not able to communicate with each other. The tomcat 6 server on which GeoServer was running and Apache server which is the default web server for my Ubuntu.

After some digging around it seems that I need to have tomcat 6 pre-installed on my machine before installing OpenGeo Suite. So having reinstalled the package, at end, the terminal gives me the following message

NOTICE: Please run /usr/share/opengeo-suite/tomcat-admin-setup.sh to complete this
installation.

So I executed the following command (from the required directory),

sh tomcat-admin-setup.sh

And I am getting this error,

tomcat-admin-setup.sh: 143: [: Ubuntu: unexpected operator
Unknown distro. The OpenGeo Suite is supported on Red Hat Enterprise Linux 5 and 6,
CentOS 5 and 6, and Ubuntu Linux 10.04 or later.

Does this mean the OpenGeo Suite is not compatible with Ubuntu 12.04??


texts:
First, I am using QGIS 2.8.2.

I have GPS waypoints and satellite imagery and am trying to draw polygons. This is easy enough, but I need to make sure each polygon is 1ha exactly. Is there a way to create polygons using the given area first?

I am creating polygons to delineate areas for community members' rights to do agro-forestry, so the final area of each polygon needs to be exact to be fair to everyone. The shape of the polygons will be rectangles, but not necessarily all right angles due to the terrain.


texts:
How to get KML attribute data on onload of html page using openlayers.

This is kml file: click here to get kml file


texts: How do I set up a one layer custom tile in Leaflet

texts: I am new to Leaflet and I have a large image which I have been attempting to break up into tiles and save them using the correct filenames so that they can be displayed using Leaflet's tileLayer.  I am trying to keep it simple first by breaking up the large image into just one layer of tiles.  But I have not been able to get Leaflet to load the tiles correctly.  Can anyone show me how a one layer custom tile should be set up in Leaflet?


texts: Creating a DEM from CAD contour polylines

texts:
Is it at all possible generate a digital elevation model raster from CAD contour polylines with elevations?

I am using straight AutoCAD 2014, not Map3D or Civil3D.

I was able to bring the dxf of the contours into QGIS, but was not able to get the elevations into the attribute table.

If I could get the elevations to transfer as an attribute, I think I could create the DEM from that, but as of now, I'm stuck at the polyline layer with no elevation attribute.

Any thoughts on how to get attributes to transfer from the CAD file, and convert the lines?


texts: Is it possible to create an average image of the analyzed period after the removal of the clouds to fill the flaws in the images?

texts:
// study area
var local = ee.Geometry.Rectangle([-42.2974, -20.5299,-40.9406, -19.1815]);
Map.addLayer(local);

var landsat8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')

// Fmask classification values var FMASK_CLEAR_GROUND = 0; var FMASK_WATER = 2; var FMASK_CLOUD_SHADOW = 3; var FMASK_SNOW = 4; var FMASK_CLOUD = 5;

var mosaic = landsat8 .filterBounds(local) .filterDate('2017-08-01', '2017-10-11') .mosaic();

var getQABits = function(image, start, end, newName) {
    // Compute the bits we need to extract.
    var pattern = 0;
    for (var i = start; i <= end; i++) {
       pattern += Math.pow(2, i);
    }
    // Return a single band image of the extracted QA bits, giving the band
    // a new name.
    return image.select([0], [newName])
                  .bitwiseAnd(pattern)
                  .rightShift(start);
};

// A function to mask out cloudy pixels.
var cloud_shadows = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 3,3, 'Cloud_shadows').eq(0);
  // Return an image masking out cloudy areas.
};

// A function to mask out cloudy pixels.
var clouds = function(image) {
  // Select the QA band.
  var QA = image.select(['pixel_qa']);
  // Get the internal_cloud_algorithm_flag bit.
  return getQABits(QA, 5,5, 'Cloud').eq(0);
  // Return an image masking out cloudy areas.
};

var maskClouds = function(image) {
  var cs = cloud_shadows(image);
  var c = clouds(image);
  image = image.updateMask(cs);
  return image.updateMask(c);
};

var mosaic_free = maskClouds(mosaic);

var visParams = {bands: ['B4', 'B3', 'B2'],min: [0,0,0],max: [2000, 2000, 2000]};

Map.addLayer(mosaic, visParams, 'With clouds');
Map.addLayer(mosaic_free, visParams, 'Cloud free');

// Average image

var median = ee.Image(mosaic_free.median())

Map.addLayer(median,
 {min:0, max: 3000, bands: 'B3, B2, B1'},
 'median composite of cloud free images');

// NDVI

var red = median.select('B4');
var nir = median.select('B5');

var calculo = nir.subtract(red).divide(nir.add(red));

//pallates NDVI

  var ndvi_palette = 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,'+
  '3E8601, 207401, 056201, 004C00, 023B01,012E01, 011D01, 011301';

var ndvi = calculo.clip(area);

Map.addLayer(ndvi, {min:-0.1, max: 1.0, palette: ndvi_palette}, 'NDVI');

// Visualization ------------------------------------------------------------------
Map.centerObject(area, 8);

After building a cloud mask, is it possible to create an average image of the period without the presence of faults to calculate the NDVI? I have adapted the link code for my study area, but I am not able to create a averege image without the presence of clouds to calculate the NDVI.


texts: Maintain Feature Data without Creating It

texts:
I work in City Gov. We utilize County Assessor data (Parcels) as they are updated from the Assessors office. Currently we store then in our Enterprise SDE system. As we get new updates, we have to shut down all users of Parcel data in order to delete the old set and overwrite with the new update.

Is there a way to update this data without having to shut down users? Some kind of versioning concept for maintainers vs. creators of gis data?

We are currently looking into TRUNCATE AND APPEND as an option.


texts: EGM2008 Global 1-Minute Geoid Undulations to raster

texts:
The EGM2008 Global 2.5 Minute Geoid Undulations is available in GIS raster format for download from https://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_gis.html
A 1-minute grid exists as well. The software and the grid can be found at https://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html./ It is the file named as
1 x 1-Minute Geoid Undulation Grid in WGS 84 - SMALL ENDIAN (Und_min1x1_egm2008_isw=82_WGS84_TideFree_SE.gz - 825 MB) - 1 x 1-minute global grid of pre-computed geoid undulations. This file has a SMALL ENDIAN internal binary representation.
or
1 x 1-Minute Geoid Undulation Grid in WGS 84 - BIG ENDIAN (Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz - 828 MB) - 1 x 1-minute global grid of pre-computed geoid undulations. This file has a BIG ENDIAN internal binary representation.
I am looking to get this 1-minute EGM2008 grid into raster format that can be used in GIS or for a download.


texts: I'm new to Collector. I have a crew of people go out into the field to do offline collection at several sites. They downloaded a base map (activating sync) then noticed after the first day of data collection that the aerial needs to be updated to cover a broader area. I can't figure out how to download a new base map for offline use. Obviously we need internet connectivity, but even with connectivity there doesn't appear to be an easy way to get a new base. Anybody had this issue and found a work around or solution? I figure it should be something very simple. Other then making a totally new map and sharing it with the field staff, I can't figure out a simple work around.


texts: Scalebar doesn't match up with actual dimensions

texts: Every time I make a map, the scalebar doesn't match with what the real dimensions of my map should be. I've been doing maps of Chicago in CRS: NAD83 and every time I put the map in the print composer and add the map and scalebar, it says that Chicago is .04m... Which is obviously not correct in the slightest. Also, I am unable to change the units in my map.


texts:
I just installed GeoNode and really love the manner in which it helps me arrange the contents related to the map but I am not sure about adding OpenLayers to it.

Well here is where I am faltering,

I am not sure where the OpenLayers code goes? I felt that I could use the code to customize GeoNode's look and feel and make use of the customization features on the Web Page, but to publish the Map you receive a embed link which only gives you the Map and not the other features.
If we have to use OpenLayers separately, like in the HTML, I am not sure, How to run that static html file, since the localhost is being used by the GeoNode.


texts: Obtaining bearings to polygon extents in ArcGIS?

texts:
In the context of a wider question, (Generating multiple visibility bearings in ArcGIS as automated procedure?) is it possible to extract bearing information from a particular point or set of points, to extreme tangents of a polygon?

For example, I have a point shapefile with a set of viewpoints (small red dots in image below). I also have a polygon shapefile containing the buffer of another set of points, which outline the exact extents, in m, of my development points (shown in blue below.)

I want to know the extents, in bearing format, of the furthest left tangent point to the furthest right tangent point of the buffer shapefile of my development points. (the red arrows below, i.e. 275 - 282 degrees.)

.

Can I do this in GIS, and if so how / what tools should I be looking at?


texts: Arcpy Error 999999: Error executing points to line. Coordinates or Measures are out of bounds

texts:
I am selecting two points from a table and then doing a points to line function. I was originally getting this error because they were multipoints with multiple locations. I ran Multipart to Single Part on that layer and now am getting the same error on only a select few points but there are no multiples.

What else can cause this error? Thanks!

temporarytodraw holds 2 points and that's the result I get. I'm sorry I don't have much more to post than this. It selects points that match certain criteria and places them in the temporarytodraw but there is an if statement check to kick out any records that match less than 2 or more than 2.

arcpy.PointsToLine_management("temporarytodraw", "offset_line_thing", "", "", "NO_CLOSE")


texts: What type of OS is best for ESRI mobile Development?

texts: I know there is going to be varied opinions on this question. We're being asked to create some mobile maps for various departments in our county. This happens to come at a time when AT&T has bought out our wireless provider and we're all due for new phones. I'm trying to figure out which phone operating system is best for ESRI mobile development. ESRI has the API out for the iOS, and Windows Mobile. And it looks as if Android is on its way. Any opinions on which is better and why would be great.


texts: Making QGIS select correct feature id in PostGIS materialized view?

texts:
Can you tell me if there is a way to make QGIS automatically select the correct feature ID when I choose a materialized view in PostGIS?

It's not a big deal, but it is nonetheless something I would like to avoid. I tried to create an unique index on the fid column, but nothing changed.


texts: Keeping geometry column from both dataframes when applying sjoin() using GeoPandas

texts:
I am trying to use sjoin() to spatially join two geo dataframs and then use both geometries for later spatial calculation, however, I noticed that from the official docs, only one geometry column from either side may be kept.
df= gpd.sjoin(gdf1, gdf2, how='left', op='contains')

Is there a way to manage that?


texts: Creating user defined input to simplify Reclassify tool and negate need for remap

texts:
I would like to build a script tool using Arcpy that reclassifies a raster layer based on a single user input parameter. For example, if the user inputs the number 7, then 0 - 7 would reclassify as 1, anything > 7 would become NoData.

The reason I need this is; I am building a site selection tool which contains many processes. In the middle of the tool there is a Reclassify process. The workflow goes something like this:

Elevation > Slope > Slope + 0.5 > Int > Reclassify > Raster to Polygon - This is the final vector that becomes part of other vector analysis to perform a site selection. It is a polygon of all areas of suitable slope.

I have been playing around with various attempts but to no avail, mostly because the Reclassify tool will not take a variable as part of the input. Everything else in the tool works fine but I don't know to have a user-defined input determine the Reclassify process.

# Import arcpy module
import arcpy
from arcpy.sa import *

# Script arguments
Elevation_Data = arcpy.GetParameterAsText(0)
Slope_Value = arcpy.GetParameterAsText(1)

# Local variables:
Slope_1 = "D:\egm722\Data\Scratch.gdb\\Slope_1"
Int_Slope = "D:\egm722\Data\Scratch.gdb\\Int_Slope"
Reclass = "D:\egm722\Data\Scratch.gdb\\Reclass"
RasPoly = "D:\egm722\Data\Scratch.gdb\\RasPoly"

# Process: Slope
arcpy.gp.Slope_sa(Elevation_Data, Slope_1, "DEGREE", "1", "PLANAR", "METER")

# Process: Raster Calculator
Slope_2 = Raster(Slope_1) + 0.5
Slope_2.save("D:\egm722\Data\Scratch.gdb\\Slope_2")

# Process: Int
arcpy.gp.Int_sa(Slope_2, Int_Slope)

# Process: Reclassify
arcpy.gp.Reclassify_sa(Int_Slope, "VALUE", "0 Slope_Value 1", Reclass, "NODATA")

# Process: Raster to Polygon
arcpy.RasterToPolygon_conversion(Reclass, RasPoly, "SIMPLIFY", "VALUE")

Here, I have attempted the variable Slope_Value to be the user-defined input as a script tool parameter type 'Long', but of course this does not work.

Traceback (most recent call last):
File "D:\egm722\elev.py", line 26, in <module>
arcpy.gp.Reclassify_sa(Int_Slope, "VALUE", "0 Slope_Value 1", Reclass, "NODATA")
File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\geoprocessing\_base.py", line 510, in <lambda>
return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 000622: Failed to execute (Reclassify). Parameters are not valid.
ERROR 000628: Cannot set input into parameter remap.

I cannot skip the reclassify process and delete values from the 'Raster to Polygon' output (using Select tool for example) because converting non-reclassified Raster layers using 'Raster to Polygon' takes far too long for it to be considered a viable tool.


texts: How to run a 2.0 ESRI DLL in a 4.0 Addin application

texts:
We have a complex set of operations running as Add-in's and services that span both ArcMap, ArcServer, and ArcGIS Explorer for Desktop.  Currently we are compiling everything under .NET 3.5 but we may need to move to .NET 4.0 to resolve an issue.

All of our code uses ProgressHelper from Explorer for Desktop 10.1 which is strictly a 2.0 framework applicaiton.  Normally I'd add the  command to the config file but the Add-In's are under ESRI control, not the .NET framework.

Is there a way to modify the ESRI configuration files so that the .NET framework knows that 4.0 is the supportedRuntime environment and allows me to load the 2.0 framework DLL?

As a testing note: ArcMap does have it's own configuration file and in that I do have the Startup section defined for both the 2.0 framework and the 4.0 framework.


texts: Obtaining Free High Resolution Satellite Imagery to use in QGIS

texts: I just downloaded QGIS 3.10 Mac version the other day and I am still very green with the program.  I was wondering what is the best way to obtain high resolution satellite imagery?  I know there are some websites you can pay for it, but I was wondering if there are any ways you can get it for free?


texts: Clickable HTML link labels for features in exported PDF

texts: Is it possible to create a map in QGIS (version 2.18.15) exported to PDF that contains features with labels as an HTML link? By clicking on feature symbol or label this link should open in browser.


texts:
I am trying to extract the 24 hour precipitation data of a given polygon by creating a dictionary . I have gotten results but i'm not sure if the parameters I have inserted in the reducer are correct.

The resolution that is given on GEE dataset for TRMM is 0.25 arc, so I believe it is not in WGS84, so i'm not sure if I should put 0.25 or to calculte it in km.
As both give me numbers, I don't know which one is correct

var meanDictionary =totalRain.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: 0.25,
  maxPixels: 1e9
});

Then, I have used the reducer sum() but i'm not sure if it will give me the sum of all the pixels I have on th eone image or for each pixel it will give me the sum of all the 8 umages I have used.

To better understand my cocnerns, I recommend to see my code:
https://code.earthengine.google.com/39edb9298606d1912331e5d7f7d7723a

the polygon I have used is located in north of Thailand, and is random, so any polygon there would work


texts: How can I add a combination of digits and alphabets using the edit symbol tool in ArcGIS. I can add a single number/alphabet for example 2 to it as shown in the screenshot, but I don't know how to add the whole combination to the symbol of that particular point which is for example 222 or 222A.


texts: Differences between SQL dialects in QGIS vs ArcView vs Oracle

texts:
I'm kinda new with SQL queries and there is something I don't get: A SQL query in QGIS is written differently than the same query in ArcView and will probably also be written differently if I want to query an Oracle db using an Oracle tool. So is it like SQL is a language and each software uses its own slang of the same language ?

EDIT:
Thanks for those answers!
Now do you guys know of a documentation somewhere that would teach me how to make a query in QGIS field calculator. I can make a simple one using the implemented functions, but what about something more complex (using case when, where, and so on) ?


texts:
I'm using QGIS 3.4.5 and created a line shape layer.
I have another line layer in my project that I want to create and trace the geometry of the other layer.
I activated the "add line object" button and also checked out the "node tool".
Snapping is enabled in settings --> options.
Using another CRS for the layers doesn't help either.
Can somebody tell me how I can activate the "trace tool" in "snapping toolbar"?
Proceeding like this, the "trace tool" isn't active.
Kind regards, Larissa


texts: In ArcGIS, I discovered the sw.lock files? What are they?

texts:
In ArcGIS, sometimes I get a sw.lock lock file with the name of my shapefile. Example : Polygon_shape.shp.sw.lock

I know about the  .lock files (sr, ed, rd, wr) but I have never seen "sw" in the past. I'd like to know where it came from since I get it using my ArcObjects application and it is intermittent.


texts: Symbol Rotation with Lines Bulk using ArcGIS Desktop?

texts:
I have feature points with specific symbol like square.

Two lines join with 1 symbol square.

I want rotate all symbols along to lines in ArcGIS Desktop 10.6.

Is there any tool or any option to fix it?


texts: Creating gpkg with many layers using ogr2ogr

texts:
I am trying to merge many shapefiles into a single gpkg using ogr2ogr (version 3.0.2).

The command I am using is:

ogr2ogr dst.gpkg src.shp -nln layerOne
ogr2ogr dst.gpkg src.shp -nln layerTwo

BUT when I run ogrinfo it says my dst.gpkg has only layerTwo.


texts: GRASS d.mon failing to connect to a monitor

texts:
I'm running GRASS 6.4.6 on Ubuntu 16.04. I can't seem to get the monitor to work.

d.mon -L

name            description                    status
----            -----------                    ------
PNG             Create PNG file                not running
png1            Create PNG file                not running
png2            Create PNG file                not running
png3            Create PNG file                not running
png4            Create PNG file                not running
png5            Create PNG file                not running
png6            Create PNG file                not running
png7            Create PNG file                not running
gism            Create PNG file for gis.m      not running
HTMLMAP         Create HTML Image Map          not running
PS              Create PS file                 not running
ps1             Create PS file                 not running
ps2             Create PS file                 not running
ps3             Create PS file                 not running
ps4             Create PS file                 not running
ps5             Create PS file                 not running
ps6             Create PS file                 not running
ps7             Create PS file                 not running
(Thu Sep 28 21:12:20 2017) Command finished (0 sec)

I thought I could go 'd.mon start=X0', but it returns:

ERROR: no such monitor 'X0'
ERROR: No such monitor as <X0>
Problem selecting X0. Will try once more
ERROR: No such monitor as <X0>

X11 is installed and working, but how do I get monitors to appear in the list above?


texts: How to create vector-based maps showing historic borders

texts:
I'm an absolute beginner in map making.

I want to create digital maps with European historic borders that all have the same projection and style. All the data should be vectorized.

My idea is to have a blank digital map on which I draw the borders. To create the borders I take an existing map that shows the historic borders.

My questions:

Is there a collection of hgis data (freely) available so that I don't have to do it myself?
Is there a good tutorial for my task?
How do I map the map showing historic borders to the blank map?
How do I copy the border map? Should I use a tablet an pin a print of the map on it?


texts: Is additional SQL syntax necessary to account for a field with a domain/coded values?

texts:
I'm attempting to create a spatial view.  The following statement returns an empty view, even though there are "Flushed" records in the MAINT_ACTIVITY field.

    SELECT h.OBJECTID, h.FACILITYID, h.SHAPE, r.Hydrant_ID, r.MAINT_ACTIVITY, r.MaintDate, r.GlobalID
    FROM wHydrant h
    JOIN wHydrant_Maint r
    ON h.FACILITYID = r.Hydrant_ID
    WHERE r.MAINT_ACTIVITY = 'Flushed'

If I remove the last WHERE clause, the query will return a view of all my hydrants.  If I change the WHERE clause to query another field, it will return correct results.  I'm only having issues with this specific MAINT_ACTIVITY field and the only major difference is that this field has a domain with coded values.  However, the coded values match the descriptions.  So, Flushed is the code and Flushed is the description.  Is additional syntax still necessary?  Why doesn't the above statement work?

*Update--When I query other fields in the Hydrant_Maint table, all the MAINT_ACTIVITY field results return as NULL, even though there are entries in the original table.


texts: Publishing in QGIS Server (WMS, WFS-T) only certain map extent

texts:
I am trying to publish a layer in WFS-T but I only want to publish a certain zone ( 10% of the whole layer ).

So I choose my extent in QGIS canvas and I go to Poryect Properties / OWS Server and I click in "Use the current view extension" and save project.

If I load the layers using WFS-T I can see the whole map and not the chosen extensión.

What am I missing?


texts: Publishing Map on ArcGIS Server Manager without using ArcGIS for Desktop?

texts:
I want to publish my maps on sever manager.

My map data is in PostgreSQL and I want to connect it without using ArcGIS Desktop.

Is there any possibility to publish service without involving ArcGIS Desktop?


texts: Why can't I load more than 100,000 Features from a WFS-Layer in ArcMap?

texts:
I want to have all the trees in Berlin as Shapefile. But I only get the first 100,000 trees.

This is the information about the WFS-Layer: http://daten.berlin.de/datensaetze/baumbestand-berlin-stra%C3%9Fenb%C3%A4ume-wfs

This is the WFS-Adress: http://fbinter.stadt-berlin.de/fb/wfs/geometry/senstadt/re_strassenbaeume/?service=wfs&request=GetCapabilities

In the Parameters I changed f.e. the constraint for max features from 100,000 to 900,0000 or 1,000,000 but the result is the same.


texts: Retrieve image coordinates (X, Y) from geographic coordinates in NetCDF file using Numpy

texts:
I'm trying to retrieve image coordinates (X, Y) from geographic coordinates in NetCDF file (Sentinel 3 data) using Numpy. The problem is that the final pixel is too close to the desired pixel but not the closest to it (about one or two pixel far away from the closest in x, y or both x,y). I'm using Numpy's norm(),  min() and where() function. The code is presented below:
from netCDF4 import Dataset
import numpy as np

# Read the NetCDF file with the geographic coordinates information
nc_coord = Dataset('coords.nc', 'r')

lat = nc_coord.variables['latitude'][:]
# Python buffer object pointing to the start of the array’s data
lat = lat.data
lon = nc_coord.variables['longitude'][:]
lon = lon.data

# define your point:
target_lon, target_lat = 23.779675, 35.33786
# lat and lon are x,y matrices, we will add another dimension z in order to stack them
lat = lat[:, :, np.newaxis]
lon = lon[:, :, np.newaxis]
grid = np.concatenate([lat, lon], axis=2)
vector = np.array([target_lat, target_lon]).reshape(1, 1, -1)
subtraction = vector - grid

# vector subtraction
dist = np.linalg.norm(subtraction, axis=2)

# query where in the 2D matrix is the closest point to your points of interest (?)
result =  np.where(dist == dist.min())

# this is your lat/lon geo point translated to a x, y position in the image matrix
target_x_y = result[0][0], result[1][0]

Can you figure out what is happening?


texts: writeRaster() not preserving EPSG information

texts:
I need to upload some rasters in .tif format to a geoserver that requires EPSG information to correctly display the data.  When I assign an EPSG using crs() to a raster that is in memory, the WKT information reflects that EPSG code correctly.  However, when I write that raster to disk and read it back into R to double check, the EPSG code is not preserved and the WKT information is different (though the crs is the same).
memoryCRS <- raster("wiltpoint.tif")
crs(memoryCRS) <- "EPSG:5070"
writeRaster(memoryCRS, "wiltpoint2.tif")
diskCRS <- raster("wiltpoint2.tif")

cat(wkt(memoryCRS))
#> PROJCRS["NAD83 / Conus Albers",
#>    BASEGEOGCRS["NAD83",
#>        DATUM["North American Datum 1983",

cat(wkt(diskCRS))
#> PROJCRS["unknown",
#>     BASEGEOGCRS["unknown",
#>         DATUM["North American Datum 1983",

all.equal(wkt(memoryCRS), wkt(diskCRS))
#> [1] "1 string mismatch"

identicalCRS(memoryCRS, diskCRS)
#> [1] TRUE

compareCRS(memoryCRS, diskCRS)
#> [1] TRUE


texts: How to modify styles extracted from KML in OpenLayers4

texts: I'm using OpenLayers4 to load and display a 3rd party's KML.  I need to retain the original styles (colors, opacity, icons, etc), but I need to change the icon sizes.  Is there any way to modify the styles which get extracted or should I focus on the sharedStyles or is there a more appropriate way to do this?


texts: How to set coordinate system on raster catalog with arcpy?

texts: When creating a raster catalog, if the optional coordinate system fields for the raster and geometry columns are left empty, they will be stored as unknown. What geoprocessing functions in arcpy should I use to set both of them?


texts:
How do I make this map remove those lines and just view the map as one?

What I mean is I wanted to remove those lines so that the outcome will just be a plain map without those lines.


texts: Create average boundary between raster classes

texts:
I have a raster with multiple fields in it representing glaciers. As the boundary between land and ice is not often strict, there are little ice fields outside the end of the glacier. I do a raster classification and get it delimited.

Now I would like to estimate mean height of this boundary, and for it I need a line feature. How can I extract an average position from it? I may be confusing, but I think a picture will explain it well:

The boundary doesn't need to be a straight line, but it needs to be a single line. So I can't just calculate circumference of classes.

I'm using ArcGIS/arcpy and open source (QGIS, Gdal).


texts: Google Earth Engine FeatureCollection getDownloadURL returns empty CSV

texts:
I am trying to create a link that lets me or users of my app download a feature collection as a CSV directly from the map. However, when using ee.FeatureCollection.getDownloadURL, the csv that it downloads is always empty. When I print the feature collection to the screen there are features and attributes there. Any ideas as to why it's outputting an empty CSV?

Here is the part of my script that is outputting the download url:

  var output_fc = ee.FeatureCollection(leng.map(function(x){
    var dict = {class_name: to.get(x), class_area: class_areas_results.get(x), tree_cover_loss_area: tcl_results.get(x), tree_cover_gain_area: tc_gain_results.get(x),
    avg_tree_area: avg_tree_area_results.get(x), tcd_00: tcd_00_results.get(x), tcd_10: tcd_10_results.get(x)};
    var feature = ee.Feature(null, dict);
    return feature;
  }));

  var download = (ui.Label('Download').setUrl((ee.FeatureCollection(output_fc).getDownloadURL({
    format: 'CSV',
    filename: 'changeAreas'
  }))))

  panel.add(download)

  print(output_fc)


texts: Meaning of different elements of delivered json-Data-Structure when geocoding with Python and OpenStreetMap

texts:
I'm trying to do some geocoding in Python with the library geocoder and OpenStreetMap:
import geocoder
g = geocoder.osm('11 Wall Street, New York')
print(g.json)

Where can I get information about the different elements of the json-Dictionary?
The json-Dictionary:
{
   'accuracy': 0.511,
   'address': 'New York Stock Exchange, 11, Wall Street, Financial District, Manhattan, New York County, New York, 10005, United States',
   'bbox': {
       'northeast': [40.7073927, -74.0108331],
       'southwest': [40.7066568, -74.0115875]
   },
   'city': 'New York',
   'confidence': 10,
   'country': 'United States',
   'country_code': 'us',
   'county': 'New York County',
   'housenumber': '11',
   'icon': 'https://nominatim.openstreetmap.org/ui/mapicons//poi_point_of_interest.p.20.png',
   'importance': 0.511,
   'lat': 40.7070653,
   'lng': -74.01117438284459,
   'neighborhood': 'Financial District',
   'ok': True,
   'osm_id': 278053530,
   'osm_type': 'way',
   'place_id': 158065100,
   'place_rank': 30,
   'postal': '10005',
   'quality': 'attraction',
   'raw': {
       'place_id': 158065100,
       'licence': 'Data © OpenStreetMap contributors, ODbL 1.0. https://osm.org/copyright',
       'osm_type': 'way',
       'osm_id': 278053530,
       'boundingbox': ['40.7066568', '40.7073927', '-74.0115875', '-74.0108331'],
       'lat': '40.7070653',
       'lon': '-74.01117438284459',
       'display_name': 'New York Stock Exchange, 11, Wall Street, Financial District, Manhattan, New York County, New York, 10005, United States',
       'place_rank': 30,
       'category': 'tourism',
       'type': 'attraction',
       'importance': 0.511,
       'icon': 'https://nominatim.openstreetmap.org/ui/mapicons//poi_point_of_interest.p.20.png',
       'address': {
           'tourism': 'New York Stock Exchange',
           'house_number': '11',
           'road': 'Wall Street',
           'neighbourhood': 'Financial District',
           'suburb': 'Manhattan',
           'county': 'New York County',
           'city': 'New York',
           'state': 'New York',
           'postcode': '10005',
           'country': 'United States',
           'country_code': 'us'
       }
   },
   'region': 'New York',
   'state': 'New York',
   'status': 'OK',
   'street': 'Wall Street',
   'suburb': 'Manhattan',
   'type': 'attraction'

}
Some elements are of course obviously clear like lat and lng.
What does a value of 0.511  for accuracy  mean?
What's the meaning of ok?
What does a value of 10 for confidence mean?
So for short: where can I get more info on this?


texts: Implementing Custom Selection Tool (Adding to current selection) with ArcObjects?

texts:
I am trying to implement a custom selection tool in ArcMap 9.3.1 that is functionally similar to the default ArcMap selection tool with "Adding to current selection" enabled. ie. the user can click to incrementally select individual features one by one, while adding on to previously selected features. When the user clicks on an already selected feature, that feature should be unselected.

How should I go about replicating and implementing this behavior for my custom tool? There are some selection-related code snippets lying around, but nothing really comes close. I was able to perform single selections by calling SelectByShape(geom, null, true), where geom = a small area around the clicked point. To support multiple selections at various locations, how then do I keep the previously selected features still selected as I select a new one, since SelectByShape(geom, null, true) only selects a single feature each time it's called.

Clearly, I will need another approach, which probably involves me keeping a growing list of currently selected features, but I am not sure if there is a more straightforward solution. Has anyone attempted this before?


texts: Help me with this code from GeoExt and OpenLayers

texts:
I haven't gone entirely through GeoExt. But I am trying to understand the below code, but I am struck. The below code is about creating popups of feature information. But the layers added are WMS Layers, How can a WMS layer contain feature information. It is not a vector layer.
     <script>
        OpenLayers.ImgPath = "../../resources/openlayers/img/";

        // the `panel` and `popup` variables are declared here for easy debugging
        var panel, popup;

        Ext.onReady(function() {

            // queryable WMS layer
            var zoning = new OpenLayers.Layer.WMS(
                "Medford Zoning",
                "/geoserver/wms",
                {layers: "medford:zoning", format: "image/gif", transparent: "true"}
            );

            panel = new GeoExt.MapPanel({
                title: "MapPanel",
                renderTo: "map-id",
                height: 300,
                width: "100%",
                layers: [
                    new OpenLayers.Layer.WMS(
                        "Global Imagery",
                        "http://maps.opengeo.org/geowebcache/service/wms",
                        {layers: "openstreetmap", format: "image/png"}
                    ),
                    zoning
                ],
                center: [-122.87, 42.34],
                zoom: 13
            });

            // create a control to get feature info from queryable layers
            var control = new OpenLayers.Control.WMSGetFeatureInfo({
                url: "/geoserver/wms"
            });
            panel.map.addControl(control);
            control.activate();

            // register a listener for the getfeatureinfo event on the control
            control.events.on({
                getfeatureinfo: function(event) {
                    // close existing popup
                    if (popup) {
                        popup.destroy();
                    }
                    popup = new GeoExt.Popup({
                        title: "Popup",
                        map: panel.map,
                        lonlat: panel.map.getLonLatFromPixel(event.xy),
                        width: 250,
                        autoScroll: true,
                        collapsible: true,
                        bodyStyle: {padding: 5},
                        html: event.text
                    });
                    popup.show();
                }
            });

        });

    </script>

The map is like this :

UPDATE
I changed the zoning layer to some other WMS layer in GeoServer. Now I can load the layer and everything else except that when I click, popup button opens but feature information is not coming.
I found in GeoServer documentation that there is GetFeature request availability for WMS layers. But there are some settings to be changed and I couldn't understand those. Please help changing those settings so that popup can show the feature information of the layer i.e., GetFeature request is accessible.
UPDATE
This is that GeoServer manual link where I read about WMS settings:
http://docs.geoserver.org/stable/en/user/webadmin/data/layers.html
UPDATE
Found it! We should use proxy settings for features unlike normal layers. So I should replace the URL with protocol and formats.


texts: PostGIS topology hierarchal topogeometries

texts:
I'm seeking guidance on how to correctly organise hierarchal topogeometries.

I'm working through PostGIS in Action's chapter on topology and, in summary, it sets up two tables with topogeometries, one called neighbourhoods (which becomes layer 1) which is based on geometries loaded with toTopoGeom, and one called cities (layer 2), which is derived from the neighbourhoods (layer 1).

The book says this about layer 2.

  "...You then define a column called topo in the cities table to define the neighbourhoods that each city is composed of.  Each city is defined by a topogeometry. In the example, you only want one city, so there's only one topogeometry"

I find this confusing. I thought the column definition was the topogeometry in which case how can each Row entry (a city) have a different one?

What say I want more than one city?  Where do the neighbourhoods go.  How do I say that 'this city has these neighbourhoods from layer 1, and this second city, has these different neighbourhoods from layer 1'?

I'm clearly missing something important, because at the moment it seems to me that every 'parent' in the hierarchy basically ends up as a single solitary row in a table.


texts: Incomplete export of data from a Postgres table

texts:
I'm having a serious data export problem with the PGAdmin widget and pgsql2shp (QGIS too). The problem lies within converting a PostGIS table to shapefile. I put in 3234 feature counts...and I get back 3221.

I get limited information from the error reports (see below...QGIS gives me the most information...but it is minimal). The errors (when available) reference ogr....but not much else. Ironically, the only way I've found to actually export the full record set is by using ogr2ogr.

I'm not sure why round-tripping with the PGAdmin tool is failing so bad. It seems to be able to import the initial 3234 - but cannot export the full set.

Here is the data (2 MB unzipped on a cloud drive).

The process I've been following..

1) Load the shapefile into Postgres using the PGAdmin tool (WORKS!)

works fine...notably it is loaded as a MULTIPOLYGON[2] EVEN THOUGH the Shapefile type is Polygon

3234 record counts

Importing with configuration: as_clean_test_ogr_export_utf8, public, geom, H:\as_clean_test_ogr_export_utf8.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=4326
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Shapefile import completed.

2) EXPORT... (FAILS!)

UPON Export (using >pgsql2shp)

H:\>pgsql2shp -f E:/as_clean_test_ogr_export_utf8.shp -h <hostname> -u <username> -p <port> -P <password> <postgrestable> "as_clean_test_ogr_export_utf8"

Initializing...
Done (postgis major version: 2).
Output shape: Polygon
Dumping: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX [3221 rows].

QGIS (Save As in QGIS 2.6.0)

Export to vector file failed.
Error: Feature write errors:
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Feature geometry not imported (OGR error: )
Only 3221 of 3234 features written.

Export Using PostGIS 2.0 Shapefile and DBF Loader Exporter

3221 record counts

Exporting with configuration: as_clean_test_ogr_export_utf8, public, as_clean_test_ogr_export_utf8
Done (postgis major version: 2)
Output shape: Polygon
Table export completed.

3) EXPORT (WORKS! WHY?!)

3234 record counts

    ogr2ogr -f "ESRI Shapefile" test_ogr_export.shp PG:"host=<host> user
    =<usr> dbname=<thatdb> password=<password>" "as_clean_test_ogr_export_utf8"


texts: Load shapefiles stored in web page in QGIS?

texts:
I prepared a map with QGIS and published it via QGIS cloud free. Some layers are quite big, and the cloud size 50MB limit has been reached. If the layer sources could be stored in my web page --- reading shapefiles directly into QGIS, instead of uploading the data into the cloud server --- that would spare storage space in my qgiscloud.com account.

I noticed that, after uploading all data to the qgiscloud server, one layer's source reads thus:

dbname='wmsdkw_hjbbxg' host=db.qgiscloud.com port=5432 user='wmsdkw_hjbbxg' key='id' estimatedmetadata=true table="mexSECCION" (wkb_geometry) sql=

How can I tweek the layer source so that it is read directly from http:/ericmagar.com/data/maps/shp/mex/SECCION.shp?

Is there an alternative way to share a large-ish QGIS map with my students?


texts: Unknown geometry type in postgis on ST_contains

texts:
I am trying to match a geometry of type Polygon to Point using ST_Contains:

UPDATE events
SET country = countries.geonameid
FROM countries
WHERE ST_Contains(countries.geometry, events.location)

The country table was created as follows:

CREATE TABLE countries (
    geonameid INTEGER PRIMARY KEY,
    geometry GEOMETRY(Geometry, 4326)
)
INSERT INTO countries (geonameid, geometry VALUES (%s, ST_GeometryFromText(%s, 4326)))

and the event table like this:

CREATE TABLE events
(
    id BIGINT PRIMARY KEY,
    location GEOMETRY(Point, 4326),
    country INT
)
INSERT INTO events (id, location) VALUES (%s, ST_GeomFromText('POINT(%s %s)', 4326)) ON CONFLICT DO NOTHING"

Unfortunately on the update query I obtain the following error:

  psycopg2.InternalError: Unknown geometry type: -671088639 - Point


texts: My legend is currently displaying both the layer name and indented data name underneath the layer name. How to display just its name shown in the table of contents?


texts: Get vector features inside a specific extent

texts:
The problem:

I have a very large vector layer with many features and a much smaller raster layer in a defined region. I want to get only those vector features which are within the raster layers extent (extent = (xmin,xmax,ymin,ymax) ).

Is there anything like a standard SQL-query (something like SELECT * FROM layer WHERE EXTENT < extent) or another command (some test?) i could use to get only those features that are within a given extent?

EDIT:
Added python code to do a bounding box intersection test for vector features using just gdal and ogr python binding


texts: I have some feature classes (7 feature classes) and I want to batch buffer in ArcGIS 10. I don't know how can I buffer at once.  Please help me.


texts: Go over array values in Google Earth Engine and perform actions based on values in another array

texts:
I am trying to perform actions in an array that are conditional to values in another array with equal length.
Example:
var values = ee.Array([827, 825, 823, 821, 819, 817, 814, 812, 810, 808, 806, 804, 802, 799, 800, 500])
var vertices = ee.Array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1])

// desired outcome: [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
-25, 0, 0, -302]

I want to subtract each value in values, for which vertices has a 1 in the same position, from the previous value in values where that is the case. An exception is the first element of values, which should be set to zero and is always 1 in vertices.
Furthermore, all values that have a 0 in the same position in vertices should be set to 0.
I have accomplished parts of what I want to do by using vertices as a mask for values and then subtracting a left-shifted version of values from a right-shifted one but cannot find a way that includes the zeros in the outcome. I looked through the documentation for ee.Array and it's functions as well but just can't figure it out.


texts:
I have an ArcPad 10.0 application for which I export tiff images from ArcMap (previously 10.1, now 10.4). The tiff images are created by rasterizing polygon layers. I then set the symbology and export to a new tiff image using the colormap of the tiff.

I generally don’t have a problem with this process, but sometimes when I add a tiff in ArcPad, I get the following message “Error opening …..tif: Error 0”. This happens only with my larger images and when I split the image in two or make the extent smaller, I don’t get the error.

I’ve run a series of processes to see if I can determine exactly what is happening: I’ve made the extent of the image smaller, resampled the image to a bigger cell size, changed the compression  and changed the pixel type from 8 bit to 1 bit (the cell values of the image is only 0 or 1).

What I’ve found is that it seems like the uncompressed size of the tiff has an influence on whether I get the error message or not. Whenever I changed the extent or resampled the image so that the uncompressed size is less than 2GB, I didn’t get the error message. Whether you put compression on the image or not does not seem to make a difference. The weird thing is however that when I change the pixel type to 1 bit I still get the same error message that I get when the uncompressed size of the 8 bit image is 2GB or more, even though the uncompressed size of the 1 bit image is much smaller.

I haven’t found any reference to a size limit for tiff’s in ArcPad. Does anybody know whether a size limit does exist and also why changing to 1 bit image didn’t make a difference? Or is there something else that is causing this problem?


texts: OpenLayers with ArcGIS REST API identify layer feature

texts:
We are using ArcGIS REST services on OpenLayers. The layers are placed on the map correctly but how to identify a layer feature on mouseover or mouse click and get the details of the particular layer feature?

Is there any function in OpenLayers library to do it..?


texts: Why geometryInstances not showing on Cesium 3D model?

texts:
I have the following code, but geometryInstances is not showing up on my 3D model.

         propertyPlacemark.push(viewer.entities.add({
      geometryInstances : new Cesium.GeometryInstance({
              geometry : Cesium.BoxGeometry.fromDimensions({

                dimensions : new Cesium.Cartesian3(8.0, 5.0, 8.0)
              }),
              modelMatrix : Cesium.Matrix4.multiplyByTranslation(Cesium.Transforms.eastNorthUpToFixedFrame(
    Cesium.Cartesian3.fromDegrees(parseFloat(block[2]),parseFloat(block[1]))), new Cesium.Cartesian3(0.0, 0.0, 0.0), new Cesium.Matrix4()),
              attributes : {
                color : Cesium.ColorGeometryInstanceAttribute.fromColor(new Cesium.Color(1.0, 0.0, 0.0, 0.5)),

              },

            })
    }));

*parseFloat(block[2])= longitude
parseFloat(block[1])= latitude

Expected outcome：


texts: display historical map without georeferencing?

texts: I have a historical map of Rome that was published as 46 separate sheets. I'd like to scan these sheets, mosaic them together, and display them as a single map. I've tried to georeference a few of the scanned sheets and it turns out that the original map is distorted with respect to actual geography, so georeferencing the sheets results in the sheets overlapping or being skewed. I am not interested in projecting the historical map onto a basemap, just stitching them together. Is there a way to mosaic these images without georeferencing, or a way of georeferencing them without causing distortion?


texts:
I have installed the latest QGIS, version 2.12.3. However, when I open it, Symantec says that it has discovered the trojan SONAR.TCP!gen4. See http://securityresponse.symantec.com/security_response/writeup.jsp?docid=2016-010506-2603-99&vid=4294923110.

Is there a trojan in this version of QGIS, or is it a mistake on the part of Symantec?


texts: How to add my own algorithms to ArcGIS toolbox

texts: I want to implement and test some algorithms for line generalization tool in ArcGIS software by Python programming, but I do not have any access to its tool bar script. My question is, how can I add my algorithms into existing ArcGIS algorithms in the line simplify tool?


texts:
I am looking for the Math and steps involved not using a tool. We have PostGIS but its an old version so we are working with the data in Python.

Given three coordinates:

-107.37, 52.46
-103.24, 52.62
-101.09, 53.97

I then join them to form a line.

I then build buffered polygons around each line segment (20NM) - resulting in two rectangles. If there was another coordinate I would have 3 rectangles and so on.

This is what the output looks like when the rectangles are assembled:

However, what I really want is to have that triangular gap - sealed up with a curve.

Like this (did it by hand):

So I have everything working just fine and dandy, but I have no idea how to properly join these rectangles so that they have a nicely fitted curve to fill in that triangular gap.

Any help would be greatly appreciated.

Edit:

Since I am assembling the polygon data myself I can easily identify the points that need to be connected. Is there a way to join two points with a curve?

Thanks!


texts:
I am in the process of writing a Plugin which should react to some data (XML-Format) sent via a socket connection. Since I have never worked with IPC in terms of python before and my knowledge of IPC is, sadly, pretty limited, I'd like to ask you for some ideas and/or help.

So far, my first shot at the problem was to implement a simple non-blocking echo-server inside of QGIS just to have something similar running. The code of the echo-server is basically taken from here and a little modified as follows:

def start_listener(self):
    PORT = 10000
    server_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    server_socket.bind((socket.gethostname(), PORT))
    server_socket.listen(5)
    # Incoming Sockets
    inputs = [server_socket]
    # Outgoing Sockets
    outputs = []

    message_queues = {}
    while 1:
        # "poll" every 0.25 Seconds
        readable, writeable, exceptional = select.select(inputs, outputs, [], 0.25)

        # Sockets ready to be read
        for s in readable:

            if s is server_socket:
                client_socket, client_address = s.accept()
                client_socket.setblocking(0)
                inputs.append(client_socket)
                message_queues[client_socket] = Queue.Queue()
            else:
                data = s.recv(1024)
                # if data received
                if data:
                    message_queues[s].put(data)
                    if s not in outputs:
                        outputs.append(s)
                else:
                    if s in outputs:
                        outputs.remove(s)
                        s.close()
                        del message_queues[s]

        # sockets ready to write
        for s in writeable:
            try:
                next_msg = message_queues[s].get_nowait()
            except Queue.Empty:
                outputs.remove(s)
            else:
                s.send(next_msg)

        # Sockets in Error-State
        for s in exceptional:
            inputs.remove(s)
            if s in outputs:
                outputs.remove(s)
            s.close()
            del message_queues[s]

The Code is launched as a Thread as QGIS will block otherwise. Ultimately the stream contains data which cause QGIS to select given features and pan to those features. My exact question now is: What is the best point of loading such a server? If I'm correct there are two options, namely startup.py or somewhere in the plugin (maybe modify __init__ ?). If I load the server inside the startup.py, is there a way to get a reference to it so that the data can be processed in the plugin-code? Or is there even a better way to solve this problem?

Another problem with the given code is, that it causes a pretty high CPU-load when a connection is handled and the load does not drop after the client has closed the connection.


texts: Python Installation and accessing Arcpy Libraries without ArcGIS Setup Installer?

texts: As part of developing scripts and to save the other services/resources accessed by ArcGIS, is there a way or separate installer to install Python and ArcPy libraries alone without installing ArcMap?


texts: ArcGIS Python - Convert feature class names to lower case

texts: Can anybody anyone point me to directions on how to convert file geodatabase feature class names to lower case by using Python?


texts: pgr_nodeNetwork doesn't place new nodes properly

texts: I wanted to use pgr_nodeNetwork tool to split lines where there's no vertex on intersection. The tool doesn't seem to work as expected, though. It doesn't place new nodes right where the two lines cross. Tolerance value doesn't affect the result. In the first attached image you can see original location of nodes and in the second one you can see the result of pgr_nodeNetwork function. This is pgrouting version 2.0 with PostGIS 2.2.0. Any ideas?


texts: 1:n relations in label expressions & general syntax

texts:
I've got a colleague who wants to create a (concatenated) label based on a 1:n relation between some shapes and related tables. They asked me if this was possible in QGIS - and I think it should.

I found this:
QGIS: relations in label expressions

However, I failed to generalise this for 1:n and m:n relations. The solution provided by @ben-w isn't using the relation, and the solution by @alexander-novikov uses relation_aggregate() which is not really suitable for simply just getting multiple character values our of a relation.

Can someone point me to to right direction, and maybe provide some syntax explanation?

(NB: Sorry not to provide a reproducible example. I not using QGIS expressions normally, I am coming from R using {dplyr}, and I don't speak SQL, nor Python. Different syntax, maybe sometimes different logic. R is no solution for my colleague, so I try to learn something new.)


texts: OpenLayers using GoogleMaps layers: issues with setCenter and transform

texts:
Past week or so I have been doing everything in just OpenLayers and has been going great: been using markers, vectors, bringing in json data via php, and some ajax.

Everything has been in one layer in regards to the LayerSwitcher() control.

Now I would like to add a different layer (so, another layer in the LayerSwitcher control) that uses GoogleMaps as the background. Whether it ends up being hybrid/street/physical is not an issue right now.

When defining layers, many OL examples I see set layers as an array in the map object. This is ok but I would rather just define them separately via var. For example:

var layer1 = new OpenLayers.Layer.MapServer( "My layer that has been working well",
"mapserv.exe url", {map: '/mapfile_dir', layers: ['ex1', 'ex2] }, {singleTile: true});

var layer2 = new OpenLayers.Layer.Google("Google Layer that is troublesome",
{type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20})

The goal:

Up to this point I've been using layer1 with vectors/markers and all has been working well. Everything that is working well in layer1 needs to work exactly the same in layer2.

So basically I'd like the functionality of layer1 (with all the working vectors/markers/ajax/json) to be the same in layer2. Frankly, the only real difference will be the map itself (whether it be Google Hybrid/Physical/Street)

I've been reading about .transform and it seems more fickle than initially anticipated, so I just tried to add the 2 above basic layers, using a setCenter as a test to see if the center would come up the same for both layers.

The code for that simple test:

var map;

function init() {
  map = new OpenLayers.Map('map', {
});
var glayer =  new OpenLayers.Layer.Google("Google Hybrid",
{type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20})
var ilayer = new OpenLayers.Layer.MapServer("Interactive",
"url//mapserv.exe", {map: '/mapfile', layers: ['layershere'] }, {singleTile: true});

map.addLayers([glayer, ilayer]);
var lonLat = new OpenLayers.LonLat(-82, 28).transform(new
OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());
map.setCenter (lonLat, 7);
map.addControl(new OpenLayers.Control.LayerSwitcher());
}

If I do just the glayer by itself, everything is fine: setCenter works. Now when I add in ilayer (as the second layer) it just comes up as a black map when I switch to it.

If I do the ilayer first like so:

map.addLayers([ilayer, glayer]);

The ilayer appears, then when you switch to the glayer the setCenter is on the coast of Africa (0,0 as we all know too well). So this would mean it's not working.

Setting displayProjection: new OpenLayers.Projection("EPSG:4326") in the map object has not accomplished anything.

So, what am I doing wrong? Like I said I'm decided to start from scratch again just to see if I could get a simple setCenter working.

In addition, I hope I don't have to do a transform on all the initial vectors and markers, but we shall see.

Any input on the matter is greatly appreciated.


texts: Enable single feature selection in QGIS using Python

texts:
I want to enable single feature selection using Python. I want to make a script in Python that will first enable single feature selection, and then it will return the id of selected feature. I have coded the second part, but how do I enable single feature selection?

layer = iface.activeLayer()
h = layer.selectedFeatures()
    for f in h:
        g =  f['ID']
        print g


texts: Error with Leaflet Tutorial on adding Mapbox tile layer

texts: I'm trying to add a Mapbox Classic base layer from my account using Leaflet's tutorial. I can't figure out what's wrong with my code:


texts: Using arcpy.mapping.ExportToTIFF() to create dataframe geotiff

texts:
Has anyone created a dataframe geotiff using arcpy.mapping.ExportToTIFF() with satisfactory results?

My output has poor quality and the georeferencing is significantly off.
This example is the geotiff at 50% transparency on top of my data. The shift between "South El Monte" is the best shift example.

Here is an example of the code:

currentMxd = arcpy.mapping.MapDocument(r"F:\Test\amp_Edit.mxd")
dataFrame = arcpy.mapping.ListDataFrames(currentMxd, "Layers")[0]
arcpy.mapping.ExportToTIFF(currentMxd, r"F:\Test\outtif.tif", dataFrame, df_export_width = 46000, df_export_height = 29000, resolution = 1000, geoTIFF_tags = True)


texts: Using the logical operator 'AND' in the query builder

texts:
I am trying to use the query builder to select numerous features in a layer.

I saw a question on here that suggests using the AND operator to do this. however, I have this to no avail.

My query returns no selected roads (see image).

Can anybody help?


texts: How to ensure centroid is calculated within polygon?

texts:
I am working with a Search Cursor approach to calculate centroids within polygon features. You can see from the screenshot that the script calculates centroids within polygon features, however, some of the centroids are located outside of the polygon features (circled in red). Is there a way to force the centroids to be calculated within the polygons? It appears as if the Feature To Point (Data Management) tool allows for this, although I am unable to use an an Advanced license for this.

import arcpy

input_fc = "C:\\temp\\geodatabase.gdb\\featureclass"
output_fc = "C:\\temp\\geodatabase.gdb\\featureclass_centroids"

cursor = arcpy.da.SearchCursor(input_fc, "SHAPE@XY")
centroid_coords = []
for feature in cursor:
    centroid_coords.append(feature[0])

point = arcpy.Point()
pointGeometryList = []

for pt in centroid_coords:
    point.X = pt[0]
    point.Y = pt[1]

    pointGeometry = arcpy.PointGeometry(point)
    pointGeometryList.append(pointGeometry)

arcpy.CopyFeatures_management(pointGeometryList, output_fc)


texts: Why does Points2One complain "Can't make a polyline out of one points."?

texts:
I am trying to connect pairs of points by lines using the Points2One tool in QGIS. I have created the data and exported it as shapefile. The data looks similar to the data in this question:
Workflow for creating line features between two coordinate pairs in QGIS

The problem is after grouping the data by id and sorting by order I get an error which says "Can't make a polyline out of one points." I am not sure what I am doing wrong here as I am new to QGIS.

Here's the attribute table:
id,       order,    longitude,  latitude
stream1,    1,     80.187006,   12.989196
stream1,    2,     78.509647,   17.441657
stream2,    1,     78.454394,   17.242307
stream2,    2,     80.157087,   13.112136


texts:
Nowadays i'm working with Mapcript / Mapserver. I would like to know how to used WFS (Web Feature Service) with ASP.NET.

I think i need to use WCF, right?
Is there any tutorial or framework?

Thanks!


texts:
I Use QGIS ver 2.14.1

I've got a little problem here
I already following how to open mdb at link below :

QGIS 2.4 access .mdb add error

But error dialog box always pop out when I click "OK" (after I Create connection and test it)
It said "Connection Failed - Check settings and try again."
"Extended error information : ORA-12541: TNS no listener in "

Why error dialog box always pop out when I test the connection ?


texts:
I am trying to transform a geometry to the EPSG:4978 projection. As the EPSG:4978 was missing in the spatial_ref_sys table, I have added it using the following query:

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 4978, 'epsg', 4978, '', 'GEOCCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["m",1.0],AXIS["Geocentric X",OTHER],AXIS["Geocentric Y",EAST],AXIS["Geocentric Z",NORTH],AUTHORITY["EPSG","4978"]]');

Source: http://spatialreference.org/ref/epsg/wgs-84-2/postgis/

When I however try to use the projection in one of my queries, it fails with:

ERROR: AddToPROJ4SRSCache: could not parse proj4 string ''

With Proj4js the empty string definition works just fine, so it seems to me that Postgis should handle the empty string case as well. Or am I misunderstanding something specifically about EPSG:4978? What is the proper way to transform a coordinate to EPSG:4978?


texts: Dataset does not exist or is unsupported SaveLayerToFile

texts:
I am trying to save a layer from an MXD to a lyrx file in a geodatabase using ArcGIS Pro.

aprx.importDocument(current_document_name)
mlist = aprx.listMaps()
infra_layers = mlist[0].listLayers()
arcpy.SaveToLayerFile_management(infra_layers[0].name, new_file_name, 'ABSOLUTE')

I am getting the following error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Layer: Dataset active_infrastructure_projects does not exist or is not supported
Failed to execute (SaveToLayerFile).

So the name exists and the layer is being imported.

Why can't it save the layer?


texts: How to setup GeoWebCache to produce same grid as GoogleMaps?

texts:
Is it possible to create gridSet in GeoWebCache which has same number of zoomlevels with scaleDenominators compatible with google maps?

There is EPSG:4326 default, but no mention of number of zoomlevels/scaleDenominators in documentation.


texts: Implement 'Zoom To Selected Features' in VB.Net

texts: Seems like this should be dead simple.  I'm writing an add-in for ArcMap 10 in VB.Net.  I need code that will reproduce the 'Selection --> Zoom To Selected Feature' menu option.


texts:
Is there any easy way to compare if features in one set are sharing vertices?

If I am cutting parcels and I want the next parcel over to have the same vertices?

I converted my vertices to points, and then used find identical to determine the mismatch.  Something similar to this but not as cumbersome.


texts: Is disabling CORS on map tiles service security risk?

texts:
I would like to know what security issues or security risk would prevent disabling Cross-Origin Resource Sharing (CORS) on map tiles service.

Most of public tiles services, such OSM, Google, Bing or ESRI, have CORS enabled.

What could be the security risk to disable CORS for a map tiles service?


texts: Save GM directions as KML - with timestamps?

texts:
In Google Maps I use Get Directions and after saving the result as a 'map' and opening it in My Places I am able to save it as a KML. All fine up to a point, allowing me to display the route in Google Earth.

But I'd also like to animate it with GE's Time Slider. Given that the GM directions assumes speeds and provides an estimated journey time, is it possible to get GM to save the KML with timestamps please?

Otherwise I must edit the 'raw' KML elsewhere and add arbitrary times myself, a relatively tedious task.

--
Terry, East Grinstead, UK


texts: Spatial Join with start and endpoints of polylines in QGIS?

texts:
I want to join 2 layers. One of them is are polylines and has two pairs of lat/long, representing the start and the end.

id    lat_start    long_start    lat_end    lat_start

The other layer is a polygon with districts.

district
district1
district2
...

What I want to do is performing a spatial join two times:

with a district_start
with district_end

I want my final line-layer to look like this

id    lat_start    long_start    lat_end    lat_start    district_start    district_end

Of course, I already tried vector --> data management tools --> spatial join and MMQGIS but they don't have the result I want. Actually they even assign districts that never work with the start nor the end.


texts: Displaying link title instead of URL in Infowindow of CartoDB?

texts:
In my data, I have a column containing some url values; where the url is unavailable, it contains the text value "Not yet available".

I want the output in the map infowindow to either display, "Not yet available" in cases where there is no url or "Link to website," where you can click through to the url.

I played with both the sql and html to try to make this happened, but was unsuccessful.


texts: QGIS field calculator and Postgres permissions

texts:
I’ve given user x the following

GRANT ALL ON ALL TABLES IN SCHEMA y TO x;
GRANT ALL ON ALL SEQUENCES IN SCHEMA y TO x;
GRANT ALL ON ALL FUNCTIONS IN SCHEMA y TO x;

However in QGIS they dont have access to check/tick the 'update existing field' button. They can chose a field to update and write an expression for its update e.g $area but the OK button remains greyed out. If I as administrator load the layer, the check/tick box appears and I can write and run expressions.

I'm sure its a permissions thing in PostGIS.

Si


texts: How can I change what port ArcGIS Server listens on?

texts:
I have ArcGIS Server installed with a web adapter on one server. It uses port 80.

I want to install another ArcGIS Server on a different server, but I want it to listen on a different port other than 80, such as 8080 or something.

How would I do this?


texts: Ranking values by shape area using ArcPy

texts:
I'm trying to make a script that performs the following: it creates a new field that will hold values based on size of the shape area (from 1 to end of table), basically ranking them.

I have a few ideas to do it, and managed to write a script, however it is not working, no errors, nothing.

import arcpy

fc = "D:/M1 Geomatique/Programmation II/Dossier/zones-de-rattachement-des-bureaux-de-vote-en-2014.shp"

try:
    arcpy.AddField_management(fc, "RANG", "INT")
    size_rank = 1
    rows = arcpy.UpdateCursor(fc, ["shape_area"], "D")
    for row in rows:
        row.setValue("RANG", size_rank)
        pop_rank += 1
        rows.updateRow(row)

except:
    arcpy.GetMessages()

Another idea would be to somehow adapt this code, but in any case it is mentioned that shapefiles do not support ORDER BY as used above by arcpy.da.UpdateCursor’s sql_clause argument.

One other idea, to use a list (ListFeatureClasses probably) to get all values, sort them by size, make a loop that takes each value and places a number based on their size. This one seems a bit complicated and I don't have a clear image on how to get through with it.

What would you recommend me to do?

I've made the modifications you suggested but it's still not working, not getting anything. Quite odd. I even made some additional modifications because it's not even creating my field:

import arcpy

arcpy.env.workspace = "D:/M1 Geomatique/Programmation II/Dossier"

try:
    arcpy.AddField_management("zones-de-rattachement-des-bureaux-de-vote-en-2014", "RANG", "INT")
    size_rank = 1
    rows = arcpy.UpdateCursor(fc, sort_fields="shape_area D")
    for row in rows:
        row.setValue("RANG", size_rank)
        size_rank += 1
        rows.updateRow(row)

except:
    arcpy.GetMessages()


texts: I'm looking to substitute my own dataset with OSM data in Leaflet Routing Machine or another routing tool. Is this possible?


texts:
How can I add an external online map in ArcMap by inserting the http address, to use as a basemap in order to georeference multiple aerial photos?

I don't really care about the precision of the georeferencing. I just want to use the template as a quick-to-use index for the aerial photos.


texts:
First things first, I am completely new to this so please bare with me if this is a bit unclear as I lack the vocabulary.
I have two raster layers (a DEM and a DSM). The DEM layer covers a large rectangular region while the DSM only covers a part of the same region (the Urban part) and the shape of the area is a complicated polygon.
I want to produce a single raster which will have all the information from both DEM and DSM. When I try to merge those two, I get a valid result (everything aligns nicely) however, it looks like the DSM is being projected on the DEM and this results in gaps.
I have high hopes that the attached image will make the situation clear. Below you can see the result of merging the DEM and DSM rasters in 3D View (top) and 2D view (bottom).
PS: I tried to use GRASS r.patch but unfortunately its UI doesn't allow to chose the order of the layers in QGIS 3.10.


texts:
I created a function to search for php to postgis and believe me that this function a file with the coordinates in the search GeoJSON format.

The problem is that when I do rendering the map and add the table when the query again takes this layer is not updated, you need to update your browser to see the change.

There any way to make it automatically change or refresh the map every time a query is made?

Here 's the code:

        var Vector= new ol.layer.Vector({
              source: new ol.source.Vector({
                projection: 'EPSG:3857',
                url: 'coordenadas.GeoJSON', // File created in php code
                format: new ol.format.GeoJSON()
              }),
              style:  new ol.style.Style({
             stroke: new ol.style.Stroke({
             color: 'red',
             lineDash:[1],
             width: 1
            })})

        });

            var map = new ol.Map({
            target: 'map',
            layers: [
                new ol.layer.Tile({
                    source: new ol.source.BingMaps({
                    imagerySet: 'Aerial',
                    key:'XXXXXXXXXXXXXX'

            ],
            view: new ol.View({
            center : ol.proj.transform([-0.33126,39.569658],'EPSG:4326','EPSG:3857'),
            zoom: 13
            })

        });
        map.renderSync(5000);

the map.renderSync () function; not how it works, I find documentation to understand its operation.


texts: How to measure sinuosity in ArcGIS Desktop?

texts:
I am trying to quantify the sinuosity of a GPS track (shows a loop pattern) and cannot seem to find any way to do this. I need to be able to quantify sinuosity for the removal of some unwanted data.

Is there a tool to measure sinuosity in arcmap, i am very unfamiliar with the codes in ArcMap 10?