texts: Are there any addresses that straddle two or more congressional districts?

texts:
I'm using a shapefile of TIGER data provided by the Census Bureau for locating an individual's congressional district. In PostGIS I'm running an ST_Intersects query in order to take a (lat, long) coordinate and return the results set.

Of course, there are many points along the polygons that intersect both, but are not contained in either, hence the choice for ST_Intersects over ST_Contains.

Are there any real world examples of addresses that actually straddle two congressional districts? This edge case makes me suddenly nervous that I could be presenting users with the wrong district.

EDIT:

It appears that congressional districts might use census blocks which divide streets down the middle, so you technically should not have an address that bounds both. A follow-on question would be whether there are any known conflicts between a service's address geocoding and census mapping.


texts: Calculating nearest distance between polygons according to attribute field in ArcGIS Desktop?

texts:
I have a two polygon shapefiles, one shapefile is a zoning and the second one is a census tract. Now I want to calculate distance from boundary of tract to closest boundary of each zone type (R, C, etc). Because of that I have dissolved all zones according to zone type and now I want to have in each census tract information about nearest distance to the all zone types (distance to C, distance to R, etc...). I am using ArcGIS 10.4 and have trid to use NEAR tool,  spatial join but I can't get a distances to all zone types,  just to the nearest one.

Also I want to have a zero distance value if each tract contain the zone type or intersect it.


texts: Georeferencing terrestrial LiDAR point cloud?

texts:
Recently I took a Faro Focus terrestrial laser scanner to a forested area and collected multiple scans for 12 fixed radius plots. Each plot had 5 GPS points collected within it in hopes to georeference the point cloud. However, I am struggling to take that geographic data and transform the unreferenced point cloud with it.

The goal was to use it alongside aerial LiDAR collected in the same area that is georeferenced with UTM coordinates.


texts: Improving display efficiency of OpenLayers 3 Cluster?

texts:
In OpenLayers3, when the number of points more than 1,000,000, the display efficiency is too low. How to solve this problem? The following code is the example in OpenLayer3, and I had changed the number of points to 1,000,000.

<script>
  var distance = document.getElementById('distance');
  var count = 1000000;
  var features = new Array(count);
  var e = 4500000;
  for (var i = 0; i < count; ++i) {
    var coordinates = [2 * e * Math.random() - e, 2 * e * Math.random() - e];
    features[i] = new ol.Feature(new ol.geom.Point(coordinates));
  }

  var source = new ol.source.Vector({
    features: features
  });

  var clusterSource = new ol.source.Cluster({
    distance: parseInt(distance.value, 10),
    source: source
  });

  var styleCache = {};
  var clusters = new ol.layer.Vector({
    source: clusterSource,
    style: function(feature) {
      var size = feature.get('features').length;
      var style = styleCache[size];
      if (!style) {
        style = new ol.style.Style({
          image: new ol.style.Circle({
            radius: 10,
            stroke: new ol.style.Stroke({
              color: '#fff'
            }),
            fill: new ol.style.Fill({
              color: '#3399CC'
            })
          }),
          text: new ol.style.Text({
            text: size.toString(),
            fill: new ol.style.Fill({
              color: '#fff'
            })
          })
        });
        styleCache[size] = style;
      }
      return style;
    }
  });

  var raster = new ol.layer.Tile({
    source: new ol.source.OSM()
  });

  var map = new ol.Map({
    layers: [raster, clusters],
    target: 'map',
    view: new ol.View({
      center: [0, 0],
      zoom: 2
    })
  });

  distance.addEventListener('input', function() {
    clusterSource.setDistance(parseInt(distance.value, 10));
  });
</script>


texts: Forcing maximum value to raster using QGIS?

texts:
I am trying to force a maximum value for a raster.  Currently the raster has values ranging from 0 to 9, I would like to force values greater than 7 to 7.

I have tried a couple ways of doing this in the raster calculator and am not having any luck.  I've also tried sieving values greater than 7 to a new raster and subtracting that new raster from the existing raster, but that returned values of 7 across the whole raster.


texts: Check validity produces error output at origin although no polygon is there QGIS

texts:
After running 'Check Validity' I keep getting an 'error output' dot at the origin for my shapefile although there is no feature there.

This is the only error, and I think may be why my intersection is not working.
Even after putting my file through mapshaper to try to simplify it, I keep getting this error output point at the origin.

I am not sure how to fix this problem, or if it even is a problem, perhaps it is just a normal artifact of the process and my error is elsewhere. Could it be that I could have a feature without a geometry at this point?

The shapefile was originally in equal area projection, then re-projected to WGS84, however I do not think this is related.


texts:
I want to use GetFeatureInfo in WMS layer.  I have got some code but it shows only a blank window, and there is not any error.

This is the code link https://jsfiddle.net/seema30/LtcvhLwf/.


texts:
I have a rivers layer which I am trying to get "Curved labels" on. When I zoom in the layer is labelled correctly with curved labels. However, when I zoom out the curved labels disappear completely.

I know other people have posted this question before and was wondering if anyone had found a solution to the problem ?.

Parallel labels works fine at all scales. I've also attempted playing with the 'Max angle between curved charachter setting' but this still dosn't solve the issue.


texts: ArcGIS Field Calculator for changeable 3 fields add from table parameter at model builder and Python

texts:
I have an Excel sheet enter it to ModelBuilder as a parameter and at next steps i need to create an expression for the fields inside the table and add the results to another field.

The problem is the fields not exist to select for expression at field calculator.

So is there any way to solve this problem and how could i write ^ inside expression.  I have ArcGIS 10.5

I tried to create a python code but it give me an error Invalid field..
I wish to know what is the problem inside the code

Table name: #it will be parameter

Fields at table need to enter at expression: Geco, Par, kela (All numbers i consider all as float)
The new field is MM

The code i tried with is:

Geco = arcpy.GetParameterAsFloat()
Par = arcpy.GetParameterAsFloat()
kela = arcpy.GetParameterAsFloat()

Geco_new = '!'+Geco+'!'
Par_new = '!'+Par+'!'
kela_new = '!'+kela+'!'

arcpy.CalculateField_management(outputfc, percentfield, '('+Geco_new+'*(2/'+Par_new+')^'+kela_new+'))', "PYTHON")


texts: Python script to batch stack Landsat images

texts:
We are trying to make a script that goes through a folder full of unzipped Landsat folders and stacks all the Landsat images it finds. However since Landsat 8 images have more bands then we want to re-order the stacking so that the thermal bands are at the end of the composite stack. The script works for just Landsat 8 images, but we tried to add in a line so that if the number of TIFs in a folder is less than 10 then it just stacks them normally and doesn't re-order them.

The script print out is as follows, which seems to suggest there is a problem with the lists for the non-Landsat 8 images having lists within them.

Stacking Landsat bands together
[u'LC80430202016124LGN00_B1.TIF', u'LC80430202016124LGN00_B2.TIF', u'LC80430202016124LGN00_B3.TIF', u'LC80430202016124LGN00_B4.TIF', u'LC80430202016124LGN00_B5.TIF', u'LC80430202016124LGN00_B6.TIF', u'LC80430202016124LGN00_B7.TIF', u'LC80430202016124LGN00_B8.TIF', u'LC80430202016124LGN00_B10.TIF', u'LC80430202016124LGN00_B11.TIF']
Stacking bands for image D:\Landsat\New\LC80430202016124LGN00.tif
[u'LC81420402016034LGN00_B1.TIF', u'LC81420402016034LGN00_B2.TIF', u'LC81420402016034LGN00_B3.TIF', u'LC81420402016034LGN00_B4.TIF', u'LC81420402016034LGN00_B5.TIF', u'LC81420402016034LGN00_B6.TIF', u'LC81420402016034LGN00_B7.TIF', u'LC81420402016034LGN00_B8.TIF', u'LC81420402016034LGN00_B10.TIF', u'LC81420402016034LGN00_B11.TIF']
Stacking bands for image D:\Landsat\New\LC81420402016034LGN00.tif
[[u'LE71420402002339SGS00_B1.TIF', u'LE71420402002339SGS00_B2.TIF', u'LE71420402002339SGS00_B3.TIF', u'LE71420402002339SGS00_B4.TIF', u'LE71420402002339SGS00_B5.TIF', u'LE71420402002339SGS00_B6_VCID_1.TIF', u'LE71420402002339SGS00_B6_VCID_2.TIF', u'LE71420402002339SGS00_B7.TIF', u'LE71420402002339SGS00_B8.TIF']]
Stacking bands for image D:\Landsat\New\LE71420402002339SGS00.tif
[[u'LM21520401977080AAA05_B4.TIF', u'LM21520401977080AAA05_B5.TIF', u'LM21520401977080AAA05_B6.TIF', u'LM21520401977080AAA05_B7.TIF']]
Stacking bands for image D:\Landsat\New\LM21520401977080AAA05.tif
[[u'LT51410402010362BKT00_B1.TIF', u'LT51410402010362BKT00_B2.TIF', u'LT51410402010362BKT00_B3.TIF', u'LT51410402010362BKT00_B4.TIF', u'LT51410402010362BKT00_B5.TIF', u'LT51410402010362BKT00_B6.TIF', u'LT51410402010362BKT00_B7.TIF']]
Stacking bands for image D:\Landsat\New\LT51410402010362BKT00.tif
Finished!

The script follows

import arcpy
import os
outws = r"D:\Landsat\New"
arcpy.env.workspace = outws

folders = arcpy.ListWorkspaces()
newrasters_l8 = []
newrasters_l7 = []

for folder in folders:
    arcpy.env.workspace = folder
    rasters = arcpy.ListRasters("*.tif")
    if len(rasters) == 12:
        newrasters_l8.append(rasters[0])
        newrasters_l8.extend(rasters[3:10])
        newrasters_l8.extend(rasters[1:3])
        name_l8 = os.path.join(outws, rasters[1].split("_")[0] + ".tif")
        arcpy.CompositeBands_management(newrasters_l8, name_l8)
    elif len(rasters) < 10:
        newrasters_l7.append(rasters[0:])
        print newrasters_l7
        name = os.path.join(outws, rasters[1].split("_")[0] + ".tif")
        print "Stacking bands for image "+name
        arcpy.CompositeBands_management(newrasters_l7, name)


texts: Tile Package being created with all zero values?

texts:
I am attempting to create a tile package of our county aerial imagery to use in a third party product. I have tried running it a few different times now with a few different tiling schemes and each time It comes out with zero values in each band for all cells in the raster, even though it takes over 10 hours to run and has a large file size (see attached image).  However, this only happens if I include larger scales such as 1:4000, 1:2000, 1:1000 (like in the picture of the tiling scheme I attached). If I only use smaller scales like 1:100,000 then the tile package is created as normal and looks fine.

Also, I have tried creating the tile package from both a .tif as well as a .sid mosaic.

Here is my exact work flow:
1. Run "Generate tile cache tiling scheme" to create a tiling scheme to use
2. File>Share as Tile Package (selecting the tiling scheme that I created in step one)

Has anyone run into similar issues or have any idea what the problem is?

FYI I am running ArcMap 10.5.1

Tiling Scheme


texts: Keeping a GeoServer image mosaic layer current as new geoTIFFs are added to the filesystem

texts:
I have an ImageMosaic layer with the time dimension enabled in GeoServer 2.8.x using postgres/postgis as the index store.  The system is running CentOS 7, and both GeoServer (jetty variant) and the database are on the same host.  The directory containing the geoTIFFs for the layer is kept current by rsync running on the server producing the geoTIFFs, and we are keeping the last 10 images.  I have been trying to find a way to keep the index info in the database current as images are added to/deleted from the directory.

Using the REST API seems like an option, and maybe one of, or the only option to do so.  While testing, I was able to delete all indexes from the database table using the DELETE HTTP method, but not recreate them; reharvest?  I think my biggest hurdle is just understanding the REST API.

Right now, at least, the less elegant method of recreating the entire granule index is fine.  Since I think some of the non-specific to my setup info is throwing me from the posts I found... I'm working with a layer called baserefl-loop.  I was able to list the granules using:

curl -v -u user:pass -XGET http://localhost:8180/geoserver/rest/workspaces/Sonalysts/coveragestores/ \
radar-baserefl-loop/coverages/baserefl-loop/index/granules.xml

Hopefully this is just something simple that I'm missing.


texts: VB Script case statment syntax for the field calculator in ArcGIS 10.2?

texts:
Can anybody help me with this VB Script code for the field calculator in ArcGIS 10.2?
Why in the final result are all the values 10?

a= 10
Select Case UCase( [PRIORITY] )
 Case "AUX 1": a = 22.5
 Case "AUX2": a = 25
 Case "PRIMARY": a = 17.5
 Case "AUX 2": a = 25
 Case "AUX 3": a = 30
 Case "AUX 4": a = 32.5
 Case "NA" : a = 5
 Case "N/A": a = 5
 Case "SECONDARY": a = 15
 Case "SECONDARY2": a = 15
 Case "Tertiary": a = 10
 Case "MALL": a = 27.5
End Select


texts: Lizmap add point to PostGIS layer - correct layer settings?

texts:
Not able to add point to the layer using lizmap.
I'w been following this guidance
https://github.com/3liz/lizmap-documentation/blob/master/source/publish/advanced_lizmap_config.rst#id28

Point is not saved & lizmap page says:
"An error has been raised when saving the form!"

Lizmap error log says:
.... 2018-04-24 12:51:46    [2]   openlog() expects parameter 3 to be integer, string given ....

Lizmap messages log says:
2018-04-24 12:51:40 ::1 default   SQL =  INSERT INTO "public"."mytablename" ("geom", "Name", "description" ) VALUES (ST_GeomFromText('POINT(21.785181524973595 63.10872517176041)', 4326), 'blaablaa', 'blaaa' )  RETURNING "id";

My setup is Win10
Lizmap 3.2Beta + QGIS Server 2.18.18 + Postgres 10.2 & PostGis
plugin Lizmap 2.4.0

Layer "id" field i set to hidden using QGIS layer properties. + tick on WMS & WFS. And "id" type is integer and type name int4

By using PGAdmin i'm able to see that this id field has a Default value defined as nextval('mytablename_id_seq'::regclass). this field is a primary key.

On Lizmap plugin this layer is added to be edited with all the options +
Popup enbabled + qgis setting


texts: Saving last selected value in Python Script Tool?

texts:
I have created a python script tool with python 2.6 and arcgis 10.0 in which the user parameter reads a field and provide list of field content , below is the list.

user selects the desired value and run the tool.

now when the user again open the same tool he has to again refine the list to choose desired value.

i want this tool to store previously selected option from the list, so it opens like this: .

how can i achieve this ?


texts: How to get the center and radius of the feature in openlayers

texts:
I created a control to draw circle and rect to the map:

        var supportControls = {
            darwpoint : new OpenLayers.Control.DrawFeature(this.drawLayer, OpenLayers.Handler.Point),
            drawline : new OpenLayers.Control.DrawFeature(this.drawLayer, OpenLayers.Handler.Path),
            drawpolygon : new OpenLayers.Control.DrawFeature(this.drawLayer, OpenLayers.Handler.Polygon),
            drawrect : new OpenLayers.Control.DrawFeature(this.drawLayer, OpenLayers.Handler.RegularPolygon, {
                handlerOptions : {
                    sides : 4,
                    irregular : true
                }
            }),
            drawcircle : new OpenLayers.Control.DrawFeature(this.drawLayer, OpenLayers.Handler.RegularPolygon, {
                handlerOptions : {
                    sides : 40
                }
            })
        };
        for ( var key in supportControls) {
            supportControls[key].handler.stopDown = true;
            supportControls[key].handler.stopUp = true;
            supportControls[key].featureAdd=function(feature){
                handlerFeature(feature);
            }
            map.addControl(supportControls[key]);
        }

        function handlerFeature(fea){
            // how to get the circle center and radius of the feature if it is a circle
        }

Now, after the draw end, I have have no idea to differentiate whether the feature is circle or rect, and how to get the center and radius if it is a circle?


texts:
I have a set of tiles for a certain area, I have to load them in the client, and draw some points,lines on the tiles.

Generally there are some opensource android library like osmdroid, mapsforge and etc to load the tiles , however the tiles I am working with are not generated by the common projection like EPSG4326 or EPSG3857.

They are generated by Arcgis Server, with a custom coordination system, I can get the tile origin, resolution of each zoom, I can load the tiles by openlayers(because ol support custom resolutions) but I can not load them by the popular library since none of them support the custom srs.

At the moment, the only choice is the Arcgis Android SDK, while it is huge, with almost 30M which I tried to avoid.

Is there other alternative solution?


texts: Added raster layers with PyQGIS do not show up

texts:
I have a list ("datList") with names of raster layers which I want to add. When I add raster layers in the Python Console everything works fine:
for da in datListe:
   layerDir = "C:/[path]/" + da
   iface.addRasterLayer(layerDir,da)

When I do the same in a processing script with...
for da in datListe:
   layerDir = "C:/[path]/" + da
   rlayer = QgsRasterLayer(layerDir,da)
   QgsProject.instance().addMapLayer(rlayer)

the layers are not added to my map. However, they seem to be in the project as
QgsProject.instance().mapLayers()

(typed in the Python Console) prints a list with the desired layers


texts:
I have a 50 MB shapefile. I am trying to import it into Illustrator. When I try printing to an svg file from the Print Composer, QGIS crashes. I am using QGIS Essen 2.14.3. I have tried lots of different things:

simplifying the geometry (QGIS crashes)
trying to print to a postscript file (will only print to a prn file, which Illustrator can't read)
installing Cute PDF (doesn't really work in Illustrator)

I am using a laptop with Windows 10, 32 bit and 4gb of RAM.

Will upgrading to a new machine with 8gb of RAM and 64 bit windows 10 make any difference?

I can't afford a new laptop, but I need to output to Illustrator so will probably have to bite the bullet.


texts: Create WFS with Mapguide Studio Enterprise 2011

texts:
i need create a WFS with MGS 2011, i used the help users mapguide for create the WFS but that didnt work, can anyone help me please?
i am new in this tecnology.


texts: Cannot decode CRS from executable JAR with GeoTools

texts:
I decode a CRS from its EPSG code using:

CoordinateReferenceSystem crs = CRS.decode("EPSG:3035");

This works like a charm, but when the same line of code is executed from an executable *.jar program, the following error occurs:

  Exception in thread "main" org.opengis.referencing.NoSuchAuthorityCodeException: No code "EPSG:3035" for authority "EPSG" was found for an object of type "EngineeringCRS".

It seems the EPSG codes can no longer be decoded. This is confirmed by:

System.out.println(CRS.getSupportedAuthorities(true));
System.out.println(CRS.getSupportedCodes("EPSG").size());

which returns

  [EPSG, http://www.opengis.net/gml, http://www.opengis.net/def, AUTO2, AUTO, urn:ogc:def, urn:x-ogc:def, CRS]

  6741

from the original code, but:

  [EPSG, http://www.opengis.net/gml, AUTO2, AUTO, urn:ogc:def, urn:x-ogc:def, CRS, http://www.opengis.net/def]

  2

once executed from the *.jar program.

I am stuck for hours with this issue...

Changing the dependency from gt-epsg-hsql to gt-epsg-wkt did not fix the problem. When unzipping the *.jar file, the EPSG database is there, but it seems it cannot be read.

The *.jar program is produced using the maven-assembly-plugin.


texts: QGIS / Spatialite / DB Browser : basic : joining tables and more

texts:
Question : Is it possible to create a table in QGIS that will be a dynamic reference for geographical data entries?

As I'm sure this will sound very general and unclear, here's the actual situation I would like to address :

let's say I have a vector (points) map of trees on a site.
Each tree is identified by a unique number, and is further identified by a code which indicates its species (to the cultivar level).
Each species has their own characteristics, such as environmental preferences (soil, sun, climate, etc.), pests, uses (fruits, timber, wildlife habitat, etc.) behavior (size, root system, blooming time, etc.), and much more

-> My goal : create a table which lists all these characteristics by species, so that when I add a tree on my vector layer, I only enter the unique number (or let QGIS do it for me), and the species code, and it dynamically retrieves all the info from the species list, and includes it in the vector layer's attributes for that tree. Of course, I'd use the same field names and type in both tables.
And in reverse, let's say I modify some characteristics for a species, and it automatically updates all the unique trees in my vector layer.

I've experimented with many things in QGIS, and one thing that bugs me is that nothing seems to "stay" or be dynamic. For instance, to make QGIS calculate a polygone surface for an 'Area' field, I need to redo the formula every time I add new entities. Can't it be dynamic?

Context :

I'm laying out the needs and solutions to map a future farm on QGIS. Not just trees but shrubs, perennials, annual crops, animals, and much more...

I have a moderate (fortunately fruitful) experience with mapping with QGIS, but it has been 100% self-taught and I'm limited by the database aspect of it, now that I want to push it further.
The last time I really dealt with databases was on Access at university, about 10 years ago, and it wasn't super deep.

I just downloaded DB Browser for SQLite (I'm on OS X), which gives me the pleasant ability to do simple edits that QGIS won't let me, such as changing the order of fields. It also seems to give me the ability to quickly mess up my data, which is why I decided to join Stackexchange and go beyond self-learning.
I'm having a tough time analyzing what is feasible in QGIS itself or requires external handling of the DB.


texts:
I have this orange Polygon, which is not saved in any layer.

How can I get rid of this and where does it come from ?

I am new to GIS.


texts:
In previous QGIS version, after load whatever layer, Data Source Manager automatically closed. However, in QGIS 3, as it can be observed at following image, Data Source Manager remains open afterward load, e.g., a raster layer. It is necessary to close it manually.

I could not find any useful in Settings Options Menu. How can QGIS 3 Data Source Manager automatically close after load a layer?

Editing Note:

I can add multiple layers at once in this way (by using previous Open GDAL/OGR supported dataset(s)):

I closed Data Source Manager quickly with Esc key. So, best option is given by @Steven Kay.


texts: Distance of points to a polygon, inside and out, without Proximity toolset, ArcGIS 10.3.1 Desktop

texts:
I have a very large data set of points representing residential properties and i need to calculate the closest distance to the boundary of one polygon from each point. Preferably a positive value for those outside the polygon and negative for those within. (Example below.)

I have been looking through various other questions similar to this but most users either have access to QGIS or the Proximity toolset within Arc. I have access to neither of these but am hoping there is an alternative way of calculating this distance?
I really do not want to sit and measure out the distance from each point manually, there are over 700 points!


texts:
I have conflict data from ACLED. The data is in CSV format. Given the nature of the data (image attached), when I am selecting the CSV option while importing the file in QGIS via "open data manager", "X" and "Y" fields do not appear as longitude and latitude respectively.
The heading of the first column in the CSV file rather appears in the "X" and "Y" fields. However, if I select the option Custom delimiters and then semicolon (I deselect the option tab), then the problem gets sorted out. "X" field is then longitude and "Y" is latitude
.
Am I doing the right thing? I want to avoid problems appearing later.


texts: I have a .adf file and want to select an elevation range

texts:
I'm super new to GIS.

I downloaded a USGS 30m resolution map of a bay, and now I want to specify all areas that are between 0m and -30m elevation. I need to find where a defined area contains 0 and -30m elevation.

The file is .adf, do I need to convert it to a .shp? I went into symbology and added all values, so I can see the sections that are between 0 and -30m, but I don't know how to select them and do anything with them.

I have ArcMap 10.2.2 for Desktop.


texts:
I'm trying to extract climate variables in QGIS for sample points with a raster dataset. The raster comes from ClimateNA (https://adaptwest.databasin.org/pages/adaptwest-climatena) and is available as in ASCII (.asc) and NetCDF (.nc) formats.

When I load either file format to QGIS, I get "CRS was undefined: defaulting to CRS EPSG:4326 - WGS 84," which is my project CRS. The raster is visible in my Layers tab but not on the map, and when I use the Point Sampling Tool to extract data for my points, the output CSV lists all my sample points but has an empty column for the variable.

Does anyone know what I need to do to fix this?


texts: Way I can find out when a shapefile was created or last updated

texts: I'm working with 100+ layers on QGIS, and I am trying to determine when each layer was last created or updated. All the layers are shapefiles. Is there a way I can find the creation/update date?


texts: Comparing areas of overlapping shapes in QGIS

texts: Given SHP files for a country's Districts in Year1 and Year2, how to output a table whose rows state the percentage composition of each Year2 District in terms of portions of Year1 Districts


texts: Choosing an aggregation strategy that will be spatially represented on a google map

texts:
I'm under an NDA so I cannot really get into specifics of what I am doing.  But, suffice to say, I am taking a large amount of spatially-acquired data and drawing this data onto a Google Map surface using color to represent a range of values.

The main issue we are having is due to the vast amount of data collected, we cannot render this data on a map in a timely manner.

I have proposed that one way to improve performance is to aggregate the data in some fashion, thus reducing the number of data points on the map.

Spatially aggregating the map points are pretty easy, but the data is a whole other matter.  We've tossed around simple statistics (min, max, median, mode) as well as some "made up" statistics (you don't want to know).

We were going to start with a simple average, but edge effects are causing some grief as the user zooms in/out of the map.

Generally, we need to maintain trends in a given area as one zooms out on the map.  In other words, areas of the map with high correlation will look proportionally similar to other correlated areas at all zoom levels.

EDIT

Imagine a rigid body with up to 64 spatially oriented points where each point has some numeric value.  At some interval (typically 200ms), a snapshot of this rigid body takes place and we record each point and map it.  This rigid body moves in a predictable manner with some variable speed.  It is possible in a given square mile that we have tens of millions of points.

Spatially, correlation only applies if that point "belongs" to the same pass of data.  So for instance, if this rigid body moves North for several meters, then makes a u-turn and moves South along side it's previous path, those points that belong to the North path are spatially "different" than those that belong to the South path even if a point on the North path are within inches of a point on the South path.  That might be a long-winded explanation for saying we have spatial/temporal correlation.

I can't really describe in great detail what the purpose of the map is due to the NDA.  Essentially, the user will look at this map and understand how good or bad things went.  They can see this data on a Google map from zoom level 15 and up.  Rendering tens of millions of points at zoom level 15 is painfully slow.

What I am looking for is a strategy to keep the same "message" at each zoom level while achieving a good performance.


texts: Broken packages and unmet dependency installation QGIS, PostgreSQL and postgis Ubuntu 14.04

texts:
I am aware that this could be a duplicate topic but I am struggling for a several days.
I upgraded QGIS from 2.8 to 2.18 and after that everything get broken. I couldn't use PostGIS extension on pgadmin3.

I figured out that there were a conflict with a libgdal1h: if installed, it removes other libraries used by PostGIS and vice versa.

I tried several tutorial about how to install QGIS, like for example:

How do I get and install the latest version of QGIS in Ubuntu?

Install QGIS 2.14 Essen on Ubuntu 16.04 Xenial

I tried every suggestion reported on line, like for example:

https://askubuntu.com/questions/351899/qgis-quantum-gis-install-fails-unmet-dependencies

https://askubuntu.com/questions/621406/dist-upgrade-in-ubuntu-14-04-uninstalled-qgis-and-results-in-unmet-dependencies

Where to find a version of libgeos-c1 compatible with libgeos-3.5.0?

QGIS install on Ubuntu 14.04 fails

Always, the reply of sudo apt-get install qgis python-qgis is:

The following packages have unmet dependencies:

python-qgis : Depends: python-qt4-sql but it is not installable
               Depends: python-qgis-common (= 1:2.14.8+20trusty-ubuntugis) but it is not going to be installed
               Depends: python-psycopg2 but it is not installable
               Depends: python-qscintilla2 but it is not installable
               Depends: python-jinja2 but it is not installable
               Depends: python-markupsafe but it is not installable
               Depends: python-pygments but it is not installable
               Depends: python-yaml but it is not installable
               Depends: python-pyspatialite but it is not installable
               Depends: libqgispython2.14.8 but it is not going to be installed
               Depends: libqgis-analysis2.14.8 but it is not going to be installed
               Depends: libqgis-core2.14.8 but it is not going to be installed
               Depends: libqgis-gui2.14.8 but it is not going to be installed
               Depends: libqgis-networkanalysis2.14.8 but it is not going to be installed
               Depends: libqgis-server2.14.8 but it is not going to be installed
               Recommends: liblwgeom-dev but it is not installable
 qgis : Depends: gdal-abi-2-1-0 but it is not installable
        Depends: libgdal20 (>= 2.0.1) but it is not installable
        Depends: libgeos-c1v5 (>= 3.4.2) but it is not installable
        Depends: libqgis-analysis2.14.8 but it is not going to be installed
        Depends: libqgis-app2.14.8 but it is not going to be installed
        Depends: libqgis-core2.14.8 but it is not going to be installed
        Depends: libqgis-gui2.14.8 but it is not going to be installed
        Depends: libqgis-networkanalysis2.14.8 but it is not going to be installed
        Depends: qgis-providers (= 1:2.14.8+20trusty-ubuntugis) but it is not going to be installed
        Depends: qgis-common (= 1:2.14.8+20trusty-ubuntugis) but it is not going to be installed
        Recommends: qgis-plugin-globe but it is not going to be installed
        Recommends: qgis-plugin-grass but it is not going to be installed
        Recommends: qgis-provider-grass but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

I have removed all the packages and the libraries so I do not have any library related to postgis, postgres and QGIS with the Synaptic Package Manager and all with the commands:

sudo apt-get purge qgis
sudo apt-get autoremove qgis
sudo apt-get --purge remove qgis
sudo apt-get update

But the problem is not solved yet.

I tried to remove the repo of ubuntugis and/or debian on the file /etc/apt/sources.list.

Nothing.

And I do not understand why, using this source http://qgis.org/ubuntugis/ I still have issues about the dipendences.

What should I do?

I am planning to format my laptop hoping this will help me with broken packages..
for the broken packages, I have already tried:

sudo apt-get update –fix-missing
sudo apt-get install -f

but nothing fixed.


texts: How to set the URL from an icon in OpenLayers 3?

texts:
I've tried to set the URL from an icon but without success. In fact, I've been searching in OpenLayers API but I haven't seen any method to set it. In OpenLayers 2 I've seen that It's possible.

I want to be able to change the URL at any time.

Do you have any idea?


texts: How to see if a function is already registered to an event on an OpenLayers map?

texts:
I am creating a data editing map using OpenLayers. I have two main states that the map editor can be in - editing existing features, or adding new ones. When a user clicks the "Edit Features" button, the editOnClick event is registered to the map. When the user clicks the "Add Feature" button, the drawFeatureOnClick event is registered to the map. I want to include a simple if statement to prevent OpenLayers from registering the same event multiple times, as well as remove the other button's event if it is active. Here's some pseudocode for the edit button:

$("#editfeatures").button({text: "Edit Features"}).click(function() {
    if (<editOnClick not already registered>) {
        map.events.register("click", map, editOnClick);
    }
    if (<drawFeatureOnClick currently registered>) {
        map.events.unregister("click", map, drawFeatureOnClick);
    }
}

I tried inspecting the map object using the Chrome debug tools, but I couldn't find an obvious way to tell if an event is currently registered on the map. How can I see what event listeners are registered on the map?


texts: Importing labels for check in/out from an SDE to file geodatabase

texts:
I am working on a large project where I am using the ArcMap 10 'distributed geodatabase' feature to check out a file geodatabase from our SDE file. I have tried several ways to avoid having to manually re-label the point features once I have checked out, but come up short every time.

I have saved the symobology to a layer file and can import that easily once checked out, but as far as the labels go, I cannot seem to easily import those. Is there any way to save the labeling with the layer, or is manually labeling the points through the properties menu the only way? Any suggestions are greatly appreciated and I will be happy to clarify if needed.


texts: Alternatives to Maplex Labeling Engine for ArcGIS Desktop?

texts:
I am running ArcGIS 10 (ArcView) with the standard labeling engine. I am very familiar with Maplex and have used it in the past but it still lacked the "feature-linked" functionality given the type of license I use. (Arc Info - Maplex combination was quite good)
Still, I found Maplex to be better than the built-in labeling engine but not stellar. Decent label masking for elevation contours and orienting labels up-slope, etc. are still difficult if not impossible to acheive even with maplex.
Area there any good alternatives?


texts: Using deactivate function of ArcGIS Python Add-in Tool?

texts:
I have not been able to figure out how the deactivate function of a Python Add-in Tool works in ArcGIS.

The help page that mentions the function, and that it 'causes the tool to no longer be the active tool' is here: http://resources.arcgis.com/en/help/main/10.1/index.html#/Tool/014p00000027000000/. I cannot find any examples of it being used.

The Tool class looks like this:

class ExampleTool(object):
    """Implementation for Example_addin.tool1 (Tool)"""
    def __init__(self):
        self.enabled = True
        self.shape = "NONE"
    def onMouseDownMap(self, x, y, button, shift):
        # My code here, then..
        self.deactivate()
    def deactivate(self):
        # Do I need to add something here..?

The help file seems to imply that when the function is called, the tool will deactivate / lose focus, but maybe I need to add something to the function to make that happen..?


texts: Conceptual advice on proximity analyses in arcpy

texts:
I have a question concerning proximity analyses in arcpy (ArcGIS 10.7). Forgive me if it takes a while to explain, but I've added a figure below to illustrate my problem.
Explanation
I'm trying to discover potential lakes (i.e. exposed subglacial overdeepenings) in the ice-free beds of ~1000 glaciers. I have the locations of all overdeepenings and the position of the glacier front at several points in time. So far, I have

created a line perpendicular to the flowline of the glacier
clipped the whole glacier shp with that line to get the ice-free bedrock
clipped the shp containing all the overdeepenings with that bedrock shp giving me all ice-free lakes

All of that works fine as long as the glacier runs in a more or less straight line. If it bends back on itself, I get problems as some parts of the upstream glacier are wrongly considered ice-free and lakes are incorrectly found (red circles in the figure below).
So far, different approaches to solve this have not worked:

I thought of trying to clip the glacier just along its main glacier bed to exclude all upper parts. However, sometimes the upper parts of the glacier converge at mountain ridges with the ends of lower tributary glacier arms.
This means, my polygon is at these points not easily separated into lower and upper parts because of the missing polygon separation (orange area in the figure). So, I think I have to use this 'endless' perpendicular line to clip the whole thing.
I have done Near analyses with different points along the upper flowline. I was hoping that some of them would be closer to my incorrect lakes than the point at the true glacier front. This way, I could have determined that the lakes were further upstream.
Unfortunately, this does not always work. In the figure, the point at the glacier front (A) is still closer to the wrong lakes than point B or C. And as the upper flowline consists of nearly 600 points, I want to avoid doing a near analysis for every point hoping that one will be closer.

Now the question
I'm looking for advice on how to solve this issue. Is there a way using arcpy or Python in general to exclude all sinks upstream from the glacier front?
Maybe a better way to separate the glacier in upper and lower parts?
Or another proximity analysis?


texts: Creating regular polygon grid in PostGIS

texts:
How to create, on the shape of a polygon, a regular grid of polygons/squares of an given size, in PostGIS?
I have thought about a function like Creating regular point grid inside polygon in PostGIS only for squares, so that squares can be 5m x 5m or even 10m x 10m. But have no idea to change this in the correct way.


texts: Generating raster with its spatial reference from array (or matrix) by using pyqgis?

texts:
I used arcgis to generate rasters from txt files by a certain python api,but i am not familiar with qgis,and now i need generate some rasters from txt files by using pyqgis.

The resulting rasters should be with their spatial reference so i can do some calculation across rasters with different spatial extents.


texts: Uncaught SyntaxError: Unexpected token : when using Google place API and Bloodhound

texts:
My code looks like this:

var placesBH = new Bloodhound({
  name: "PlaceNames",
  datumTokenizer: function(d) {
    return Bloodhound.tokenizers.whitespace(d.name);
  },
  queryTokenizer: Bloodhound.tokenizers.whitespace,
  remote: {
    url: "https://maps.googleapis.com/maps/api/place/textsearch/json?query=180 St Kilda Rd, Melbourne VIC 3006&key=***&" + "callback=?",
    filter: function(data) {
      console.log("data: " + data);
      return $.map(data.results, function(result) {
          return {
            name: results.formatted_address,
            lat: results.geometry.location.coordinates.lat,
            lng: results.geometry.location.coordinates.lng,
            source: "PlaceNames"
          };
        });
    },
    ajax: {
      dataType: 'jsonp',
      jsonpCallback: 'callback',
      complete: function(jqXHR, status) {
        console.dir("jqXHR: " + JSON.stringify(jqXHR));
        console.dir("Status: " + JSON.stringify(status));
        $('#searchicon').removeClass("fa-refresh fa-spin").addClass("fa-search");
      }
    }
  },
  limit: 10
});

The result returned by google seems correct but I got the following error:Uncaught SyntaxError: Unexpected token :

I tried different option: with and without dataType: 'jsonp', with and without jsonpCallback: 'callback', with and without  + "callback=?" at the end of the url, ...

From the console, I got the following:
jqXHR: {"readyState":4,"status":200,"statusText":"load"}
Status: "parsererror"

Still got the same error :(


texts: How to overcome invalid input geom and self-intersection when intersecting shapefiles in R

texts:
I downloaded two shapefiles from IUCN (available here and here) to determine the range overlap between two species. When I bring the shapefiles into R, I am able to plot them on a map together, however, when I try to calculate percent overlap, I get this error:

 p <- readOGR(getwd(), layer="species_7928")
 w <- readOGR(getwd(), layer="species_21314")
 EPFU<- gArea(SpatialPolygons(p@polygons))
 intersected.area.tabr <- gArea(gIntersection(SpatialPolygons(w@polygons),
 SpatialPolygons(p@polygons)))
 (intersected.area.tabr/EPFU)*100

  Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id,drop_lower_td,
  unaryUnion_if_byid_false,  : TopologyException: Input geom 0 is
  invalid: Ring Self-intersection at or near point -82.517899999999997
  8.2930355900000006 at -82.517899999999997 8.2930355900000006

I've tried the suggested workarounds with gBuffer (width=0) and gUnaryUnion, but these also return errors.
I was also able to take a separate map provided by a friend from a different project and each of these files individually worked when I tried to calculate the intersection.
I can open these files together in ArcMap. I've run the check geometry tool there and the output table was empty. Is this occurring because there are multiple elements (discontinuous portions of the map)? I do not necessarily have to do this in R, if there is a straightforward way to do the same in ArcMap, but this is my first experience with any sort of GIS data.


texts: How to add a title in a GetLegendGraphic Request

texts:
The GetLegendGraphic request works well. I display the legend at a fixed position in my viewer using OpenLayers. But I would like to know if there is any way to add some custom text to the legend.

The legend that geoserver produces is like this: http://docs.geoserver.org/latest/en/user/_images/samplelegend.png

But I need to produce a legend with a title: http://minus.com/m4AP102Jj

Does anyone know if the GetLegendGraphic request in GeoServer, or the SLD style have a way to add this?

UPDATE (23/08/12)

For now I´ll use a litle trick to add a title by adding in my SLD another custom rule above all others that not match any other rules:

Example:

<sld:Rule>
  <sld:Name>[NAME]</sld:Name>
    <sld:Title>[TITLE]</sld:Title>
    <ogc:Filter>
      <ogc:PropertyIsEqualTo>
        <ogc:PropertyName>[PROPERTY_NAME]</ogc:PropertyName>
        <ogc:Literal>[SOME_VALUE_THAT_NOT_MATCH_OTHER_RULES]</ogc:Literal>
      </ogc:PropertyIsEqualTo>
    </ogc:Filter>
    <sld:PolygonSymbolizer>
      <sld:Fill>
        <sld:CssParameter name="fill">[MATCH_BACKGROUND_COLOR]</sld:CssParameter>
      </sld:Fill>
    </sld:PolygonSymbolizer>
</sld:Rule>


texts: Is there a way to filter the colors out of an image?

texts:
I am using ArcGIS 10.2 to show a facility's size and position in a project area. I was not given any shape files or survey files of the area. I only have a deed (.pdf) with a B&W sketch of what the area of the project is. I created a JPEG from the deed image and have used the georeferencing tool to lay it into the correct position in my map as a raster. I would like to show only the black portions of the jpeg and filter out the white background. Is there a way to do this?

Here is what I have right now. I would like to get rid of the white, completely leaving only the text and black lines. I realize that I could just trace the image--creating a new feature--but for my purposes the preservation of the original writing and lines from the deed is necessary.


texts: leaflet layer control checkbox options min/max

texts:
I have a straight forward leaflet map that is on a webpage. I am using the angular-leaflet-directive and all works as expected except for the layer control box.

When only one baselayer is specified I get no control object (expected behavior).
When two or more baselayers are specified I get the control layer with toggleable options  (expected behavior).
However, I also get a min and a max checkbox that I did not configure, create or in any other way define.

I cannot find any reference to these in the directive, the leaflet library, any plugins or even find a mention of it when Googling.

Does anyone know what this is and where it is coming from??

The map otherwise works exactly as expected (with legends, markers, etc.).

Here is the layer object source:

angular.extend($scope, {

      pageClass: 'maps-data',

      texas: {
        lat: 31.555502,
        lng: -98.959761,
        zoom: 5
      },

      defaults: {
        // tileLayer: "http://{s}.tile.opencyclemap.org/cycle/{z}/{x}/{y}.png",
        // maxZoom: 14,
        path: {
          weight: 10,
          color: '#800000',
          opacity: 1
        }
      },

      legend: {
        position: 'bottomleft',
        colors: [ '#904', '#ff0000', '#28c9ff', '#0000ff', '#ecf386' ],
        labels: [ 'Texas', 'Drainage', 'Recharge', 'Artesian', 'Authority Zone' ]
      },

      markers: {
        eaa: {
          lat: 29.555502,
          lng: -98.959761,
          message: "Edwards Aquifer Authority",
          focus: false, // true means the message box will be visible initially.
          draggable: false
        }
      },

      layers: {

        baselayers: {
          mqosm: {
            name: 'Map Quest OSM',
            url: 'http://{s}.mqcdn.com/tiles/1.0.0/sat/{z}/{x}/{y}.jpg',
            type: 'xyz',
            layerOptions: {
              subdomains: [ 'otile1', 'otile2', 'otile3', 'otile4' ],
              attribution: '&copy; OpenStreetMap contributors',
              continuousWorld: true
            }
          },
          osm: {
            name: 'Open Street Map',
            url: 'http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
            type: 'xyz',
            layerOptions: {
              subdomains: [ 'a', 'b', 'c' ],
              attribution: '&copy; OpenStreetMap contributors',
              continuousWorld: true
            }
          },
          landscape: {
            name: 'Landscape',
            url: 'http://{s}.tile.thunderforest.com/landscape/{z}/{x}/{y}.png',
            type: 'xyz',
            layerOptions: {
              subdomains: [ 'a', 'b', 'c' ],
              attribution: '&copy; ThunderForest terms',
              continuousWorld: true
            }
          },
          bwmap: {
            name: 'B&W Map',
            url: 'http://{s}.www.toolserver.org/tiles/bw-mapnik/{z}/{x}/{y}.png',
            type: 'xyz',
            layerOptions: {
              subdomains: [ 'a', 'b', 'c' ],
              attribution: '&copy; OpenStreetMap contributors, CC-BY-SA',
              continuousWorld: true
            }
          }
        },

        overlays: {
          majorAquifers: {
            name:'Texas Major Aquifers',
            type:'geoJSON',
            url:'../../data/geojson/NEW_major_aquifers_dd_reduced100.geo.json',
            layerOptions: {
              style: {
                'color': '#000',
                'fillColor': '#09c',
                'weight': 1.0,
                'opacity': 0.7,
                'fillOpacity': 0.7
              }
            },
            pluginOptions: {
              cliptiles: true
            }
          }
        }
      }
    });


texts: Visualizing part of the attribute table in QGIS

texts:
I have a shapefile with a long attribute table made of different columns. Some columns have several levels (es. "Year": 2009, 2010, 2011,... or "Season": autumn, spring, summer, winter) and other columns have just numerical values related to the environmental parameters (es. "Temperature", "humidity", etc).

I would like to display only values of Temperature in autumn, I have tried different expressions in 'Categorized' or 'Rule-based', but without success. How can I make it?


texts: Updating identify tool every time feature is selected with PyQGIS3

texts:
I want to use the "Identify Feature Tool" in code (PyQGIS3). What I've found in the internet doesn't use the tool, people usually print in console point features and ids. I also found a tool that would open a window form with the point information but this method won't get attached to the docks or update itself.

I am providing examples of the methods I found below:

METHOD 1 (Printing features on console):

atriVector = self.iface.mapCanvas().currentLayer()
featget = atriVector.getFeatures()

atributes = [i.attributes() for i in featget]

print("{}".format(atributes[0][0]))

METHOD 2 (opening a feature form):

atriVector = self.iface.mapCanvas().currentLayer()
f = atriVector.selectedFeatures()[0]
self.iface.openFeatureForm(atriVector, f)

This method as I said before won't update itself, it would open a new window for every selected feature.

These two methods are the ones closer to what I am expecting (which is the actual tool from QGIS interface)

Here is an screenshot of it:

I want to use this tool with PyQGIS, but I have only found previous methods which are not close to what I expect.

See that in the image there is a selected feature and in the box the information of the feature is shown.


texts: OpenLayers 3 create bounding extent which includes two points

texts:
I'm working with openlayer3 in my ionic application and am stucked with the problem with creating a bounding extent to get a few points in the viewport.

My code looks like this:

var destLoc = [6.18581,50.88093];
var currentLoc = [-122.29925626,37.47836852];

var ext = ol.extent.boundingExtent([destLoc,currentLoc]);
map.getView().fit(ext,map.getSize());

But this code gives me a extent of the currentLoc and not the whole extent which includes both points.

What am I doing wrong?


texts: Snapping feature nodes doesn't work after feature creation (filter involved)

texts:
I've experienced a strange behaviour during editing mode:

I have set a filter to 2 features (polygons), part of on of my PostGIS layers: "name"='A' OR "name"='B'
After that I have created a new feature with the name='C'
After saving the layer edits, the new feature logically disappears due to the filter
So I have extended my filter to the new feature (all 3 features appear now in the canvas): "name"='A' OR "name"='B' "name"='C'

When I now want to edit the new feature it is not possible to snap to any node of the new feature.
'Identify Feature' and the 'Moving Feature' tools are working.

Strange, after moving the feature the node snapping works.

EDIT

When I define the filter "name"='A' OR "name"='B' "name"='C' a priori, before I create the new polygon with the name='C', it works.

Is this an expectable behaviour or a kind of bug?

I'm working with QGIS 2.8 on a Linux Mint virtual machine.


texts: Using Spatial Join, Dissolve and Aggregate in ArcGIS Pro ModelBuilder

texts:
I am creating a workflow using ModelBuilder in ArcGIS Pro 2.5.2 that will take updated parcel polygons and create/update the Zoning_Parcel (parcels symbolized at the parcel level) and Zoning_Districts (parcels that are dissolved by zoning and then aggregated into large districts that are used for smaller scale viewing).
It is important to note that every time new parcels are given to me I use the most current version of the Zoning_Districts to perform the spatial join at the beginning of the workflow.
Right now I have a model that can perform the spatial join and dissolve, but I am stuck at how to iterate through the ZONING field, creating batch selections of each of the 17 zoning types and then outputting them as separate features and also as a single merged feature.
I am confused on what Iterator to use and how to create the batched selections sets and aggregate them to crate the districts.
The snips reference the inputs and outputs of the model.

 NEW PARCELS (INPUT)
 EXSISTING ZONING_DISTRICTS (INPUT)
 NEW ZONING_PARCELS (OUTPUT)
 NEW ZONING_PARCELS_DISSOLVED (OUTPUT)
 NEW ZONING_DISTRICTS (OUTPUT)


texts:
I am fairly green with GeoTools, and having a hard time converting a raster (GridCoverage2D) into a vector format.

In the hours I've spent searching this site and Google, there seems to be multiple ways this could be done: PolygonExtractionProcess, RasterToVectorFactory, RasterToVectorProcess to name a few.  However, regardless of what JAR's I reference in my project, I can't seem to find any of these classes.

So, I am hoping someone can tell me the easiest way to perform raster to vector with GeoTools 10.


texts: arcpy update cursor: can't connect to a database table

texts:
I am using  updateCursor to retrieve some rows from a non spatial table using a connection of type "Spatial Database Connection" (That creates a .sde file) from ArcCatalog.

However, it doesn't work, I have the following error:

  (RuntimeError('ERROR 000972: Table must have Object IDs in order to be
  sorted.',)))

Even if my table is not spatial (there is no objectid field).

Also, when I try to connect to another non spatial table from another database, It works and I can retrieve the rows.

Does anyone have a solution?


texts: Masking areas outside of polygon using GDAL?

texts:
I am trying to "mask" raster cells that are outside of a polygon layer, setting them to "NoData". Also, I need to use GDAL, rather than R which I normally use, because the raster is too large to fit in memory.

So far, my workflow is -

## Burn -9999 into areas outside of polygons
gdal_rasterize -i -b 1 -b 2 -burn -9999 -l pol pol.shp modis_south_copy.tif

## Set -9999 to NoData
gdal_translate -of GTiff -a_nodata -9999 modis_south_copy.tif modis_south_copy1.tif

The first expression sets all cells that are outside of polygon pol.shp in raster modis_south.copy.tif to -9999. The second expression creates a new raster modis_south_copy1.tif where the value -9999 is replaced with NoData.

Here is a graphical illustration of the process:

I start with a raster named modis_south_copy.tif, here is an image of the first four bands -

After running the first expression with gdal_rasterize, areas outside pol.shp are replaced with -9999 -

Finally, after the second expression with gdal_translate, the copy modis_south_copy1.tif has -9999 replaced with to NoData -

The input files modis_south.tif and pol.shp are available here -
https://www.dropbox.com/sh/b9wtqgpylohuerx/AACyzhzfeGqgthGxMe1ZfW1pa?dl=0

Can I somehow specify that I want to -burn NoData instead of -burn -9999, in gdal_rasterize and, if so, what is the right syntax? (-burn NoData is just a guess and doesn't work)


texts: Reprojecting ECW raster in FME gives black border?

texts:
Using FME i want to reproject an ECW image to EPSG:3857 from EPSG:28354.

I intend to sit this image(and many more) over a large mosiac.

Using CSMapReprojector this re-projects with a boarder. I've read other discussions about how to clip and set the nodata value but none seem to work.

Can anyone please suggest the right method for this

Example below is something i've tried (also with the RasterNoDataSetter but to no success). I think it has something

I think i'm failing due to the varying RGB values in the black collar area as below


texts: For a project, I try to visualize an image of the world map in QGIs print layout. However, in QGIS, the round earth is flattened in a 2-D flat map. How can I include the curving of the earth, so that a 2-D map will be more representable? My layouts currently look like the arching is straightened, thereby stretching out countries.


texts: Polygons connecting to unknown point off QGIS map that was not designated?

texts:
I'm taking a methods course where we are currently learning how to take points in QGIS with an RTK GNSS system (Emlid Reach RS: https://emlid.com/reachrs/).

All of our points and lines show up OK, but when we attempt to make polygons, QGIS closes the polygon, but then attaches said polygon to some unknown point that we never desginated that goes way off of our map.

Does anyone have any suggestions on what may be causing this issue so we can fix it for next time?


texts: Get layer's "Feature Type Details" using geotools

texts:
I'm trying to do something like this:
Check Geometry Type of a Shapefile Using GeoTools?

... but getting from a layer (Shapefile source) in Geoserver.
In short, how can I get the data from this screen (Feature Type Details) using the Geotools API?

I have only this, but can't go further:

GeoServerRESTStoreManager storeManager;
storeManager = new GeoServerRESTStoreManager( new URL(RESTURL), RESTUSER, RESTPW);
RESTDataStore rds = reader.getDatastore(workspaceName, storeName);


texts: Large offset between two identical vector layers

texts:
I want to have a dxf layer to overlay another layer. The dxf layer is already saved to a shapefile with the needed crs. The layer has a large offset and nothing happens when I set the crs to the needed RD NEW/ EPSG 28992, it has a large offset. The feature is far away from the place it should be and it is rotated. I think this is due to the fact that there was no geometry assigned when it was drawn as a DWG/DXF.

How can I rotate the shapefile layer and move it to the right place?


texts:
I read everywhere that Leaflet is internally using a spherical web mercator (EPSG:3857), but whatever I try I get either Lat/Lon (EPSG:4326) or strange coordinates I dont understand:

map.getBounds().getNorthEast();          // LatLng(52.5262, 13.44907)
map.getPixelBounds().getTopRight();      // Point(2253845, 1375425) ???
map.getPixelOrigin();                    // Point(2252165, 1375425) ???

Let's try to project the spherical mercator:

L.Projection.SphericalMercator.project(
    map.getBounds().getNorthEast()
);                                       // Point(0.23473, 1.08117) !?!

All right, there is a known bug in Leaflet 0.7.3, so let's add the earth radius:

L.point(
    L.Projection.SphericalMercator.project(
        map.getBounds().getNorthEast()
    ).x * 6378137.0,
    L.Projection.SphericalMercator.project(
        map.getBounds().getNorthEast()
    ).y * 6378137.0
);                                       // Point(1497143.40914, 6895833.38915)

So, really, is it so hard to get the default projection coordinates of my map? Why do I have to call project() at all? What did I miss?

I am using leaflet 0.7.3 and did not change the CRS, using default EPSG:3857. Any ideas?

How to access Leaflet EPSG:3857 coordinates?


texts: adding geometries to a source, target table for multiple rows in pgrouting

texts:
I am having an issue trying to perform the pgr_dijkstra algorithm and getting back a line in my table for multiple rows.  I have a table that has a columns of source, target node id's but I cannot figure out how to feed them into the pgr_dijkstra function properly.  My table basically looks like this, but I need to fill in the routeline column with the shortest path line connecting the source and target.

ID | source | target | routeline
1  |  30    |  217   |  geometry(linestring)
2  |  217   |  2143  |
2  |  15    |  43    |

I would eventually like to make a linestring out of several of the rows from the roueline column, as I am basically trying to follow survey respondents around a network to all the stop locations they visited, but I can't seem to find a pgr function for multiple ordered nodes.  For now, i just figured I would connect the two nodes and then connect the lines together later.  Please let me know if you have a better approach!  I have been following the tutorial here : http://workshop.pgrouting.org/chapters/geoserver.html
My problem is when I run my code:

select a.uniqueid, a.source, a.target , st_makeline(route.geom) as routeline
from
     analysis.routing as a,
    (select geom from pgr_dijkstra('ways',a.source,b.target)
order by seq) as route;

I get the following error, and can't figure out how to pass the source and target from each row into the function.

ERROR:  invalid reference to FROM-clause entry for table "a"
LINE 1: ...uting as a, (select geom from pgr_dijkstra('ways',a.source,b...
                                                             ^
HINT:  There is an entry for table "a", but it cannot be referenced from this part of the query.

Does anyone know what is going on, I can't figure out how to reference my source from the reference table.  Eventually, I need to insert this into an update statement such as

update analysis.routing set
routeline = (select ***the querry above with results as linestring)


texts: Are there WMS servers providing maps of UTM zones, actually in the zone projection?

texts: There are lots of EPSG:4326 (lat lon) services out there but if one has to use an UTM projected map as base layer to display utm scenes, are there any WMS endpoints publicly available for that?


texts: Merging attributes into shapefile (US County-level data) using R

texts:
I'd like to make a plot of the US, with each county colored according to a scalar value (the average summer temperature).
I'm using R.
I have a shapefile with all the counties:
> shape = st_read(paste0(my_dir, "cb_2018_us_county_5m/cb_2018_us_county_5m.shp")
> head(shape)
Simple feature collection with 6 features and 9 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox:      xmin: -120.0726 ymin: 30.28043 xmax: -83.34348 ymax: 39.37874
geographic CRS: NAD83
  STATEFP COUNTYFP COUNTYNS     NAME                       geometry
       39      071 01074048 Highland MULTIPOLYGON (((-83.86976 3...
       06      003 01675840   Alpine MULTIPOLYGON (((-120.0725 3...
       12      033 00295737 Escambia MULTIPOLYGON (((-87.62999 3...

And a csv file with temperature data:
> temp_data = read_csv(paste0(this_dir, "temperature_data.csv"))
> head(temp_data)
  GEOID    State  County Temperature
1001     ALABAMA AUTAUGA          25
1003     ALABAMA BALDWIN          31
1005     ALABAMA BARBOUR          15

I feel like I should be able to match these up - but I don't see any identifying column that matches across both files. (The county names get close, but the US has a lot of Jefferson and Madison counties)
How can I merge the temperature data into the shapefile for plotting?


texts: Cast geography Polygon to MultiPolygon; function st_multi(geography) does not exist

texts:
I'm trying to cast a geography Polygon to a MultiPolygon resulting from a 'geographic' intersection with a circle on Earth of a 10'000m radius, but PostGIS is raising this error when using ST_Multi on that resulting Polygon:
ERROR:  function st_multi(geography) does not exist
LINE 2:  ST_Multi(
         ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
SQL state: 42883
Character: 15

The message is not 100% clear to me, because ST_Multi does exist as I successfully used it on many other queries on that same database.
So, I have the feeling that it is simply not able to work on a 'geography' polygon (or maybe the error lies somewhere else?).
Here's a minimal reproducible code snippet (postgis 3.1 / postgres 13.1):
SELECT CAST(
  ST_Multi(
    ST_Intersection(
      ST_Simplify(
        ST_MakePolygon(
          ST_Force2D(
            ST_GeomFromGeoJSON(
              '{"type":"LineString",
              "coordinates":[
                  [8,48,45],
                  [8.1,47.2,45],
                  [8.2,47.3,45],
                  [8,48,45]
                ],
              "crs":{"type":"name","properties":{"name":"EPSG:4326"}}}'
            )
          )
        ),
        0.001
      ),
      ST_buffer(
        ST_SetSRID(
          CAST(ST_MakePoint(8.1, 47.2) AS GEOGRAPHY), '4326'),
        10000
      )
    )
  ) AS GEOMETRY
);

How would you cast a Polygon to a MultiPolygon in that particular case?


texts: Setting up MapServer WMS and use it with OpenLayers

texts:
I have set up MapServer WMS on the localhost by the following map file:

MAP
  NAME "Prvi_WMS"
  SIZE 600 300              # 600 by 300 pixel image output
  PROJECTION
       "init=epsg:4326"
  END
    #EXTENT -180 -90 180 90
    WEB
    IMAGEPATH "/ms4w/tmp/ms_tmp/"
    IMAGEURL "/ms_tmp/"
    METADATA
      "wms_title" "My Global Map WMS Server"
      "wms_onlineresource" "http://localhost/cgi-bin/mapserv.exe?mode=map&map=c:/ms4w/Apache/htdocs/MapFile06_wms.map&"
      "wms_enable_request" "*"
      "wms_srs" "EPSG:4326"
    END
  END

  LAYER

    NAME countries
    TYPE POLYGON
    STATUS DEFAULT
    DATA "/ms4w/data/countries_simpl.shp"
    MINSCALE 1000
    MAXSCALE 1000000000
    PROJECTION
       "init=epsg:4326"
    END
    EXTENT -180 -90 180 90

    CLASSITEM 'NAME'
    CLASS
      NAME 'Bulgaria'
      EXPRESSION 'Bulgaria'
      OUTLINECOLOR 100 100 100
      COLOR 255 255 150
    END

    CLASS
      NAME 'All Countries'
      EXPRESSION ("[NAME]" ne 'Bulgaria')
      STYLE
        OUTLINECOLOR 100 100 100
      END
    END
    METADATA
      "wms_title" "World Boundaries"
      "wms_onlineresource" "http://localhost/cgi-bin/mapserv.exe?"
      "wms_srs" "EPSG:4326"
      "wms_enable_request" "*"
    END
  END
END

then, I checked GetCapabilities and GetMap requests and the results seem fine. However, I have problems getting the map inside the OpenLayers. Here is how I define my WMS layer:

    <title>Karta</title>
    <link rel="stylesheet" href="openlayers/theme/default/style.css" type="text/css">

Primer prekrivanja slojev in izbire podlag

The WMS layer shown in OpenLayers is completely white. What am I missing?


texts: Reading Coded Domains from ArcSDE 10 RDBMS using SQL/ADO.net

texts:
I need to read ArcGIS Domain Coded Value pairs (value/description) from SQL Server / Oracle in SQL (i.e. not using ArcObjects). This is an XML field. I can get the domain definition thus

select gi.definition from dbo.gdb_items gi, dbo.gdb_itemtypes git
where git.Name='Coded Value Domain' and gi.Name='MAJORTOWN'

(where MAJORTOWN is my domain name)...just looking for the next, XML-parsing part.

I want an IDictionary of the value,description.

Anyone done this?


texts:
Is there a way to use an ArcGIS Geocoding Server with QGIS?

I recently lost access to my ArcMap License, so an alternative would be very helpful.

The service in question would be:
https://tnmap.tn.gov/arcgis/rest/services/LOCATORS

My overall goal would be to use this like I would use it in Arcmap. i.e. setting an accuracy threshold, and feeding it a list of addresses and receiving the associated points.

I'm not partial to QGIS if there is another free alternative to make use of the geocoding service. I would prefer a GUI, but if there is a decent guide to a python script or something similar that could work too.

Edit: As mentioned above, I must be able to use the service linked above. The reason is that this service has the most accurate addresses for Tennessee, and has limits that can accommodate my needs (up to 250,000 addresses).


texts: Alternative to bring OpenStreetMap Basemaps into ArcMap

texts:
The standard ArcGIS OpenStreetMap basemap tiles have been performing quite slow as of late. I was curious if there was a (preferably free) alternative service to bring OpenStreetMap basemap tiles into ArcMap.

Most people here vastly prefer the OSM Maps to ESRI's version of Streetmap as a basemap. I was researching on OpenStreetMaps wiki and could not find one that was working in ArcMap sufficiently as of yet.


texts: ArcSDE database connection strings for ArcCatalog

texts: Anyone know where the Database Connection strings are located in Program Files for ArcCatalog?  I have all my Database Connection strings fixed and all my scripts use them.  My computer went down and I have a chance to grab files off the hard drive.  Instead of me re-creating all my database connections..I think I've seen a .txt file that has the connection parameters and connection string names.  I figure this would be under C:\Program Files\ArcGIS....some where here.... but I can't find them.  Anyone know?  Thanks.


texts: Query for finding intersection points between adjacent features within a buffer tolerance - PostGIS

texts:
I've developed the queries to find out intersection points (within 0.2 buffer) between line and polygon and also within adjacent polygons by using PostgreSQL 10 and PostGIS 2.4 database.
Line data is around 200,000 and Polygon is around 80,000 rows. The query took long time (aprx 120mins) to produce the result.

Line with Polygon:  The situation is:

Query is:

SELECT distinct ST_AsText((ST_Intersection(ST_Buffer(a.geom, 0.2, 'endcap=round join=round'), ST_ExteriorRing(b.geom)))),a.elemid,b.elemid FROM "Line" as a,"Polygon" as b WHERE st_Intersects(a.geom, b.geom)

Result is:

Polygon with Polygon:  The situation is:

Query is:

SELECT distinct ST_AsText((ST_Intersection(ST_Buffer(ST_ExteriorRing(a.geom), 0.2, 'endcap=round join=round'), ST_ExteriorRing(b.geom)))), a.elemid, b.elemid FROM "Polygon" as a,"Polygon" as b WHERE st_Intersects(a.geom, b.geom) and a.id != b.id;

below is the sample output with buffer tolerance

I have applied the buffer for the geometries from table 'a' and finding the intersection with the geometries from the table 'b'. Also there is no major difference from the above two queries except the a.id != b.id case.

How can I determine an optimal method to speed up the execution.


texts: How would one write a ArcGIS JSAPI layer to access WMS service via a Tile Server?

texts:
I want to request small WMS tiles instead of a single image for the entire map viewport because I know the request is going through a Tile Service (TileCache).  OpenLayers handles this with the WMS layer by setting singleTile=false.  JSAPI has the ability to set the maxHeight and maxWidth via the resourceInfo options attribute.

But the JSAPI does something unexpected but very practical when these values are set.  Instead of getting multiple images to fill the map it still gets a single image and stretches it to fill the map (see http://jsfiddle.net/LA2aH/ for an example).

What I appears to be missing is an option to get multiple images.

The reason I'm asking about this is I have a TileCache service which is accessible via a WMS service but the requests must have specific BBOX values.


texts:
I'm trying to modify the scale calculation in OpenLayers. Looking at the code, I couldn't understand what is going in this piece of code:

scale = (this.map.getGeodesicPixelSize().w || 0.000001) *

                inches["km"] *OpenLayers.DOTS_PER_INCH;

especially what is getGeodesicPixelSize().w and 0.000001 representing?

Bellow is the full block of code

updateScale: function() {

    var scale;

    if(this.geodesic === true) {

        var units = this.map.getUnits();

        if(!units) {

            return;

        }

        var inches = OpenLayers.INCHES_PER_UNIT;

        scale = (this.map.getGeodesicPixelSize().w || 0.000001) *

                inches["km"] *OpenLayers.DOTS_PER_INCH;inches["km"]

    } else {
        scale = this.map.getScale();
    }

    if (!scale) {
        return;
    }

    if (scale >= 9500 && scale <= 950000) {
        scale = Math.round(scale / 1000) + "K";
    } else if (scale >= 950000) {
        scale = Math.round(scale / 1000000) + "M";
    } else {
        scale = Math.round(scale);
    }

    this.element.innerHTML = OpenLayers.i18n("scale", {'scaleDenom':scale});
},


texts: Symbolize OSM Data Tool - ArcMap Error - The data necessary to complete this operation is not yet available

texts:
I am receiving an error using the ArcMap Symbolize OSM Data tool and  I cannot figure out the cause.
It seems to run flawlessly on the symbolize portion, then on Executing (Combine Layers): GPCombineLayers throws the error The data necessary to complete this operation is not yet available. (Exception from HRESULT: 0x8000000A). This is on data processed using the "Load OSM File" tool into Geodatabase.
Here is the  log below from the tool
Executing (Select Data (21)): SelectData c:\esrimain\arcgis\desktop107\desktop10.8\ArcToolbox\Toolboxes\Polygons.lyr "Geological (Polygons)"
Start Time: Tue May 11 07:57:59 2021
Succeeded at Tue May 11 07:57:59 2021 (Elapsed Time: 0.15 seconds)
Executing (Apply Symbology From Layer (20)): ApplySymbologyFromLayer "Geological (Polygons)" "c:\esrimain\arcgis\desktop107\desktop10.8\ArcToolbox\Toolboxes\Polygons.lyr\OSM Downloaded Polygons\Geological (Polygons)"
Start Time: Tue May 11 07:58:00 2021
Succeeded at Tue May 11 08:03:55 2021 (Elapsed Time: 5 minutes 54 seconds)
Executing (Combine Layers): GPCombineLayers 'Waterway (Polygons)';'Place (Polygons)';'Historic (Polygons)';'Power (Polygons)';'Aeroway (Polygons)';'Barrier (Polygons)';'Building (Polygons)';'Man made (Polygons)';'Shop (Polygons)';'Amenity (Polygons)';'Highway (Polygons)';'Railway (Polygons)';'Aerialway (Polygons)';'Leisure (Polygons)';'Military (Polygons)';'Tourism (Polygons)';'Natural (Polygons)';'Land use (Polygons)';'Boundary (Polygons)';'Geological (Polygons)' "OSM Downloaded Polygons"
Start Time: Tue May 11 08:03:55 2021
Succeeded at Tue May 11 08:03:55 2021 (Elapsed Time: 0.12 seconds)
Succeeded at Tue May 11 08:03:55 2021 (Elapsed Time: 13 minutes 41 seconds)
Executing (Combine Layers): GPCombineLayers 'OSM Downloaded Points';'OSM Download Lines';'OSM Downloaded Polygons' E:\OpenStreetMap\OSMSouth
Start Time: Tue May 11 08:03:55 2021
The data necessary to complete this operation is not yet available. (Exception from HRESULT: 0x8000000A)
Failed to execute (Combine Layers).
Failed at Tue May 11 08:03:55 2021 (Elapsed Time: 0.22 seconds)
Failed to execute (OSMSymbolizer).
Failed at Tue May 11 08:03:55 2021 (Elapsed Time: 17 minutes 31 seconds)
Failed to execute (OSMSymbolizer).
Failed at Tue May 11 08:03:55 2021 (Elapsed Time: 17 minutes 31 seconds)


texts: Should i use jp2 or cloud optimized geotiffs to efficiently access small windows from Rasters on S3

texts:
I am building a system where we have a lot of large rasters (Sentinel-2 bands) stored in S3. A lot of this data is stored in our own buckets so we can store whatever format we find most usable.

I need to efficiently and often access small windows (often less than 500x500 pixels) from those rasters. Since this happens a lot I need to be able to download those small windows fast and without tranfering more data than i need.

For this i am using vsis3 which solves this problem in a nice way, just as explained in the answer here.

I am using GDAL 2.2.4 from conda-forge which seems to use OPENJPEG 2.3.0.

The problem

I see 2 choices of how i can store my data.

As cloud optimized geotiffs
As JPEG 2000

The nice properties of option 2 is that the filesizes are a lot smaller, and that all data which we have not processed on our own can be bulled directly from amazons own S3 bucket

However the problem with option 2 is that it seems to download much more data to access the window than if i am using option 1.

I can see that the jp2 files for sentinel-2 images by default has tilesize 1024. So i created a cloud optimized geotiff with a tilesize of 256. This performs much better (in terms of how much data i have to download), so i expected the tilesize to be the reason. However I then tried to make a cloud optimized geotiff with a tilesize of 1024 and again it performs much better than the .jp2 file.

Here is a very crude visualization of the data transfer for each file

This show the data transfer required to fetch a 100x100 pixel window from a single band raster from each file type.

Now here is the question

Why do I have to download so much more data, when i try to access the same window from a jp2 file than from a geotiff file? It does not seem to be the tilesize, so what is the extra data i am downloading, and can i somehow avoid it?

I am just using the OpenJPEG driver, can this be the problem? and will a proprietary driver solve the problem i am describing?

Or do i simply have to bite the bullet and use cloud optimized geotiffs to access the windows faster, with the cost of some extra file size?


texts: ArcView: Trim polyline, that is shorter than specific dangle length (alternative to tool: Trim Line)?

texts:
I am looking for a way to trim dangles from a polyline shape, that are shorter than a specific length like the Trim Line Tool for ArcEditor and ArcInfo, but with ArcView-License.

I already discover, that the ET geowizards-Toolbox has a similar tool called "Clean Dangling Nodes". Unfortunately this Tool also create new features where other features are closer than the specified dangling tolerance. You can see this in the following picture:

Result (orange) of the tool Clean Dangling Nodes compared to the input data (blue)

Where the tool works like I want to is marked with green circles. But there are several situations where the tool creates new features (what I don't want) that are marked with red circles.

Any help is appreciated!


texts: Minimum distance from a point to nearest line among multiple lines shapefile

texts:
I am new to QGIS. I am looking for any plugin or any manual method in QGIS for finding minimum distance (perpendicular distance) from a point which is part of one shapefile to nearest line among many lines which is part of another shapefile.

MMQGIS is not useful for this purpose.


texts: Generating Point from Centroid of Polygon falls outside of it?

texts:
I have a feature class of polygons within a geodatabase and I'm trying to generate points based on the centroid of the polygon using the Feature to Point tool in ArcGIS 9.3.

However, I end up with several points that fall just outside of the polygon.

Is this the wrong tool to be using for this application?

Also,  am I making the wrong assumption about this tool generating points from the centroids of polygons?


texts: Wrong output wth group by and aggregate function in spatialite

texts:
I have a grid and for every cell I want to set the column "landuse" according to the largest intersecting polygon.

SELECT
  grid.PKUID,
  landuse.landuse,
  MAX(AREA(INTERSECTION(grid.Geometry,landuse.geom)))
FROM grid, landuse
WHERE Overlaps(grid.Geometry,landuse.geom) OR WITHIN(grid.Geometry,landuse.geom) OR  INTERSECTS(grid.Geometry,landuse.geom)
GROUP BY grid.PKUID

Without the max() the table would look like this

PKUID | landuse | AREA(INTERSECTION(grid.Geometry,landuse.geom))
168     mais      0.14
168     wheat     0.04
...

and this is what i get with the above query:

PKUID | landuse | MAX(AREA(INTERSECTION(grid.Geometry,landuse.geom)))
168     wheat     0.14
...

so the maximum value is selected but it doesn't tale the corresponding landuse. I PostgreSQL i would need to write the group part like this

GROUP BY grid.PKUID, landuse.landuse

but when I do this in my Spatialite query the aggregate function doesn't work and the result looks like i wouldn't have used it.

Any suggestions how to solve this?


texts:
I updated QGIS to version 3.6 using OSGEO4w

When I try to run any GRASS script, I get the following error message :

  This algorithm cannot be run :-(  The specified GRASS 7 folder
  "C:\OSGEO4~1\bin\bin" does not contain a valid set of GRASS 7 modules.
  Please, go to the Processing settings dialog, and check that the GRASS
  7 folder is correctly configured

It looks like the path doesn't exist, but I couldn't find a place to change it; it is not an option in the processing settings any more like in versions 2.x.

Is there a way to run those processing scripts in QGIS 3.6 and if so, how could I tweak the settings needed?

Edit : Here is what my processing options looks like


texts:
I have 61 pts and numerous rds. I would like to find the shortest distance to the nearest road using st_shortestline. I need 61 unique line segments, one for each point to its nearest road. Instead I am getting 61 line segments from each point to a particular road (not the nearest). How can I adjust my query to produce the proper results?

SELECT DISTINCT ON (p.id) p.id, p._unique, r.rd_name, ST_ShortestLine(p.geom,r.geom) into shortest_line FROM pois_sample p, allrds_sample r;

Outcome:

Desired Outcome:


texts: Searching value in string and replacing it with double quotes in ArcPy?

texts:
Adds a new field Match_PINCODE
Searches for a 6digit number in a string from Address field like "aaaaa,dsjfhdjfhd768543"
set "Match_PINCODE" field with the 6digit row by row.
Replaces 6digit number with ""

The 3 steps take almost less than a minute but calculatefield is taking a lot of time.Does anyone have better approach? so that i can fetch 6 digit number from one field update it into another field and finally replace 6 digit number with doublequotes.

import arcpy,os,re
import arcgisscripting

ints = []
cnt=0

gp=arcgisscripting.create()
gp.toolbox="management"
fc=arcpy.GetParameterAsText(0)
fc_fld=arcpy.GetParameterAsText(1)

try:
    arcpy.AddMessage("Adding Field ")

    gp.addfield(fc, "Match_PINCODE", "TEXT", "#","#",10, "#","NULLABLE","#", "#")

except:
  print "Error occured while adding Fields"
  print gp.GetMessages(2)

try:
    rows = arcpy.SearchCursor(fc,"","",fc_fld)
    for row in rows:
          cnt=cnt+1
          value=row.getValue(fc_fld)
          p=re.findall(r'\d{6}\b',value)

          for s in p:
            ints.append(int(s))
            finalval=int(int(s))
            if str(finalval) <> "":
                whereclause=fc_fld + " LIKE " + "'" +"%" + str(finalval) + "%"+"'"
                cursor = arcpy.UpdateCursor(fc,whereclause)
                row = cursor.next()
                while row:
                    row.setValue('Match_PINCODE', str(finalval))
                    cursor.updateRow(row)
                    row = cursor.next()

            arcpy.SelectLayerByAttribute_management(fc,"NEW_SELECTION",whereclause)

            arcpy.CalculateField_management(fc,fc_fld,'!'+fc_fld+'!.replace('+"'"+str(finalval)+"'"+',"")',"PYTHON_9.3")

            arcpy.SelectLayerByAttribute_management(fc,"CLEAR_SELECTION")
            arcpy.AddMessage("Row Count:" + str(cnt))

except:
  print "Error"


texts: Converting js time stamp to work in CMV application that read field values

texts:
This might be a somewhat basic question but I cannot seem to get the logic to work. I am having a search widget search through dates to find similar records of what an asset was last serviced in the field. I noticed that when I read the raw data output from the application, it is not in a usable format.

The input from my layer is in the formatt "mm/dd/yyyy hh:mm" but the raw output from my application in js is a string of numbers. I want to convert this string back into a real date so the user can search for a year, such as "2012" or whatever they want.

Here is the code block I am working with...

           attributeSearches: [
                {
                    name: 'Search For Hydrants',
                    searchFields: [
                        {
                            name: 'Hydrant Test',
                            label: 'Year Tested',
                            expression: ('TESTDATE LIKE \'[value]%\')',

                            placeholder: 'Enter a year',
                            required: true,
                            minChars: 4
                        },
                    ],

Anyone know the proper way to format the expression?


texts: How to edit table in locked ESRI StreetMap data?

texts:
Preface: I'm newer to GIS, but know the basics of ArcMap. I'm using ArcMap 10.1 on a 64bit.

I've been working with the new 2014 Release 3 of North America StreetMap Premium (TomTom) for ArcGIS, but am unable to edit the data table in the Streets layer (I'm able to edit Highways, etc.). I've tried to uncompress file gdb, but it's failed twice, possibly because of how large the dataset is.

I wanted to try to open up the Streets file in Microsoft Access and query the data that way, but it can't be opened in whatever format it is (I've tried creating shapefiles, as well as the original .gdb).

Anyone have suggestions or has worked with this set of data?


texts: Display feature length in legend in ArcMap

texts: Using ArcMap 10.2, how can I show the total feature length instead of feature count in the legend?


texts: How do I obtain layout coordinates of features using arcpy?

texts:
I just read the topic about Convert point XY to page units XY using arcpy? but I didn't find a solution to my problem.
I have a point layer file, I know coordinates in map view, I'd like to know coordinates in layout view. I suppose the main difference with the case of a dataframe is that the point is not a graphic element itself in layout view, so the referencing of the point is more tricky. Is there anyone who solved this (or similar) issue?
ps. the dataframe where I'm working is rotated, the solution proposed in the Convert point XY to page units XY using arcpy? post doesn't work for rotated dataframe. Does anybody knows what is the base point for rotation of dataframes in ArcGIS?


texts: Data defined position of table in print composer (QGIS 2.18)

texts: As the screenshot below shows: the data defined position of the text field is working, but that of the table doesn't. Can anyone confirm this behavior?


texts: PyQGIS - Adding layer from PostGis database on QGIS interface

texts:
I have the same problem, bellow my code if some one can help me please.
from qgis.core import *
from qgis.core import QgsProject
from PyQt5.QtCore import QFileInfo
from qgis.core import QgsVectorLayer, QgsDataSourceUri

def run_script(iface):
    uri = QgsDataSourceUri()
    uri.setConnection("localhost", "5432", "Base_test", "user", "****")
    uri.setDataSource("public", "mytable", "geom","")

    vlayer = QgsVectorLayer(uri.uri(), "mytable", "user")
    QgsMapLayerRegistry.instance().addMapLayers([layer])

When I run this code I haven't un error but i see nothing on my QGIS interface.
I'm using, QGIS 3.12.

from qgis.core import *
from qgis.core import QgsProject
from PyQt5.QtCore import QFileInfo
from qgis.core import QgsVectorLayer, QgsDataSourceUri
from qgis.utils import *
import processing

def run_script(iface):
    uri = QgsDataSourceUri()
    uri.setConnection("localhost", "5432", "Base_test", "Nassim", "*****")
    uri.setDataSource("public", "BPE", "geom")

    layer = QgsVectorLayer(uri.uri(), "BPE", "Nassim")
    QgsProject.instance().addMapLayer(layer)

This is my code whith the solution given by @Vincent Bré, but it still doesn't work


texts: Correlation between a raster and a point dataset

texts:
I have an interpolated surface (in raster form) with temperature values and I also have the measurements of temperature in certain locations that belong to the same region (with the raster). I want to plot in a scatterplot the relationship between the raster (estimated values) and the measurements (observed values) and compute the correlation between the two datasets and also the line that best describes the relationship among the two datasets.

Any idea how can I do this in ArcMap or QGIS?


texts: Extracting raster data to multiple polygons using ArcGIS Desktop?

texts:
I have a reclassified raster dataset (10 classes) and I want to find out the area of each class contained within multiple (hundreds-to-thousands) polygons.

I have all the polygons in the same shapefile and have used 'Extract by Mask' which has given me the raster for all the polygons, but now I want to do this for each individual polygon. I have been selecting each polygon individually and running 'Extract by Mask' which is working but will obviously take forever with 1000's of polygons.


texts: How can I edit gpx tracks, property timestamp and waypoints?

texts:
I want to reconstruct the historical movements of a military regiment. I am looking for a publisher, in the end I wish a gpx file. A priori there is full, yet I find none that corresponds: I need to display a map to determine where to put my points (with dates and hour-> timestamp) and sometimes some waypoints (for battle etc. ...).

Gpx files seems very appropriate, but what editor use ? I tried garmin basecamp, but I can't edit timestamp. I looked josm, but I do not feel that I can edit the tracks without losing the timestamp when I export ... well, I do not find (or I search very badly: /) an editor allows to create manually (visually) plots, edit timestamp and waypoints.

Can anyone guide me?

I want a editor to display a map on the screen and click at specific positions to create a point and give this point a timestamp and export to GPX file. In the end, I use the gpx file with osm and D3, but it's not my question. I not use a gps device, I don't want to import gpx file in a editor because I want to create track and waypoints myself from archive military documents.

I just want to transcript this kind of military log (dates, places and events) to a map. It's very long to do, so I prefer to have a GUI to input each event at once.

Exemple of a log (sorry it's in french).


texts: Repairing data source with joins to Excel without losing symbology?

texts:
I have a set of MXD files, that do not use relative paths, and related data that have been copied to a new computer.  Normally it's not a big problem to repair the data sources in the project, but in this case some of the layers are joined to an Excel file.  When I try to repair the data sources any layer that is joined to an Excel file looses it's symbolization.

Has anyone run into this problem and figured out a fix?


texts: Using raster calculator in QGIS gives solid black block, why?

texts:
I did an arithmetic average between rasters with raster calculator.
then to the overall average I should add 100.
But I have a problem: adding the result to the qgis I don't see anything.
But seeing the result in the output folder I have a black box.
The output is in GeoTiff format and I also get a tif.aux file
I tried to recalculate the statistical values of the data, For example, under Load Min/Max values, select Min/Max, then select Actual (slower), then press Load.
But with no result


texts: Create "mean" line from multiple lines using QGIS

texts:
I have multiple lines features (also have them in multipoints...) representing similar path (mountain trails, in this case) with various GPS precision, resulting in many lines close to each other, but not perfectly overlapping.

For the purpose of this project, I'm looking for a way to calculate a "mean" line from all of them and generate a resulting line of the most probable position of the trail.

How would you proceed, using QGIS or any other tool (I thought of OGR...)?


texts: PPP corrected base station coordinates relation to points

texts:
I have a set of Emlid Reach M2 L1/L2 receivers. They are attached to Harxon GPS 500 antennas and I am using a set of LoRA radios for communication between the two units. I live in Canada, so no free CORS stations are available and NTRIP networks are paid subscriptions.
I am using the PPP method through NRCAN to set my base on an unknown point and get a precise location from that occupation period (say 6-8 hours). This method works well to correct the base coordinates.
Now, if I had a lot of time, I could set my base for the observation period, wait 24-48hrs, process the data, re-occupy the same point and set the base station coordinates to manual which would effectively give me accurate rover positioning through RTK. But my workflow is to set the base, collect my points, do my other work and leave at the end of the day.
So I need to correct my points coordinates in post. Am I right thinking that I can do the following:

Get the PPP coordinate and calculate the offset in X, Y and Z from the base coordinates associated with my points
Apply this offset to my points

I am looking for a cleaner way to do this than by manually calculating the point positions which I fear can lead to mistakes. I have access to Global Mapper, QGIS, RTKLIB. Do not have access or funds for EZSurv and the like.


texts: How to create points around edge of viewshed results?

texts:
Is there a tool or module, or multiple tools, in ArcGIS that will allow you to automate the creation of points around the outer edges of a viewshed?

What I am trying to do is gather data about the elevations on the farthest parts of a viewshed. If I have points, then I can extract the elevation values to the attribute table and then do rise/run calculations from the original point to possibly create a horizon diagram in something like excel.

Is this possible?


texts: Python: Grabbing X,Y coordinate from matching fields

texts: Here is my table. What I am trying to do is write a python script to go over the Zip_Crid values. If there are matching values, I need to grab x1, y1, x2, y2 and put it in a new field. Haven't really worked with python before, where should I begin?


texts: How to merge dbf files in Python or arcpy?

texts:
I have over one hundred dbf files to merge. Each of files have two columns: ID and another field which differs from one file to another.  I want to merge those fields based on ID such that each raw corresponds to an unique ID containing all merged fields. For example:

file 1:
ID FirstName
1  Alan
2  Bell

file 2:
ID LastName
1  R
2  A

merged file:
ID FirstName LastName
1  Alan      R
2  Bell      A

Solution:

I have a series of files named as 'rawfile200101', 'rawfiles200102',...,'rawfile201104'.  Each of files have two columns: ID and another field same as the file name.  I want to merge those fields based on ID such that each raw corresponds to an unique ID containing all merged fields.

Different from Stata, the "merge" process is called "Join Field" in ArcGIS.  Since I have lots of these files, I will do this in batch.  Take the following as an example

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True

yrmn = ['200102', '200103', '200104']

for list in yrmn:
    inFeatures = "I:/rawfile200101.dbf"
    joinField = "VALUE"
    joinTable = "I:/rawfile{}.dbf".format(list)
    fieldList = ["rawfile{}".format(list)]

arcpy.JoinField_management (inFeatures, joinField, joinTable, joinField, fieldList)


texts: Looking for a deep learning platform for image recognition, preferably in QGIS

texts:
As part of a project I am working on, I am collecting an image training dataset consisting of various images of agricultural fields (details spared for confidentiality). The end goal of the project is to build a convolutional neural network (CNN) that will be able to classify high-resolution aerial images based on the training images provided.
I am new to deep learning and CNNs, but I was wondering if there is a recommended QGIS (or other GIS) package that will allow orthomosaics images to be classified in this manner? I was hoping to use an IDE such as Tensorflow or PyTorch; is there a way to use these IDEs in conjunction with perhaps a QGIS deep learning plugin? Any recommendations? I'm just trying to get started in this area.