texts: After turning census block data into a raster, how do I redistribute the density to reflect land cover types

texts:
I plan on rasterizing (terminology?) census block data relating to density to a regular 30 m by 30 m grid for use with NLCD data.  I then hope to distribute the density of the original census block among each result set of grid cells within the original polygon according to the land cover type dictated by the NLCD dataset.  For example, let's imagine a string of city blocks with a park, all found within one census block.  If within this census block the average density is 1000 people/km^2, we ultimately would want the resulting grid cell containing the park to have a lower population density than the grid cells containing the surrounding apartments.  This information may be captured by the NLCD dataset, with the park being classified as '21' (open but developed land) versus 24 for apartment buildings which could be classified as '24' (high density development).  Is there a method, or set of methods, within ArcMap that would allow me to:

Redistribute the density based on the NLCD value
Perform this redistribution only within the original census block polygon

I don't know if this is possible, but I figured I would ask.

Edit

Here's what I have been trying:

I converted the census block data to a raster, snapping it to the NLCD land cover data and ensuring the output grid cells were the same size.  Each cell contains the same info as the original census block (housing umbers, population, etc). I have been trying to figure out a way to edit the new grid so that if the NLCD grid cell contains a land cover type that doesn't support a business or housing (park, forest, farm, etc.) that it's housing and population values area set to 0.  My next move was going to be figuring out how to subtract the non-residential/non-business grid cells area from the area of the original census block, and force a recalculation of the density.  Converted the altered census block layer to grid cells again, with each having the new density value, and once again set the cells' housing, population, and density values that don't support a business or housing (park, forest, farm, etc.) to 0.

Kind of circular, and I'm not sure of all the steps (raster calculator for editing individual grid cells based on another raster's grid cell value?), but it seems to logically work out if I can get the parts to work.

On the other hand, @dof1985's answer seems to simplify the first part of my workflow (by directing me to the correct toolset), but I will still have to nullify the grid cells where parks, etc., do exits after the density has been recalculated for the census data.


texts:
pgr_dijkstra returns a dashed line. I tried with different databases and tables and always returns a dasher line.
Some explanation?


texts: How to get routes and download them from OpenLayers?

texts:
I am trying to generate a route in openlayers with a start and end points and some points during track.

I read about openroute: http://wiki.openstreetmap.org/wiki/OpenRouteService#RouteLink

But I don't know how can I include it on my webpage. My idea is about to have a button which when you press it, it generate the route between all markers in map.

I want that possibility of download it would be available.

In openroute webpage you can introduce points and it generates the route but I don't find an api to see how I can introduce it on my webpage.

Do you know other way to generate route in openlayers?


texts: Is there a better way to use the operator to select a group in an expression in Select by Attribute?

texts:
Is there a better method to do a group like in a row or in colum that you want to select from in using the Select by Attribute?  I have a point shapefile that I have has photo index. The field has all the numbers that looks like this 32-50,32-51,32-52, so on. This field Photo is in string/text data type.

I have a trouble using the expression but this one is the only way I can do

  "PHOTO" IN ( '11-94', '11-95', '11-96', '11-97', '11-98', '11-99',
  '11-100', '11-101', '11-102')

Here is the screenshot I have here


texts: How to connect existing postgis tables with Django models

texts:
I have some existing tables (shapefiles) and I want to connect each table with a django model. I imported the tables with shp2pgsql tool:

shp2pgsql -I -s 2100 -W 'iso-8859-7' world.shp public.world | psql -h localhost -d django_db -U django_user

Then I ran:

python manage.py inspectdb

and I copied the model of that table in my models.py :

class World(models.Model):
    aa = models.IntegerField()
    yparith = models.CharField(max_length=254)
    name = models.CharField(max_length=100)
    dep = models.CharField(max_length=50)
    arfakellou = models.CharField(max_length=50)
    sxolia = models.CharField(max_length=250)
    parakatath = models.CharField(max_length=250)
    yp = models.CharField(max_length=20)
    geom = models.MultiLineStringField(srid=2100)
    objects = models.GeoManager()

    class Meta:
        managed = True
        db_table = 'geoloc_world'

    def __str__(self):
        return self.name

    def __unicode__(self):
        return self.name

After that I ran:

python manage.py makemigrations
python manage.py migrate

Then when in django shell I type:

World.objects.all()

the output is an empty list.

I tried to use LayerMapping data import utility but it raises the following error:

 LayerMapError: Could not retrieve geometry from feature.

My load.py:

import os
from django.contrib.gis.utils import LayerMapping
from .models import World

geoloc_mapping = {
    'aa' : 'aa',
    'yparith' : 'YpArith',
    'name' : 'Name',
    'dep' : 'Dep',
    'arfakellou' : 'ArFakellou',
    'sxolia' : 'Sxolia',
    'parakatath' : 'Parakatath',
    'yp' : 'Yp',
    'geom' : 'MULTILINESTRING',
 }

geoloc_shp = os.path.abspath(
os.path.join(os.path.dirname(__file__), 'data',                                                                                'world.shp'),
)

def run(verbose=True):
     lm = LayerMapping(
     World, geoloc_shp,
     geoloc_mapping, transform=False,
     encoding='iso-8859-7',
     )
lm.save(strict=True, verbose=verbose)

Any suggestions ?


texts: How to Format QGIS Field names on Editor Form

texts:
I am using PostGIS table in QGIS with custom QGISFieldWidgetFactory. However, I wanted to rename the field names in QGIS Edit Form from say first_name to First Name using Python. Currently, I was able to format all labels/field names using the following method.

class WidgetWrapper(QgsEditorWidgetWrapper):
     ....
    def format_form(self):
        title = format_name(
            self.entity().short_name
        )

        title = QApplication.translate(
            'WidgetWrapper',
            '{} Editor'.format(title)
        )

        # Set title and format labels for QGIS form
        if self.parent.parent() is not None:
            self.parent.parent().setWindowTitle(title)
            label = self.parent.parent().findChildren(QLabel)[-1:]
            if len(label) == 1:
                text = label[0].text()
                formatted_text = format_name(text)
                label[0].setText(formatted_text)

      ....
   def initWidget(self, widget):
       ....
       self.format_form()

This is not a clean way and it doesn't work in QGIS 2.8. Also on QGIS 2.14, it can't work for the last field name. More over, it is not good performance wise as format_form() is called for every field in the layer which is not desirable.

Is there a fix for this?


texts: Errror using the Simplipy QGIS plugin: 'QgsRasterLayer' object has no attribute 'selectedFeatureCount'

texts:
I am adding a raster in QGIS 2.14.1 and get this new error:

AttributeError: 'QgsRasterLayer' object has no attribute 'selectedFeatureCount'

I am not programming, just using a QGIS menu.

I am using Linux Mint 17.3. update and upgrade made recently without effects.

This error started when I installed the Simplipy plugin. Now is disabled. I would like to test it.

What could I do?

This is the complete output:

An error has occurred while executing Python code:
AttributeError: 'QgsRasterLayer' object has no attribute 'selectedFeatureCount'
Traceback (most recent call last):
  File "/home/cubarro/.qgis2/python/plugins/simplipy/simplipy.py", line 274, in refresh_input_layer_list
    if selected_layer.selectedFeatureCount() > 0:
AttributeError: 'QgsRasterLayer' object has no attribute 'selectedFeatureCount'
Versión de Python: 2.7.6 (default, Jun 22 2015, 18:03:54) [GCC 4.8.2]
Versión de QGIS: 2.14.1-Essen Essen, exported
Ruta de Python:
/home/cubarro/.qgis2/python/plugins/processing_gpf
/usr/share/qgis/python/plugins/processing
/usr/share/qgis/python
/home/cubarro/.qgis2/python
/home/cubarro/.qgis2/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python2.7
/usr/lib/python2.7/plat-i386-linux-gnu
/usr/lib/python2.7/lib-tk
/usr/lib/python2.7/lib-old
/usr/lib/python2.7/lib-dynload
/usr/local/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages
/usr/lib/python2.7/dist-packages/PILcompat
/usr/lib/python2.7/dist-packages/gst-0.10
/usr/lib/python2.7/dist-packages/gtk-2.0
/usr/lib/pymodules/python2.7
/usr/lib/python2.7/dist-packages/ubuntu-sso-client
/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode
/home/cubarro/.qgis2//python
/usr/share/qgis/python/plugins/fTools/tools
/home/cubarro/.qgis2/python/plugins/DigitizingTools/tools
/home/cubarro/.qgis2/python/plugins/mmqgis/forms


texts: How can I get the profile graph on ArcMap to display specific points?

texts:
I currently have 2 layers:

1) DEM of existing terrain
2) Polygon in vector format of future water retention basin

I need to draw cross sections across this basin for eventual cut and fill calcs. I already know how to create a profile of the existing terrain based on my DEM. However I don't know how I can have the profile curve display the extents (or at least start and end points) of the basin on the graph. IS there any way to incorporate vector data into a profile drawn of a DEM? Or to have the profile curve display the specific points that are the  start and end points of the basin?


texts: Adding WMS and WFS layers from GEOSERVER in arcgis viewer for flex 3.0

texts: I am interested in adding feature layers from geoserver in arcgis viewer for flex config file and consequently display them on the web map.IS that possible? If yes, could someone share how they did that and perform querying on the features (points and polygons)?


texts: How can I get list of counties in a route from OpenStreetMap?

texts:
I have a route created by MapQuest Open Directions Service. I would like to get a list of all counties that route runs through. Does anyone know if this is possible?

I have tried figuring out if I can do it with Overpass or Nominatim, but  haven't found any way to do it yet.


texts: Disable Certain MapTips in Graphics Layer for ArcGIS Silverlight

texts:
I am using the ESRI Silverlight API to allow users to use the Identify functionality on a map.  The graphics layer on this map is set to use Map Tips, and a pin will be applied to the graphics layer when a user clicks on the map.  However, in specific cases, I want the map tips to be disabled.

I am using the MouseEnter event on the graphics layer to catch each graphic as the user rolls over them.  Once inside, I can determine whether or not I want the Map Tip to show up based on the graphics' attributes.

Is there a way for me to disable the Map Tips inside of the MouseEnter event?  I don't want to disable them for the entire graphics layer.


texts:
I setup ArcServer 10 with postinstall and restart and I am getting the following error when I try accessing ArcServer Manager through URL:

Forbidden

You don't have permission to access /arcgismanager on this server.

I can admire that i am a member of agsadmin group.


texts: Create raster shaded relief map in GDAL from ASCII files

texts:
I have just installed GDAL ver 1.11.4 and Python ver 2.7. I am trying to create a relief map using Ordnance Survey Terrain 50 ASCII data using gdal_merge.py.

Being new to GDAL, I don't know how to handle the many .asc files in the dataset. What are the tools available in GDAL to help transform these into a format like GeoTiff and merge the many separate files into one larger raster?


texts: Where can I deploy a JavaScript app for collaborative testing?

texts:
I'm developing a new ArcGIS Server JavaScript application, and I'd like to show the progress to the client, so he can provide feedback as I develop it.

The application consists of an Index.html file, a *.JS file containing my code, a *.CSS file for the styles and a *.JSON file containing the map configuration. I don't have access to a public-facing web server in my organisation.

Are there any sites (preferably free) where I can upload my files, so he can run and test the application?

I looked at github but from what I can gather it's a repository, rather than a place where my application will actually run (eg, test the Identify tool and verify that it works) - please let me know if that's not correct.

Thanks


texts: How to preserve original resolution of scanned 1-bit PDF raster with gdal_translate?

texts:
I am using gdal_translate to convert non-spatial PDFs to non-spatial TIFFs (intended for subsequent manual georeferencing).

This is my gdal_translate command:

gdal_translate -of GTiff -sds -co NBITS=1 --config GDAL_PDF_DPI 300 %input_pdf% %output_tif%

Here is what it looks like at ~8x magnification (note the "nibbling" around the edges of text and the complete loss of many "speckles"):

Input (rendered by Adobe Reader, disregard the slight blurring of edges):

Output (rendered by ArcMap):

If I do not specify the output bit depth (by omitting the NBITS option), I can tell that GDAL is re-rasterizing and resampling the input 1-bit image into a blurred grayscale image:

Ideally it would just copy each pixel in the input straight to the output.

How can I retain the original resolution without information loss in the output? Should I use something else, like ImageMagick, Poppler or GhostScript?

All this re-rasterization GDAL is doing is painfully slow, I might add, so I think I will try one of the alternatives for that reason alone.


texts: Export more than one object with catalog_apply

texts:
I am looking for advice on how to export more than one type of object when using catalog_apply() from the lidR package.
For example, in the code below the custom function fn_bufferTile() computes a lidar tile extent for a buffered lidar tile.  I would like to return or write both the buffered original las tile and the tile extent SpatialPolygons object. I am not sure what mechanism to use to pass back multiple objects, or how to simply write outputs directly from within a chunk?
fn_bufferTile = function(chunk)
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)

  ext_plyi = as(extent(chunk), "SpatialPolygons")

  # perhaps there is a way to generate custom paths here?
  #writeOGR(ext_plyi, "somePath1.shp")
  #writeLAS(las, "somePath2.las")

  # this works but write only the shapefiles
  return(ext_plyi)

  # doesn't work and fail because writing list is not supported
  #return(list(las, ext_plyi))
}

# set things up
ctg = readLAScatalog("folder/")
opt_output_files(ctg) <-  "folder/extent_{ORIGINALFILENAME}"
opt_chunk_buffer(ctg) <- 100
opt_chunk_size(ctg) <- 0
ctg@output_options$drivers$Spatial$extension <- ".gpkg"

#run function
out <- catalog_apply(ctg,  fn_bufferTile)

Not terribly important, but related, I am also wondering if the exported extent polygons can be written as separate layers to the same geopackage (e.g. the sf GeoPackage writing approach with append = T) instead of having to write to separate .gpkg files for every exported spatial object?


texts:
geoserver.log gives the following warning a bunch of times when loading my Image Mosaic layers

WARN [geotools.jdbc] - Could not find mapping for 'rast', ignoring the column and setting the feature type read only

I'm not using ImageMosaic JDBC, but I am using Image Mosaic with the time dimension attribute which generates a POSTGIS table pointing to the locations of the .tif files

fid, the_geom, location, ingestion
1,  010300...,  arnoldred_bloom_prism_2001.tif, 2000-12-31 17:00:00.000000

Do I need to add a rast mapping since the raster is not actually stored in the database. If not is there a way to get rid of the warnings?


texts: RuntimeError: ERROR 000732: Input Raster: Dataset ___ does not exist or is not supported

texts:
I am really new to ArcPy.

Here is my script:

# Name: Times_Ex_02.py
# Description: Multiplies the value of the input raster by a constant
#Requirements: Spatial Analyst Extension
# Import system modules and check extension
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

# Set environment settings
env.workspace = "C:\Users\srchang\Desktop\outputshp"

# Set local variables
inRaster1 = arcpy.Raster("tifdodge10m")

# Execute Times to multiply inRaster1 by constant
outTimes = Times(inRaster1, 50.5443234836703)

# Save the output
outTimes.save("C:\Users\srchang\Desktop\outputshp\outputtimes")

which returns the following error: ERROR 000732: Input Raster: Dataset tifdodge10m does not exist or is not supported.

I am thinking it is something to do with my path but i'm not sure


texts:
I have a ArcPy script to calculate Moran's I several times (50 to be precise). However, I am not interested in the Moran's I value, what I need is the z-score returned by the tool. Even though the Geoprocessing results window dispalays the z-score, I am not able to get the score through the script.

I have been through the documentation on ArcGIS resource center and the example given by them is:

arcpy.SpatialAutocorrelation_stats("olsResults.shp", "Residual","NO_REPORT", "GET_SPATIAL_WEIGHTS_FROM_FILE","EUCLIDEAN DISTANCE", "NONE", "#","euclidean6Neighs.swm")

The result returned by this script is the Moran's I index value.

Is there any way I can get the z-score besides looking at the Geoprocessing results window or the Reports?


texts:
I'm running QGIS 2.6 on OSGeo-Live 8.4 Virtual Machine based on Lubuntu 14.04.

I've set up a GRASS (GRASS GIS 7) environment in QGIS, and I have a merged .tif file of LiDAR images open in QGIS, so I want to use the grass tools to read the LiDAR image. Both the "grass tools" menu option and symbol are not highlighted and I can't select them (but some Grass tools are available such as add Grass raster/vector layer). What am I missing? How can I access the Grass tools?


texts: using css sprites in Openlayers 2 instead of icons

texts:
I'm using Openlayers 2 with Gmaps and a vector layer to display features with icons.

I'm using this code to set the icons on the features:

feature.style = {
            externalGraphic : iconUrl,
            graphicWidth : 30,
            graphicHeight : 20,
            fillOpacity : 1
        };

The problem is that I'm using the same icons in the site with css sprites and it would better to use them with openlayers rather than requesting the same icons from the server.

So I'm looking for a way to give each feautre a css class rather than a url for the icon.


texts: How to translate (reposition) a raster in Python?

texts:
I tried to use QGIS GeoReferencer to translate a raster layer, but somehow failed. It just goes to somewhere unexpected.

Now I decide to do it myself, with python. I think some libraries like GDAL must have done it. I only want to do an x-offset and y-offset with the tif, with no scale and rotation. Do you know what module I could use?


texts:
MODIS-NDVI product MOD13Q1 shows range from -1986 to 9976 but NDVI ranges from -1 to 1. What these value means? What do I need to do as corrections?

I have attached an image showing MODIS-NDVI ranges in ArcGIS.


texts: Covering a Polygon with Minimum number of circles with fixed locations

texts:
In QGIS, I have a vector polygon of a fixed size and shape. I also have a set of circles of a fixed radius and fixed location, which cover part of the polygon. What I would like to do is remove as many of the circles as possible while leaving as few gaps as possible, a gap being an area not covered by the circles. No circle should have overlap with another circle of more than 90% (or 80% or 70%- I am trying to run a few scenarios) of its area. I would also like to treat circles which go off the edge of the polygon in the same way as overlap, i.e. a circle with more than 90% of its area not covering the polygon would be eliminated. Inevitably there will be times when two or more circles could be chosen, I would like to prioritize them based on specific attributes. Here is a picture of the complete set of circles(pink) and the polygon I need them to cover (blue).

I am new to GIS and any extra elaboration on how to implement is greatly appreciated.


texts: Retrieve OBJECTID from query result from GISServer

texts:
Is there any way to obtain the OBJECTID value from returned result from GISServer query?

For example:

var searchString = "OBJECTID = 0001 OR OBJECTID = 0002";
var queryTask = new esri.tasks.QueryTask(url);
var queryParams = new esri.tasks.Query();
queryParams.outSpatialReference = new esri.SpatialReference({ wkid: xxxx });
queryParams.returnGeometry = true;
queryParams.where = searchString;
queryTask.execute(queryParams, resultsCallback, errorCallback);

If the returned result from above queryTask.execute() has only one feature, how do I know it is OBJECTID 0001 or OBJECTID 0002?

Of course..I can do the query one by one, but I am just wondering if there is better way to do it.


texts: QGIS / GRASS GIS How to make offset / parallel lines with Z values

texts: Please how to make something like v.parallel (GRASS function in QGIS) or offset (QGIS) from 3D line with keeping Z value (with z geometry). My vector 3D line is center of river and i need to make offset lines on both sides - result from functions above is perfect but it generates lines without Z geometry.


texts: ArcGIS Pro: Create polygons from outlines of unclosed lines

texts:
First a disclaimer: I have very little experience with ArcGIS, I usually use QGIS and only recently started using ArcGIS Pro.

I have a River dataset from the Natural Earth Site, then made its symbology thickness vary depending on its "strokeweig" attribute.

My goal is relatively simple: Creating a polygon that follows the outlines of this river line with variable thickness.

At first i saw two main ways of doing it:

1) Simply creating a polygon and drawing it (using the Freehand tool) so that it follows the outline of the rivers, however this is time consuming, and the fact the size of the river change when zooming prevents accurate drawing.

2) Exporting the river alone as a raster, then creating a polygon out of that raster, however I would have to reasign coordinates for this new polygon, which would bring some inacuracies.

-- Most of the Internet tutorials i saw recommended using the "Construct Polygon" tool, however it doesn't work on unclosed lines, and it only creates polygons out of the few rivers that join back.  (Exemple here: https://i.imgur.com/CrTByTD.png
https://i.imgur.com/ipt5hK8.png )

Is there a way to directly (as automatically as possible) create a polygon out of the shape/outlines of these river line (with their symbology) without either having to draw polygons by hand or exporting it as a Raster image.


texts:
I'm using QGIS to analyze some image of the NOLA area pre- and post-hurricanes. I'd like to calculate NDVI and a couple other indices for the images, but I want to exclude pixels with clouds in them from my final analysis. I already have an fmask band for each of the images, however I'm not sure how to incorporate that into the index calculation so that those pixels are excluded. Is this something I would need to do after calculating the indices? Ultimately I'm going to compute difference layers and use multidate classification for the change analysis.

I have access to QGIS, gdal, and otb, but not too much experience with the latter two.


texts: How to use the "Iterate Rasters" tool in ModelBuilder to link "Raster to Other Format (multiple)" tool to the "Extract by attributes" tool in ArcGIS?

texts:
I wish to link the "Raster to Other Format (multiple)" tool  to the "Extract by atttributes" tool. The "Raster to Other Format (multiple)" tool  only allows you to specify an output workspace, but the "Extract by atttributes" tool requires a specific raster to be named as input.

I am attempting to use the "Iterate Raster" tool to provide the specific raster needed as input to the Extract tool but cannot get the 2 tools to link.

Any suggestions?


texts:
I have a polygon layer and a polyline layer. I got the centroids of the polygon layer and used the Near Tool to find the distance from the closest point from the centroids to the polyline layer.

One thing I tried doing was to create a buffer from the points based on the Near value then finding the intersection but that did not work well.

One suggestion I got but was unsure what he meant and am unable to follow up with is:

  Use the Near function (not Generate Near Table nor Point Distance) with your original point layer that you added the coordinates to and specify another point layer as your near points. When the function is done, two variables will be added to the original point layer that contains the x and y coordinates of the point nearest to it, those will be your end x and y.

I think what confuses me is that I am not sure if I am supposed to find the nearest point using the polyline layer rather than the point layer.

Question:
How can I find the X/Y coordinates of nearest point on the polyline in order to create a line between that set of X/Y coordinates and a centroid point?

The end result is supposed to be:
The end result is that I want to make a line between the centroid and the nearest point on the polyline layer and connect it to the polyline layer. I tagged this Network Analyst because eventually I want to make a connected network.


texts: Modifying non-spatial data from an opengeo server

texts:
I'm developping an web-application using opengeo server. The data is stored into a PostGIS database, and some tables contains non-spatial data (that is related to the geographical data).

Is there a way to access the data in the non-spatial tables (to populate Ext.data.Store) and modify it?

I've tried many things from now, but none of them have succeed.

To read the data, I have for example set up a WFS layer, and then a WFS store in the web application. It gives me a read access, but no write access. Even if I get the record from the store, modify it, mark it as dirty and commit the changes, I see that the message that is sent to the server does not contain any modification (no updated records) - maybe because the store is not designed to work without layer.
The REST api does not provide access to the data.

An option could be to use web frameworks, but in that case, which ones are the best to use with Opengeo? I prefer to avoid this solution because it would require to manage an other server application (and it will complicate the management of the creditentials, since authentication is required to access the data).

Thanks!


texts:
I need to query all values in raster A which intersect with raster B.

Found this, but it's not accurate, and I'm getting ERROR:  column reference "rast" is ambiguous.

WITH first AS (
   SELECT id, (ST_Intersection(rast, rast, 0)).val
   FROM file_caolvkddzzlgmvyeofrh
   INNER JOIN file_xbwnbydpmwesokzmrcht ON ST_Intersects(rast, rast, 0)
)
SELECT id, COUNT(val), SUM(val), AVG(val), stddev(val), MIN(val), MAX(val)
FROM first
GROUP BY id;

This is probably very wrong - but basically, I need to figure out how to get all values of raster A that intersects with raster B.


texts: Rasterizing Shapefile, using other shapefile as Mask

texts:
I have a shapefile (ShapeY) of European extend and want to rasterize each time only a certain area of the shape (ShapeX) (i have the shapes indicating the prevailing extend) with python and QGis.

However so far i have not found a way to write the code in a manner that it already uses Shape X as mask.

My code reads as following:

err = gdal.RasterizeLayer(target_raster, (3,2,1), source_layer1, options=["ATTRIBUTE=%s" % RASTERIZE_COLOR_FIELD1])

The only solution i found was clipping the resulting layer afterwards with Shape X, but this takes significantly much time. I assume there is a smarter way of initially setting the mask while rasterizing?


texts:
I'm trying to convert a polygon to raster. This error pops up:

  ERROR 6: At least one of -3d, -burn or -a required.

What am I missing?


texts: Are Zonal Statistics possible in PostGIS2?

texts:
Now with postgis2 which supports rasters by default is it possible to do a zonal statistics analysis?

I've google it but I haven't found anything solid? Is there a tutorial outthere to get me started?

Can anyone give me an sql example on how to do it?

EDIT :

An updated (simplified) query according the blog Aragon mention :

CREATE TABLE sum_pop3 AS
 SELECT gid, SUM((ST_SummaryStats(ST_Clip(rast,1,geom))).sum)
 FROM perez_grid, ls_den
 WHERE ST_Intersects(geom,rast)
GROUP BY gid;


texts: Creating or updating geometry polygon column with new point using PostGIS

texts:
I have PostgreSQL tables with the following columns:

Table (1)

| Column        | Datatype                |
| ------------- |:-----------------------:|
| group         | TEXT                    |
| geom          | Geometry (Polygon)      |

Table (2)

| Column        | Datatype     |
| ------------- |:------------:|
| group         | TEXT         |
| lat           | NUMERIC      |
| long          | NUMERIC      |

I want to, for every row in table (2),  create or update the geometry columns in the first table (1) using the lat and long where the group column names are matching. Apologies if something like this has been answered before but I could not find a good example. I've tried queries like below (I know its wrong) but I can't seem to figure out how to use the PostGIS functions to get the behavior I want...

UPDATE TABLE_1
SET geom = st_setsrid(
    st_makepolygon(
        st_makeline(
            geom, st_makepoint(longitude, latitude)
        )
    ),
4269)
from (
     select
        long as longitude,
        lat as latitude,
        group
      from TABLE_2
) as subquery
where group = TABLE_1.group;

From discussion in comments:

Preferably, the geom column should be a polygon (if there are 4 or more points), if less than 4 it should just be a collection of points. Essentially I am trying to create my own geometries from some latitude and longitude data I have. Also I am using: WGS 84, EPSG 4326


texts: Geometry Type error preventing saving of layer after editing - geometry type is correct under layer information

texts:
I have a problem with a linestring layer that doesn't want to be saved after editing.
The layer itself shows 'LineString' under the information tab. When I add a line with the tool by snapping the line along some points, it doesn't let me save the layer. I get the error 'geometry type not compatible with the current layer'. I checked the points to see if they are 3D, but they are not according to the information tab.
I have uploaded a line that comes up in the 'check geometries' as invalid geo type to this link.
This one line also shows 'LineString' under the information tab.
Has anybody seen this before and has a solution? I thought I had found a workaround by just exporting the layer in editing mode, but it destroys, I think, the invalid lines in random ways. For example cuts them in weird places, and in random places.
I do have to note that not too long ago the information tab would tell me a wrong geometry type, not sure why.
I'm on a Windows 10 machine with the LTR Hannover QGIS.


texts: E: Unable to locate package geonode: Ubuntu 20.04 LTS

texts:
I am trying to follow quick installation guide for GeoNode on Ubuntu 20.04 LTS.
https://docs.geonode.org/en/2.8/tutorials/install_and_admin/quick_install.html
sudo add-apt-repository ppa:geonode/stable
sudo apt-get update; sudo apt-get upgrade; sudo apt-get autoremove
sudo apt-get install geonode
However, when I try to run last line, it shows me error:
E: Unable to locate package geonode
If you look at the screenshot: there are already errors while adding repository, such as "E: The repository 'http://ppa.launchpad.net/ubuntugis/ppa/ubuntu focal Release' does not have a Release file." and others...
What could be the solution?


texts:
In an attempt to join two rasters with a very small offset, the join on ST_ClosestPoint geometry column is logical to me:
SELECT s.geom, s.val, c.val
FROM s
JOIN c ON ST_ClosestPoint(s.geom , c.geom) = s.geom

But this does not return any results.
(PostGIS v2.2.1 and PostgreSQL v9.5)


texts: Determining WSEN of desired geographic area?

texts:
I'm looking to create a map of Afghanistan and I am in the middle of resampling binary geographic data into topojson format. The trouble is that my research of WSEN borders has returned nothing and I do not know how to adapt this section of the code to my area of interest.

For instance, in the makefile/d3.js tutorial I'm trying to learn from uses this syntax to display France:

# boxing:
crop.tif: ETOPO1_Ice_g_geotiff.tif
    gdal_translate -projwin -005.50 051.30 10.00 041.00 ETOPO1_Ice_g_geotiff.tif crop.tif
    # ulx uly lrx lry  // W N E S // -005.50 051.30 010.00 041.00

The author must assume the users have this as background knowledge, but sadly I do not and I could not find any more information on how WSEN borders work. To me, a rectangle whose vertices are coordinate pairs would make sense, but perhaps there is a reason for using the WNES border format?

How do I compute WSEN borders for the general case and what tools do I need for the task?


texts: Seeking tool for simplifying HUGE KML file?

texts:
I have a KML file which shows up fine on Google Earth but is missing lines on Google Maps. I have realized that this is because those lines have geometry that is too complex. The map in general probably has about 10x more points than it needs.

I have found a few sites online which will simplify the KML file but this file seems too big for that (8 Mb).

I read that Leaflet has a library called simplify which can help with this but I am new to JavaScript and would prefer a simpler alternative.


texts: Select features within an envelope in PyQGIS

texts:
I am trying to select (point) features within an envelope (buffer) of the start point of a line feature.

I have

line_layer = QgsProject.instance().mapLayersByName('DSC_Roads')[0]
chainage_layer= QgsProject.instance().mapLayersByName('DSC_Roads_10m_Chainage')[0]
for feature in line_layer.getFeatures(): # for each road select line and point features
    road = feature['RoadMntnc']
    if road == "Bamboo_Creek_Rd": #temporary for testing
        request = QgsFeatureRequest(QgsExpression('"RoadMntnc"= \'{0}\''.format(road)))
        l_sel = line_layer.getFeatures(request)
        c_sel = chainage_layer.getFeatures(request) #make sure points are on the same road
        for l in l_sel:
          geom = l.geometry().asPolyline()
          g = QgsGeometry.fromPointXY(geom[0])
          buffer = g.buffer(5,4) #10m buffer with 4 segments
          print (buffer) #test to see if the buffer works
          processing.run("qgis:selectbylocation",buffer,c_sel,['within'],0,0)
          chainage_selection=layer.selectedFeatures()
          print (chainage_selection) #test to see if selection worked

With this I am getting an error as there's probably something wrong with the selectbylocation that's based on https://gis.stackexchange.com/a/259441/2891

Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 12, in <module>
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\tools\general.py", line 106, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/OSGEO4~1/apps/qgis-dev/./python/plugins\processing\core\Processing.py", line 130, in runAlgorithm
    if context.feedback() is None:
AttributeError: 'int' object has no attribute 'feedback'

--- NEW CODE based on Joseph's answer

I tried the following but still get the same error

for l in l_sel:
    geom = l.geometry().asPolyline()
    g = QgsGeometry.fromPointXY(geom[0])
    buffer = g.buffer(5,4) #10m buffer with 4 segments
    buff_layer = QgsVectorLayer('Polygon', 'buffer' , 'memory')
    prov = buff_layer.dataProvider()
    feat = QgsFeature()
    feat.setGeometry(buffer)
    prov.addFeatures([feat])
    buff_layer.updateExtents()
    QgsProject.instance().addMapLayers([buff_layer])[0]
    processing.run("qgis:selectbylocation",buff_layer,c_sel,['within'],0,0)
    chainage_selection=layer.selectedFeatures()
    print (chainage_selection


texts: TauDEM fails to create stream network shapefile (-net option) when running StreamNet method on Ubuntu

texts:
When I try running TauDEM's "StreamNet" method including the optional -net flag to create the stream links / channel network as a shapefile, it fails with a lot of errors like:
ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_SetFieldDouble'.
ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_SetFieldDouble'.
ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_SetFieldDouble'.
ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_SetFieldDouble'.
ERROR 10: Pointer 'hFeat' is NULL in 'OGR_F_SetFieldInteger'.
More than 1000 errors or warnings have been reported. No more will be reported from now.

It looks like there have been some instances of related issues being reported and their fixes being introduced in that code base over time.
Thinking maybe the issue was it was looking for an existing dataset to write to, I tried using QGIS to create an empty shapefile hoping TauDEM might access it and write to it but that failed. (I wasn't surprised, but it was worth a shot.)
The rest of the method outputs are rendered, only the shapefile output fails.
Does anyone know if a different OGR file type will work properly, or if there is a trick to resolve this issue?


texts:
I need a printable administration map of a country. I would like to use it in Corel Draw for some colour manipulations and poster printing.
I have found the map on http://www.diva-gis.org/ but it is in shapefile format.

How can I convert it to any acceptable vector format like svg/cdr/eps/pdf?


texts:
I am having GPS and SFM (Structure for motion) trajectory in UTM coordinates. I wanted to convert the both trajectory to NED (North East Down).

Method 1

In First step: I converted the data from UTM ---> Lat, Long, Height

In second step: Lat, Long, Height ----> NED.

Method 2:

As another method I took the first data from SFM in UTM coordinate as reference. Then I subtracted the first data [easting, northing, height] from all other points. So, I got the data in meters.

The same data point [easting, northing, height] again subtracted from GPS data points. The GPS data is in now meters.

The plots looks as following:

Which method is correct for conversion of data to NED?

Because method 2 data is matching what the plot should look like in UTM!

Whereas method 1 might have precision problem due to two conversions, that's why the both data are close and seems overlapping each other. But that is not the case with original data.

What is the best method to convert UTM to NED frame?


texts: How do I search for a coordinate in QGIS?

texts: I have a very simple longitude and latitude I need to locate on the map and cannot find this for the life of me, sorry for the dumb question.  So basically I need to find 36.336429, -119.314682 on the map, I don't need a pin or anything special, I just need QGIS to point me to it.


texts: How to surround a clipped raster with transparency? (gdalwarp)

texts:
I am trying to perform the following task, but it does not seem to work:

I have got a raster file (raster.tif) containing an aerial photo, which I want to clip using a shapefile (clip.shp) that contains a polygon. The result I want is raster file (clipped.tif) containing the aerial photo inside the clip polygon and the rest being transparent, so I can overlay my clipped raster on an OSM map.

Basically I did this:

gdalwarp -cutline clip.shp raster.tif clipped.tif

The resulting clipped.tif shows the correct clipped aerial photo, but surrounded by a black background covering the extent of the original raster.tif file.

I have played around with the dstnodata and the dstalpha options, but with no luck. I am not exactly professor in raster processing, so the alpha channel stuff is new to me.

Any good ideas of how to perform such an operation?


texts: How to time function calls in OpenLayers?

texts:
We have an existing application that has a function that wraps the OpenLayers.Layer.Vector constructor, which gets its layer info from a KML source over http.  Once this function completes, I need to run another function which will read some of the layers data from the map object.

Normally in JavaScript I would of course use a callback when executing a remote AJAX call like this.  But in this case -- perhaps because I'm still pretty new to JavaScript -- I don't know how to delay the execution of my function until the Vector constructor is done.  There doesn't seem to be any ability to pass a callback to the OpenLayers.Layer.Vector constructor. Any suggestions?  Thanks.

[Addendum]

How do I add the listener to the layer before the KML is loaded?

var newLay = new OpenLayers.Layer.Vector(name, {
    displayInLayerSwitcher: true,
    visibility: visible,
    styleMap: MapLayers.styleMaps.polygonMap,
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
        url: url,
        format: new OpenLayers.Format.KML({
            maxDepth: 2,
            extractAttributes: true,
            extractStyles: useStyles
        })
    })
});

function loadEndListener() {
  alert('loading ended);
}

var events = new OpenLayers.Events(newLay);
events.on({"loadend": loadEndListener});


texts: Cast 3D geometry as 2D for ogr2ogr GeoJSON

texts:
I'm trying to take a shapefile that can be downloaded here (specifically the shapefile titled PVS_19_v1_vtd_21031.shp), convert it from EPSG:4269 to EPSG:3857, and output the GeoJSON into CSV with only XY coordinates. Currently, all steps work, except I'm still receiving a trailing 0 coordinate. Here's my command:

ogr2ogr -f csv -dialect sqlite -t_srs epsg:3857 -sql "select AsGeoJSON(geometry) AS geom, * from PVS_19_v1_vtd_21031" PVS_19_v1_vtd_21031.csv partnership_shapefiles_19v1_21031/PVS_19_v1_vtd_21031.shp -dim 2

and here's some sample output:

{"type":"Polygon","coordinates":[[[-83.66666399999999,36.71549199999999,0],

How can I drop the trailing 0 coordinate (I assume it's intended to be a Z coord?)


texts: How to programmatically create nested layer groups?

texts:
Using PyQGIS via the QGIS 2.4.0 Python console and the editor, I can create parallel groups for storing some layers with:

MyFirstGroup = iface.legendInterface().addGroup('Group No. 1')
MySecondGroup = iface.legendInterface().addGroup('Group No. 2')

The variables MyFirstGroup and MySecondGroup will store the GroupIDs (apparently 1 & 2).
But how can I create the groups in a nested way so that MySecondGroup is a sub-group of MyFirstGroup?

Looking at the API docs, I find the addGroup function twice - one time expecting the parent group with a QTreeWidgetItem and one time via a int parentIndex.

When using the latter form as

MySecondGroup = iface.legendInterface().addGroup('Group No. 2', parentIndex = 1)

then I get a TypeError: QgsLegendInterface.addGroup(QString, bool expand=True, QTreeWidgetItem parent=None): 'parentIndex' is not a valid keyword argument

When only using

MySecondGroup = iface.legendInterface().addGroup('Group No. 2', 1)

or

MySecondGroup = iface.legendInterface().addGroup('Group No. 2', MyFirstGroup)

then groups are not created nested but in parallel.

What is the correct way to achieve this or where do I get the corresponding QTreeWidgetItem object for MyFirstGroup from?


texts: ArcMap "Could not load data from the data source" Error

texts:
I used ArcMap 10.4.1 to connect to a postgre SQL database generated by osmosis on a server. I could see the feature rendering (e.g., points and lines) without a problem. However, when I tried to examine the attribute tables, some could be displayed while the others couldn't. For those tables that couldn't be loaded in ArcMap, I received the following error message:

However, in ArcCatalog, all tables can be viewed except this one:

Table Schema

Error Message in ArcCatalog

What could be the possible causes of this error?


texts: How to find chainage of points using LRS

texts:
I have a Shapefile of road that I have converted from a KML file (created on Google Earth) and a txt file containing lat long (from GPS) of points with some attributes. The points are not exactly on the line layer but are very close.

How can I get the chainage of all points from the first point using LRS?

There is an answer given earlier but I have not understood what are the lines route field, points route field and measure field in my layers. I tried to read several notes but to no avail.

Can someone kindly explain? I can send my layer files if needed.


texts:
if I'm programming a Python script, and in my script, I import ArcPy.

Is that downloaded with ArcMap?

Also, in that case, if I give the script to someone who doesn't have ArcMap, and thus don't have ArcPy, does that mean the script won't run as they're trying to import something that doesn't exist on their computer?

So if I'm sending a script to someone, would I have to send all of the libraries I'm importing with it?


texts: How to Use IWorkspaceFactoryLockControl in Java

texts:
I have tried to convert the interface but failed.

IWorkspaceFactory2 pWorkspaceFactory;
pWorkspaceFactory = new CadWorkspaceFactory();
ipWsFactoryLock = (IWorkspaceFactoryLockControl)pWorkspaceFactory;

It's failed when I use java.
I also tried the way to use the keyword NEW like this:

IWorkspaceFactory2 pWorkspaceFactory;
pWorkspaceFactory = new CadWorkspaceFactory();
ipWsFactoryLock = new IWorkspaceFactoryLockControlProxy(pWorkspaceFactory);

And I failed again.How can I use this method to release FeatureClass in Java


texts:
I have an sdc file with water bodies, i.e. oceans-rivers-lakes. I want to combine it with a shapefile with the boundaries of countries all over the world and produce two measures: 1) The average distance from waterways for each coutnry, and 2) the average fraction of land that is within e.g. 100km from water.

I work with ArcGIS 10.1 for Desktop.

Any suggestions on how to proceed?


texts: Handling large datasets with Erase/Dissolve in ArcGIS (Dissolve alternative?)

texts:
I have what should be a simple problem: erase a large multipart polygon layer from a larger polygon (see below for example). However, I believe the Arc Erase tool is trying to erase each part of the multipart from the larger polygon. I'm currently pushing 24 hours processing with ample storage and 12gb ram...

I have tested the Erase vs. Dissolve --> Erase workflow on a smaller dataset, and Erase runs substantially faster with the single dissolved polygon (~2 minutes for Erase with multipart, ~1 second with single feature in test dataset). Unfortunately, Dissolve also takes around 2 minutes with the same dataset.

So, is there an alternative to Dissolve that will return the outer extent of my feature class? I can not simplify the boundary with a convex hull or minimum bounding extent, I just need the feature mask. I'm open to python options as well.

UPDATE: Unfortunately, none of the provided answers gave me any great improvement in processing. I ended up buffering the inner dataset by 1 cm (smallest integer value available in Arc Buffer tool, and difference between inner and outer layers is ~150m) and erasing with that single polygon dataset.


texts: Why am I getting "sequence size must match size of the row" error?

texts:
I'm trying to use an Insert Cursor to insert a point geometry into an existing point feature class.

I'm finding the coordinate and appending them to a list called rowValues, ex:

rowValues.append("((({0}, {1}), 34))".format(X, Y))

This is giving me a list of value in this format:

((3656320.51246, 1613747.99484), 34)

Here is the cursor I am using:

cursor = arcpy.da.InsertCursor("Points",("SHAPE@XY","QuarterSections"))

When I run this code:

for row in cursor:
    cursor.insertRow(row)

I get this follow error:

Traceback (most recent call last):
  File "Quarter.py", line 110, in <module>
    cursor.insertRow(row)
TypeError: sequence size must match size of the row

Can I not use the string formatting? I've been trying to figure this out for awhile now, but I'm stumped.


texts: Choosing only one contour line to label in QGIS?

texts:
Is it possible to select only one contour line for labelling?

I have numerous labels for each line and I'll prefer to label only the main ones, such as 1000 m and 1500 m


texts: Switching from "layer selector" to buttons to toggle between layers in CartoDB

texts:
As of now, we have a map that you can toggle between layers using the "layer selector".

We want to instead have buttons on the map that, by clicking on one, you can switch between layers.  In our's, each layer has a different year's data with different lines and we want buttons for each of the three years.

We looked at this tutorial:
https://docs.cartodb.com/tutorials/toggle_map_view/
It appears to be out of date.  It directed us towards the "basic interactivity" tutorial for under "CartoDB.js is from the ground up".  We want to be able to style the map by the end user in the same way that the tutorial suggests.  The tutorial uses a template to edit that seems straight forward enough.  We want to edit our map similarly.  However, we can't find the javascript that we should edit.  If we copy the cartodb.js file under "publish", it gives some but not all of our customizations in an html file and nothing shows up if we try to bring it from a text editor into safari.

How can we create buttons and, more importantly, WHERE do we do that?


texts:
I have a legacy workflow that I would like to optimize and automate (at least partially).  The current workflow goes like this.  I receive a MS Access .mdb file containing a single table.  This table is then loaded into an an existing Access db that contained some macros ,which I will reference later, and renamed BTMUA.  To the table BTMUA several fields are manually inserted into the middle of the table.  Additionally two fields data types are changed from text to date to more accurately reflect the data they contain.    Once these changes have been completed a series of queries, via two macros, are run to parse, trim, and concatenate data from the original fields into the new fields.  Specifically Field A id parsed into fields A1 A2 using a left and right functions, likewise field B is parsed into field B1 and B2.  Any whitespaces are removed from fields A1, A2, B1, and B2 fields are removed.  Finally fields A1, A2, B1, B2, and C are concatenated in one of 7 different ways depending on which fields actually contain data.  Here ends the data modification and my use of Access.

I then use FME Desktop load this data into an existing table in my MS Sql based SDE GDB (after first truncating the table).

I would like to redesign this workflow and remove Access from the equation although I cannot change the fact that the data will originate in a .mdb file.  I would also like to automate the process as much as possible so it can run either with a few mouse clicks or on a fixed schedule.

Basically I am trying to figure out if it is best to use FME to load the data from the .mdb into my GDB and use either ArcCatalog and Model Builder or MS SQL Server management Studio an T-SQL to transform the data.  Or is it more efficient to utilize FME Desktop to do the data trans formations.


texts: Calculate Distance from Point outward via polylines?

texts: I am trying to figure out how best to calculate a 1.5 and 2.5-mile driving distance from a single point. Let's say we have a Fire Station connected to a road network and you want to find the 1.5 and 2.5-mile maximum driving distance along the road, is this possible?


texts: Sharing QGIS projects containing PostGIS layers

texts:
How do we share QGIS projects that includes PostGIS data?  If I create a project loading PostGIS data with reader permissions, nobody can modify the data.  If saved layers had been loaded from PostGIS with an administrator or editor's credentials then anyone opening that project can modify the data.
We could create projects for individual users with PostGIS data loaded with that individual's credentials which would give the required permissions and protection but if we share that project with co-workers then they are effectively accessing the data via the first individual's credentials.
I'm sure others must have come up against this but I couldn't find a solution.
We are using QGIS 2.18, PostgreSQL 9.0.23 and PostGIS 1.5.


texts: How to create a grid of 10degx10deg polygons for a shapefile using GDAL

texts:
I want to create a shapefile using the GDAL library. And I want this Shapefile, when created, to have a grid of square polygons. And where each polygon cell is 10 degrees by 10 degrees, ultimately I want this grid to cover the entire world. As well, I also want the shapefile to contain a decimal column 'LAND', which should be populated with 0.0 for all rows/cells.

Having confusion here because I've never used GDAL before and just sifting through the documents is confusing me. I think I understand how to create a shapefile something like the following:

from osgeo import ogr
from osgeo import osr
import gdal

# specifying that I want to work with a shapefile
DriverName = 'My Shapefile'
driver = ogr.GetDriverByName(DriverName)

Filename = 'points.shp'

Just a little lost here. Not sure exactly how this would be done through GDAL.
I know if I am creating the world, the extent of the world is -180 deg to 180 deg longitude and -90 deg to 90 deg for latitude. And that I should create a nested loop somewhere to loop the outer of longitudes and the inner of latitudes, stepping a value of 10 each time for spatial resolution for the grid.

It's just the whole-set up and the library functionally is confusing not sure what way to go about all this.

The general structure of the code I'd think would be something like:
- Import GDAL and other modules
- create a spatial reference
- Create an empty shapefile
- create a layer in the shapefile and assign it the spatial reference
- create a set of point geometries as a ring
- add those point geometries in an order to a polygon geometry
- create a feature and add the polygon Geometry to it
- put the feature in the layer

I examine the GDAL library website and found this:
link: https://pcjericks.github.io/py-gdalogr-cookbook/

and then more specifically a section that said "Create a New Shapefile and Add Data"
link: https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data

I looked at the following but they did not help so much:
link: Creating square grid polygon shapefile with Python?

and

link: Creating grid polygons from coordinates using R or Python

I am trying to piece all this information together in order to solve my task but it is proving difficult.

Please help how would I create a shapefile being a grid of 10x10degree polygons?


texts:
I cannot find any popular libraries for parsing (reading into memory) of LAS or LAZ  (lidar data) files in NodeJS (desktop Javascript)

https://github.com/yuyou/liblas4node/blob/master/lib/liblas.js

This depends on liblas which is deprectaed/replaced by PDAL, for which, no js bindings exist. Are there any other libs out there? npm search is a little clusmy, since 'las' is found in many completely non related packages


texts: Plotting XY on ArcMAP from csv file using arcpy

texts:
I have a .csv file with XY coordinates that can be converted to a .dbf file and loaded in ArcMAP. How do I plot the coordinates in that file in ArcMAP using arcpy.

Without arcpy, we could right click the loaded .dbf from table of contents and click "Display XY Data".


texts: Creating Direction data in QGIS between two layers?

texts:
Can I do it in Qgis, that there are two layers, a target layer and an input layer, and I ascertain with an algorithm so that the closest points on the input layer are on the target layer in what direction they are?

Measured in degrees, with designation east, west, north, south, north east, south east, northwest, southwest?

Creating a table that contains this data?


texts: Calculate new Lat/Lon from initial Lat/Lon plus Cartesian X,Y

texts: For small distances (below 10 km), assuming Earth is flat and I have initial Lat0/Lon0 and I have X(km), Y(km) distance From Lat0/Lon0, how can I calculate new Lat/Lon for the added X,Y?


texts:
I am working in QGIS 2.18.18. I have converted maps into PDF in a couple of prior projects, but the current map I am working on will not complete all the layers in the finished PDF.  Specifically, the OSM layer.  The OSM layer is visible across the bottom of the PDF map and barely makes it half way up.

What could be causing this??


texts: Changing the color and near ocean with QGIS

texts: The border of region near ocean looks looks very fat. How can I change color of border lines near ocean to make it thicker? I want to change color and thikness of line only near oceans


texts: map.remove throws type error after removing teh layer in ArcGIS JS API 4.3

texts:
I am trying to remove a layer from the mapView.

var LayertoRemove=map.findLayerById(layerID);
map.remove(LayertoRemove);

The layer is found and is actually removed, but after it's removed it continues to throw TypeError: this.controller.destroy is not a function, and the rest of the script is not run. when I remove these two lines of code, the script works fine.

the error in the console is pointing to MapView.js of the API and not to any line in my script. Any suggestions on what is going wrong here?

Update1: The FeatureLayers that I am trying to remove are not hosted layers, but generated from client-side geometries (graphics are created from two csv files taken from the client);  I suspect that the error is happening because the setting of one of the properties of the layer -- which I need to set when the layer is created.

Update2: the full error log is:

 Uncaught TypeError: this.controller.destroy is not a function   MapView.js:510
    at Object.d.detach (MapView.js:510)
    at Object.update (MapView.js:25)
    at k._animationFrame (init.js:329)

Update3: I am now almost sure that the issue is related to my featurelayer properties, as when I try map.remove with a hosted featurelayer, it works fine. The featurelayers are generated as follows:

function createLayer(graphics,renderer,type,fields,title){
          lyr = new FeatureLayer({
            source: graphics, //an array of esri/Graphic
            fields:fields,
            objectIdField: "ObjectID",
            renderer:renderer,
            spatialReference: {wkid: 4326},
            geometryType: type,
            title:title,
            id:title
          });
          map.add(lyr);
    }


texts: Customizing Europe Lambert Conformal Conic to match the Mediterranean Basin extent

texts:
This is a follow-up question to Which projection when dealing with least-cost path estimation across the Mediterranean area?. The Europe Lambert Conformal Conic has been suggested by @Fezter.

Shall I customize the projection modifying the two standard parallels in order to best match my area of interest?

In fact, I am not interested in the whole Europe, rather in its southern-most part, i.e. the Mediterranean basin (let's say from Gibraltar to Cyprus [west-east], from Corsica to Lybia [north-south]).


texts: Add text to shapefile in specific locations

texts: I have downloaded a grid with the latitudes and longitudes of the world but the numbers (e.g 0°, 30°, 60°, 90°) are not shown in the shapefile. Is there a way to add this text manually and then merge it with the shapefile layer?


texts: Why does Buffering in QGIS fail with this error: Argument 1 has unexpected type 'NoneType'?

texts:
I get this error on both 2.10 & 2.12 on two different machines running Ubuntu 14.04 with Python 2.7.6.

When trying to run the fixed-distance buffer on a shapefile with polylines I get the following error:

QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType' See log for more details

With full trace-back:

Uncaught error while executing algorithm
Traceback (most recent call last):
Traceback (most recent call last):|
File "/usr/share/qgis/python/plugins/processing/core/GeoAlgorithm.py", line 230, in execute|
self.processAlgorithm(progress)|
File "/usr/share/qgis/python/plugins/processing/algs/qgis/FixedDistanceBuffer.py", line 72, in processAlgorithm|
dissolve, segments)|
File "/usr/share/qgis/python/plugins/processing/algs/qgis/Buffer.py", line 81, in buffering|    outFeat.setGeometry(outGeom)|TypeError: QgsFeature.setGeometry(QgsGeometry): argument 1 has unexpected type 'NoneType'|

The first time I opened the "Processing" log to see more about this error, I noticed a problem with the SAGA installation but fixing that doesn't seem to solve this problem.


texts: localization by trilateration ( weighted nonlinear least squares )

texts:
I am writing an app for beacon localization and I am looking for a formula which,
given beacons of known positions and distance from target, finds the target position.

eg.

given beacon1 located at (x1,y1), distance from target=d1,
beacon2 located at (x2,y2), distance from target=d2....
find target postion

I have searched a question which seems to be useful:
Trilateration algorithm for n amount of points

but for the formula in the above website:

fit = NonlinearModelFit[data, Norm[{x, y} - {x0, y0}], {x0, y0}, {x, y}, Weights -> 1/observations^2]

I don't have an idea to implement this in the android app.
Also, it said it is using weighted nonlinear least squares fitting. I have searched for it and it is full of matrix and complex steps so I feel lost to implement it.


texts: QGIS "join attributes by location" fails

texts:
I am trying to join attributes from a vector polygon layer to a vector points layer wherever a point falls inside one of the polygons.

I tried to use the QGIS 'join attributes by location" tool, setting the points shapefile as the target layer, polygons as the join layer, and 'within' as the geometric predicate. The operation fails immediately upon starting. The log contents:

Algorithm Join attributes by location starting...

unknown See log for more details

The points shapefile is fairly small (~900 points) while the polygon layer contains many features (~3M), although I don't think this should be an important factor. Interestingly, I have been able to run the operation in reverse (switching the target and join layers and using 'contains' as the geometric predicate) successfully, however this is not what I need to do.

I have also tried using "Add polygon attributes to points" from the SAGA toolbox without any luck. (I will post the log contents if it helps)

Is there a possibility that there is something wrong with the geometry of this shapefile? I have also tried to save to a new vector layer, or to copy the points to a fresh layer, with no change in result.

Running QGIS 2.18 on Mac OS X El Capitan


texts: Reading raster cell values on specific PostgreSQL points with rasterio and updating PostgreSQL table with those values?

texts:
I have 9 UK wide raster files (TIFF) and about 3 millions points on a PostgreSQL table. I need to add, for each point on the pg table, the values of the 9 raster cells where the point falls in.

I am trying with the following piece of code:

pg_table = '.....'

#Here I create 9 new fields in the pg table to store the 9 values
query = """ALTER TABLE """ + pg_table + """
           ADD COLUMN U_R_20 numeric,
           ADD COLUMN U_R_75 numeric,
           ADD COLUMN U_R_100 numeric,
           ADD COLUMN U_R_200 numeric,
           ADD COLUMN U_R_1000 numeric,
           ADD COLUMN U_C_75 numeric,
           ADD COLUMN U_C_100 numeric,
           ADD COLUMN U_C_200 numeric,
           ADD COLUMN U_C_1000 numeric;"""
cursor.execute(query)
con.commit()

raster1_path = '....tif'
raster2_path = '....tif'
...
raster9_path = '....tif'

river_1 = rasterio.open(raster1_path)
...
river_9 = rasterio.open(raster9_path)

#Create list with all paths of tiff files
rasters = [river_1, ... , river_9]

#Select all midpoints and find out, one by one, value of raster cell where midpoint falls in (for all 9 rasters)
selection = etl.fromdb(con, "SELECT feature_id, ST_AsText(geom) FROM {0}".format(pg_midpts))
for (id, geom) in selection[1:len(selection)]:    #First element is only the attribute name
    #Convert the 'geometry' string of characters into an actual shapely point (type shapely).
    pnt = shapely.wkt.loads(geom)
    #Create list to store cell values where midpoint falls in
    values = []
    #Iterate through all 9 rasters to extract their values on current point
    for i in range(0, len(rasters)):
        #xxx.sample() takes an iterator over x, y tuples, so:
        for val in rasters[i].sample([(pnt.x, pnt.y)]):
            values.append(val[0])    #function sample returns a numpy array with only one element. We extract it with [0]

    #Insert values from rasters into table
    query23 = """UPDATE """ + pg_midpts + """
                 SET U_R_20={0}, U_R_75={1}, U_R_100={2}, U_R_200={3}, U_R_1000={4},
                     U_C_75={5}, U_C_100={6}, U_C_200={7}, U_C_1000={8}
                 WHERE feature_id = {9};""".format(values[0], values[1], values[2],  \
                                                   values[3], values[4], values[5],  \
                                                   values[6], values[7], values[8], id)
    cursor.execute(query23)
    con.commit()

This code seems to be working fine. However, it takes about 48 hours.

Is there any way to optimise what I am trying to do? Or maybe a completely different approach?

PS: By the way. ArcMap has a tool to do that (extract multiple values), but it doesn't work at all, it just messes up the whole table.


texts: Make Route with Python and network analyst in ArcGIS 10

texts:
Good Day,

I'm having trouble debugging my python code to create a network analyst layer in ArcGIS 10.0. Basically when run from a script tool the code completes successfully without actually doing anything or returning any sort of useful information, and when run from the embedded python window in ArcGIS 10.0 i get the following error.

"Runtime error : Failed to execute. Parameters are not valid. ERROR 000192: Invalid value for Time Attribute Failed to execute (MakeVehicleRoutingProblemLayer)."

The script is supposed to solve the route and create a layer file for the output
Here is my code below:

#Import system modules
import arcpy
from arcpy import env

try:
    #Check out the Network Analyst extension license
    arcpy.CheckOutExtension("Network")

    #Set environment settings
    env.workspace = r'C:\Layout3.gdb'
    env.overwriteOutput = True

    #Set local variables
    inNetworkDataset = r'C:\Layout3.gdb\NtsFDS\NtsFDS_ND'
    outNALayerName = "TestResult"

    impedanceAttribute = "CABLE_COST"
    distanceAttribute = "TRAVEL_TIME"
    timeUntis = "Minutes"
    distanceUntis = "Meters"
    inOrders = r'\Layout3.gdb\NtsFDS\Orders'
    inDepots = r'C:\Layout3.gdb\NtsFDS\Depots'
    inRoutes = r'C:\Layout3.gdb\NtsFDS\Route'
    outLayerFile = env.workspace + outNALayerName + ".lyr"

    outNALayer = arcpy.na.MakeVehicleRoutingProblemLayer(inNetworkDataset, outNALayerName,
                                                         impedanceAttribute,
                                                         distanceAttribute)

    outNALayer = outNALayer.getOutput(0)

    subLayerNames = arcpy.na.GetNAClassNames(outNALayer)

    ordersLayerName = subLayerNames["Orders"]
    depotsLayerName = subLayerNames["Depots"]
    routesLayerName = subLayerNames["Routes"]

    candidateFields = arcpy.ListFields(inOrders)
    orderFieldMappings = arcpy.na.NAClassFieldMappings(outNALayer, ordersLayerName,
                                                       False, candidateFields)
    arcpy.na.AddLocations(outNALayer, ordersLayerName, inOrders, orderFieldMappings,"")

    depotFieldMappings = arcpy.na.NAClassFieldMappings(outNALayer, depotsLayerName)
    arcpy.na.AddLocations(outNALayer, depotsLayerName, inDepots, depotFieldMappings, "")

    arcpy.na.AddLocations(outNALayer, routesLayerName, inRoutes, "", "")

    #Solve the VRP layer
    arcpy.na.Solve(outNALayer)

    #Save the solved VRP layer as a layer file on disk with relative paths
    arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"RELATIVE")

    print "Script completed successfully"

except Exception as e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "An error occured on line %i" % tb.tb_lineno
    print str(e)


texts: Clustering data on map using cartodb.createLayer?

texts:
We are using Google MAP API V3 and CartoDB.

I am trying to achieve whole UK wide properties on the map with points and I have plotted all the points using cartodb.createLayer as per its Geo coordinate its working fine but I want to cluster all these points together on the map.

How can I cluster the data on map using cartodb.createLayer?


texts:
I'm trying to create a script for a main script to reference as a set of commonly used terms. For example, we currently have in one of our main automation scripts

LAT = ['Lat', 'lat', 'LAT', 'Latitude', 'latitude', 'LATITUDE',     'SITE_EST_LATITUDE']

which is then referenced later as

for n in LAT:
with open(report_list) as f:
    reader = csv.reader(f, delimiter=',')
    for row in reader:
        for field in row:
            if field == n:
                lat = field

My company handles telecommunications and reports are regularly pulled from multiple sources which all use different terminology for referencing latitude/longitude. Ideally what we'd like to do is have our script reference another script which would contain all of these definitions (examples include latitude, longitude, state names, city names, etc...). So there is a lot of information -- thus making it a bit easier to just have a script reference the library of terms as opposed to typing them out in every new script we create.

Bottom line, I'm not super advanced with Python so I'm not sure how to create a new script that houses these terms (preferably in a CSV) and can be used as a reference in other scripts.


texts: How to generate a correct HTTP response for a pbf within AWS API Gateway?

texts:
I am attempting to implement this approach to serving vector tiles from PostGIS using a lambda function and an AWS API Gateway Method: https://www.sparkgeo.com/blog/vector-tile-server-using-postgis. And for the most part, I have it working.

But it doesn't seem to be in the correct format and not a usable vector tile. this is what I am getting back:

And this is what I receive back from https://a.tiles.mapbox.com/v4/mapbox.mapbox-terrain-v2/13/1310/3168.vector.pbf

I tried switching the content type to application/x-protobuf, but that had no effect.  I think I need to do something additional to the byte array to make it into a consumable vector tile but the blog post was light on details regarding the API Gateway.

But this is what I have setup:


texts: Turn dissolved polygons to singlepart in geopandas?

texts:
I have performed a dissolve but I want to turn it to singlepart.I have searched in the docs but haven't seen something that does it.
What do you think?

Actually, I want to make it singlepart and then create a column that will have something like id so later that I will be doing an intersection I will be using this id to remember the polygons that were as one and do a condition that keeps the max area.

The important thing here is how to turn it to singlepart and maybe how to approach the column with the id thinking.

data = dat[dat['TYPE'].str.startswith(('01'))]
dslved=data.dissolve(by=TYPE')

# here needs to be the singlepart process

interscted=gpd.overlay(dslved,file_m,how='intersection')


texts: QGIS 2.14 map rendering issue when symbolising on attributes from a join

texts:
I have a spatial table that has a join to non-spatial table. I then symbolise using the joined attribute from the non-spatial table.

If the join is non cached, the rendering on QGIS freezes. I am using a PostGIS view for the non-spatial table, but I got the same results using a CSV joined the same way.

However if I cache the join the map will render; however there is no way to refresh the join, except going into the layer properties and doing a 'fake edit' on the join.

Is this a bug or are there any particular settings/limitations on joining spatial/non-spatial tables in QGIS 2.14?

Note: I didn't have this rendering issue in 2.12 or 2.8


texts: Resolving topology errors after merging polygons

texts:
I have two SpatialPolygonDataFrames, and I have two merge them. SPDF1 contains watershed are falling within the country, SPDF2 has watershed are falling outside the country, so the both should match properly.

maptools::unionSpatialPolygons

After merging, I have several types of error.
1. I have inner rings of 0 area. I can remove them based on the assumption that I only need the first (most outer) ring of the SpatialPolygon.
2. I have "Peninsula"s also of 0 area.

How can I remove these "Peninsulas" automatically (with R code)?


texts: Convert between different GML like formats

texts:
We have two different XML formats representing logically Geography Markup Language (GML), but they are our proprietary formats. Think of it like two separate groups looked at different formats using GML and each one of them decided to modify them in their own way to build up a new format.

Now we want to be able to merge them to single format which will be identical to one of the two input formats.

Do tools exist for this task like supplying input/output samples and it will prepare the specific conversion?
If no - What are the best practices here? Because the boot-force method is just to build simple C/Java program building custom XML format.


texts: Automatically add date/time and sequential number to new GPS points in the field

texts:
We are using a tablet running QGIS 3.4.9, connected to a bluetooth GPS device, to record data in the field. We create GPS points using the "Add Point" button in the GPS panel. After tapping on this button, we can fill in attributes manually.
Is there a way to set it up so that the date/time and a sequential index number is added automatically to each new GPS point? (Preferably without having to use Python.)
For the sequential index number, we cannot use rownumber because the points are added in a seemingly random order. For example, the last point added is usually not the last row in the attribute table.


texts: Labels overlapping annotations in ArcMap 10.6

texts:
I came across this issue regarding labels overlapping annotations in ArcMap 10.6 (picture left) where they don't overlap in 10.5 (picture right)... What am I missing? I simply opened the mxd in the newer version so the settings are all the same (I would have guessed, but obviously 10.6 ignores some of them):

I use the Maplex Label Engine and already tried to add another polygon feature class serving as a mask for my annotations, giving it a feature weight of 1000 (both interior and boundary) which should block all labels from being displayed in that area. But that is simply being ignored!


texts: Get coordinates of nodes in a polyline using ogr2ogr

texts:
I'm trying to figure out how to extract coordinates of first and last node of a ESRI shape polyline using ogr2ogr.

In QGIS I use xat(n) and yat(n), so what I'm looking for is the same result using ogr2ogr in the terminal.

The reason I want this is to first calculate the angle between the points (in QGIS I use (atan((xat(0)-xat(1))/(yat(0)-yat(1)))) * 180/3.14159 + (180 *(((yat(0)-yat(1)) < 0) + (((xat(0)-xat(1)) < 0 AND (yat(0) - yat(1)) >0)*2))) and then store these values in an attribute column and turn the line object into a point object.

I'd like to do this for each feature in the shapefile.


texts: Fully automating a workspace to create maps with different file names but same information?

texts: I have a lot of maps (over 600) to create. They all contain the same column titles and the layouts will be the same. I don't want to have to go through and repeat the same process everytime, I would like to just have to open the file and have it map on its own. I have created a workspace to open the other tables I need for every map, but I need to know if it is possible to write a program or code with mapbasic to ask the user to specify what table to open, then using that finish the program by joining the tables and mapping. The first and last parts I can do. I just need to know how to create a dialog box to have the user select the table, or if it is even possible.


texts: Finding Normal to a given point in ECEF coordinates

texts:
I have a point of reference called Ground Reference Point (GRP) given in Earth Centred Earth Fixed (ECEF) coordinates as given below:

(1.4892e+6, 6.1102e+6, 1.0637e+6) (all distances in meters)

I want to find a normal to the earth ellipsoid at the above mentioned GRP again in ECEF coordinates. What is the process to calculate it?


texts: Join/Relate in ArcGIS 9.3.1 ArcEditor through a python script - is it possible?

texts: I'm working on a ArcGIS 9.3.1 ArcEditor and I want to write a python script that will do a Join or Relate. Problem is, that join table\join field are restricted to ARCINFO (higher license), and relate isn't doable by code (or so it seems from the ARCGIS help). As upgrading the license isn't an option, Does anybody have a way I can do a join between two shapefiles in ArcGIS or Opensource tools?


texts: MMQGIS:Create:Hub Distance: Output shapefile in wrong coordinate system

texts:
I have a QGIS-project that is completely in EPSG:31467 - DHDN / Gauss-Kruger zone 3. It consists of random points (camera traps) surrounded by a polygon (forest edge). I already converted the polygon into lines in order to measure the distance between each camera trap and its nearest forest edge. The tool works well apart from the problem that the output shapefile surprisingly is in EPSG:4326 - WGS 84, but still displayed correctly (connecting lines between points and outline).

After then correcting the CRS of the new shapefile (WGS84) to the original project CRS (DHDN 3) the new layer ends up far away from where it should be.

Any suggestions how to get it right?


texts: Transfer geometry of a feature class into a table

texts:
I have a table that has a field that matches a field in a feature class.  the table has M records for 1 row in the feature class.  I would like each row in the table to have the geometry of the feature class.  This would result in stacked geometries with different attributes.

If this were a 1:1 relationship, I would join the table to the FC and then export for the result I want.  if I join the FC to the table the attributes line up, but there is no geometry included in the table.

Any Ideas on a tool or workflow for this?

I am working with 10.3 standard.


texts: ST_Length (GEOM, true) for calculating distances in Geographic Reference System

texts:
I use (1) ST_Length (ST_Transform (geom, 3395)) to get distances in meters from points in Geographic Reference Systems EPSG: 4326. The thing is I found that  with  (2) ST_Length (geom, true)  you can get also the distance in meters but using the spheroid.
If I am right, is the option (2) more accurate if the GIS project use a Geographic Reference System for Global data, especially near to the poles?
ST_Length reference: http://postgis.org/docs/ST_Length.html
EPSG:3395 reference: http://spatialreference.org/ref/epsg/3395/


texts: Google Earth Engine print result for a collection

texts:
I have the next code and I would like to print the pixel value summarize of "surf_refl_01" for all the collection instead to do it one image by one.

Is it possible?

The code:

// Region of interest
var region = table2;

var collection = ee.ImageCollection('MODIS/MOD13Q1')
.filterDate('2000-01-01', '2017-12-31')
.filterBounds(table2);

print(collection);

var image1 = ee.Image('MODIS/MOD13Q1/MOD13Q1_005_2000_02_18');

var clip = image1.clip(region);

var sumB3 = clip.select('sur_refl_b01').reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: table2,
  scale: 250,
  maxPixels: 1e9
});

var sumB4 = clip.select('sur_refl_b02').reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: table2,
  scale: 250,
  maxPixels: 1e9
});

Map.addLayer(clip);

print(sumB3)
print(sumB4)


texts: Leaflet - Leaflet.MarkerCluster.LayerSupport add/remove layers using external checkboxes?

texts:
I'm trying to implement the Leaflet.MarkerCluster.LayerSupport plugin following this example, however while it seems easy to add and remove data using buttons, there is nothing here to indicate one can use the same format with an external checkbox input.

I've been trying to take apart the L.controls inside the map for using externally but it seems that's rather difficult as its built around leaflet.js code. Does anyone have any suggestions for how the same functions for Leaflet.MarkerCluster.LayerSupport called upon by a button with "data-group="2" data-op="checkOut" (for example) could be applied to an external checkbox?


texts: Are there any web services that return an elevation raster given a bounding box?

texts:
I have an application that uses raster elevation data extensively in its computations.  I would like to have a tool that allows users to define the minimum and maximum latitude-longitude coordinates (as a bounding box) and download publicly available elevation data from a web server as a single TIFF file or some other raster format that supports elevations.  I don't know whether this would be a Web Map Service or a Web Coverage Service.  I have searched online for something like this and have not been able to find anything.  Of course, I know about the National Map for US data.  I am also familiar with Google Earth and have used NASA World Wind in the past.  I am also familiar with ESRI's World Elevation Services.

Does anybody have details of a web server that does something like this?  Free would be nice, but maybe there's something available for a fee based on how much the server is used.  If not, any suggestions of how I could set up a server that does something like this myself if I download all the elevation data I need?


texts: Invalide charater error associated with "Select-Case" vba

texts:
Anyone see a problem with the following code.  I keep getting a 999999 - invalid character error.

Select Case [code]
  case 5112: x = "AP030"   ' trunk

  case 5113: x = "AP030"   ' primary

  case 5113: x = "AP030"   ' primary

  case 5114: x = "AP030"   ' secondary

  case 5115: x = "AP030"   ' tertiary

  case 5121: x = "AP030"   ' unclassified

  case 5122: x = "AP030"   ' residential

  case 5123: x = "AP050"   ' living_street

  case 5124: x = "AP050"   ' pedestrian

  case 5131: x = "AP030"   ' motoway_link

  case 5132: x = "AP030"   ' trunk_link

  case 5133: x = "AP030"   ' primary_link

  case 5134: x = "AP030"   ' secondary_link

  case 5135: x = "AP030"   ' tertiary_link

  case 5141: x = "AP030"   ' service

  case 5142: x = "AP050"   ' track

  case 5143: x = "AP010"   ' track_grade1

  case 5144: x = "AP010"   ' track_grade2

  case 5145: x = "AP010"   ' track_grade3

  case 5146: x = "AP010"   ' track_grade4

  case 5147: x = "AP010"   ' track_grades5

  case 5151: x = "AP050"   ' bridleway

  case 5152: x = "AP050"   ' cycleway

  case 5153: x = "AP050"   ' footpath

  case 5154: x = "AP050"   ' path

  case 5155: x = "AQ150"   ' steps

  case 5199: x = "AP050"   ' unkknown

end select


texts: Why openlayers show wrong measures of line from ECW geoserver?

texts:
I have loaded the following raster into geoserver and added measurement control.
The thing is that the measure gives me wrong result (it shows the car length is 49m)
Any idea what can cause this ?!

For the code , I used the generated code from "Preview Layer" in geoserver and just added the dynamic measure control to it.


texts:
I have a fairly simple leaflet project in RStudio

library(rgdal)
library(leaflet)
library(htmlwidgets)

#load data/layers
m = leaflet() %>% {add polygons....}
saveWidget(m, file="HTML_File.html", selfcontained = F)

The resulting HTML page displays nicely, but I would want to respond to a click event by adding a pop-up that provides the lat / lon of the point clicked, as well as the attribute of the top-most layer.

I have no idea how to go about it. I noticed that the leaflet.js script that gets exported has code to handle double clicks. I tried to hack that code to provide an alert, along the lines of

alert("Lat, Lon : " + e.latlng.lat + ", " + e.latlng.lng)

but clearly I don't know what I'm doing.


texts: Converting x,y on mercator to WGS84 latitude and longitude

texts: How do I convert a particular point on a Mercator map (points in a Mercator map svg) to WGS84 coordinates? Is it simply a direct translation of the x,y values, using a multiplier which matches the map width & height to -180..180 and -90..90?


texts: Seeking shapefiles for real world places

texts:
Possible Duplicate:
  Are there any free administrative boundaries available as shapefiles?

This is a follow up to my previous question How to highlight the area in google map on mouse hover. In my application i need to highlight the area/place/region when mouse hovers. @simo pointed me out to try kml shapefiles to generate the polygons on mouse over.

Is there a open wiki where  i can download the shapefile mappings of the real world boundaries.

I have searched but couldn't find any solid data.  But i find this wikimapia kml integrated google map. This allows the wikimapia places/kml shapefiles to be rendered on top of google maps. This is pretty much what i need. But instead  the whole wikimapia layer like this

i only wanted to highlight the area on mouse hover (without the explicit blue boundaries in above pic), like this

How can I do this?


texts: Altering number of ranges / classes in Rx shapefile?

texts:
Is there a way in QGIS to import a shapefile that is a variable rate dry fertilizer prescription that has a large number of rates (can be in excess of 50) and reduce the number of rates to less than 10?

For example, rates of 50 to 65 lbs/A can be reduced to 57. Rates of 66 to 80 can be reduced to 74. Rates of 81 to 95 can be reduced to 87. These are just examples, but hopefully it further details what I am trying to do. I want to export the resulting shapefile with the reduced number of rates into the Mapstar software to create a PMH file that can be imported directly into the aircraft for application.


texts:
How to create intensity image using point cloud in FME software? I already generate using  image rasterizer transformer but the quality still need to be improved. Is there any other transformer that can give better result?
I attached sample between intensity image generated from TerraScan and FME. Intensity created from FME is not very smooth compare to intensity generate from Terrascan. I also having problem with contrast balance which I think can be solved at FME.
]2


texts:
I've decided to work on GIS implementation in transportation.

I've these small tiles of images, with world files.

I want to combine them to create a base map.

I don't know coding/programming, so please suggest some other other method to combine these tiles.


texts: Passing variable from ArcPy script to ModelBuilder?

texts:
I have a model that passes a variable to a python script. The script uses that variable to query some directories and return a list of paths. What I want to know is:

How can I then use each of those paths as the parameter in a subsequent operation?
Is the bridge between python and ModelBuilder a one-way street?


texts: If a ship gave a mayday/SOS call out to a coast guard and sent the lat/lng coordinates where they needed to be rescued, is WGS84 implied or assumed?

texts: Lat/Lng coordinates don't necessarily mean anything without knowing if it's WGS84 vs. something else. If a ship/crew/captain were to send an SOS call out with just coordinates, do they have to say WGS84?


texts: I'm a landscape architecture student who is starting to learn QGIS. How would I be able to manually input a scale I want in inches to feet when exporting for print? For example, if I would like to print 1" = 100' would I be able to do that?


texts: Programmatically removing layer from QGIS composer legend without removing it from map canvas?

texts:
I am working on a QGIS Python plugin which includes automated map production.  I use a series of map templates, each of which has a main map, an overview map (with a layer set different to the main map), and a legend which is to be linked to the main map.

In QGIS 2.2 and 2.4 this is very easy to achieve:

First load the main map layers into the canvas then:

mainMap.storeCurrentLayerSet()
mainMap.setKeepLayerSet(True)
legend.setComposerMap(mainMap)

However in QGIS 2.6 this does not work - no layers are listed in the legend.  Also the legend properties panel shows the 'Map' property of the Legend as 'None' rather than 'Map 0' (this was also the case in QGIS 2.2 and 2.4).  I have tried several methods to get the legend to show the correct set of layers, all of which have failed.  My most recent approach was to load all layers for both maps into the map canvas, and then try to remove the overview map layer names and symbology from the legend:

model = legend.modelV2()
for r in range(0, model.rowCount()):
    for c in range(0, model.columnCount()):
        if model.index(r, c).data() == "Locality Map Layers":
            model.removeRows(r, 1)

legend.updateLegend()

This removes the Locality Map Layers group from the legend, but it also removes them from the main QGIS map canvas, which I do not want to do.  It seems that the map canvas and the legend are both using the same model.

Is there another way to use pyqgis to remove a layer from the composer legend without also removing it from the map canvas?


texts: Should I use slippy tile folders to display a map when having a large set of information?

texts:
I have about 1 million polygons divided into 3 ShapeFiles, each of them having 2GB, a total of 6GB. If I generate 256x256 tiles from zoom level 8 to 17, I will have 25 millions of tiles (image files) and about 200GB of folders and images.

Afterward, I can display the tiles with leaflet using Slippy Tile format: http://server.com/z/x/y.png. The frontend app would just request image files to a static bulk instead of using any Map Server (GeoServer, Mapnik or MapServer).

Is this a good strategy when it comes to a large set of spatial information as I have? Or should I use something more robust like GeoServer or mapnik?

If I have 10k users a day, will I have good performance accessing files directly?

UPDATE
I've updated the question to clarify it.


texts: Disabling/hiding customize option from ArcMap menu bar

texts:
I need disable or hide the Customize option. I have listen about a mxd file configuration, but isn't clear to me.

Is there any way?

I'm using ArcGIS 10.4


texts: Segmentation Violation in ArcView 3.1 while performing UpdateLegend?

texts:
I am loading a ".avl" legend file into ArcView 3.1 and I am plagued with a recurring Segmentation Violation Error when I execute the updatelegend command in the theme.

It looks like the legend file was successfully loaded (evidenced by the updated legend thumbnail on the theme name in the view), but the actual image in the view does not update. If I double-click on the theme and click the "apply" button in the legend editor, then it works no problem, but I can't get the code to work!

The sample code is:

theView = av.GetActiveDoc
x = theView.GetThemes.Get(0)
theLegend = x.GetLegend

chk = theLegend.Load("C:\ArcView\Legend1.avl".AsFileName , #LEGEND_LOADTYPE_ALL)

if (chk) then

  x.SetLegend(theLegend)
  x.UpdateLegend

end


texts: Using ordnance survey data to make a D3 map?

texts:
I would like to display a map of the UK Westminster Parliamentary Constituencies in D3 with a width of 980px.

I have downloaded and extracted the ordnance survey ESRI shapefiles from here: https://www.ordnancesurvey.co.uk/opendatadownload/products.html.

I have converted the .shp file to GeoJSON and then to TopoJSON:

ogr2ogr -t_srs WGS84 -f GeoJSON input.json bdline_essh_gb_0516/Data/GB/westminster_const_region.shp
topojson -o output.json -s 7e-7 --id-property NAME constituencies=input.json

The file output.json is about 110K. Here's some sample output:

{
  "type": "Topology",
  "objects": {
    "constituencies": {
      "type": "GeometryCollection",
      "geometries": [
        {
          "type": "Polygon",
          "id": "Aberavon Co Const",
          "arcs": [
            [
              0,
              1,
              2,
              3,
              4,
              5
            ]
          ]
        },
...
      ]
    }
  },
  "arcs": [
    [
      [
        4058,
        1538
      ],
      [
        0,
        6
      ],
      [
        -1,
        0
      ],

    ]
  ],
...
  "transform":{
    "scale":[
      0.09800980098009801,
      -0.18306403629989057
    ],
    "translate":[
      0,
      1830.4572989626058
    ]
  }
}

My D3 code looks like this:

var width = 980;
var height = 980;

var projection = d3.geo.albers()
.center([1.3, 55.4])
.rotate([4.4, 0])
.parallels([50, 60])
.scale(100)
.translate([width / 2, height / 2]);

var path = d3.geo.path()
.projection(projection);

var svg = d3.select("#data-vis")
.append("svg")
.attr("width", width)
.attr("height", height);

queue()
.defer(d3.json, "output.json")
.await(ready);

function ready(error, uk) {
  if (error) throw error;

  svg.append("g")
  .attr("class", "constituencies")
  .selectAll("path")
  .data(topojson.feature(uk, uk.objects.constituencies).features)
  .enter()
  .append("path")
  .attr("d", path);

}

The output looks like this:

I've tried these commands for creating my map data:

ogr2ogr -f GeoJSON input1.json bdline_essh_gb_0516/Data/GB/westminster_const_region.shp
ogr2ogr -t_srs WGS84 -f GeoJSON input2.json bdline_essh_gb_0516/Data/GB/westminster_const_region.shp
ogr2ogr -t_srs EPSG:4326 wgs84.shp bdline_essh_gb_0516/Data/GB/westminster_const_region.shp
ogr2ogr -s_srs EPSG:27700 -t_srs EPSG:4326 input4.json bdline_essh_gb_0516/Data/GB/westminster_const_region.shp
ogr2ogr -f GeoJSON input4.json wgs84.shp

...and then these:

topojson -o output1.json input1.json
topojson -o output2.json constituencies=input1.json
topojson -o output3.json -q 1e4 constituencies=input1.json
topojson -o output4.json -q 1e3 constituencies=input1.json
topojson -o output5.json -s constituencies=input1.json
topojson -o output5.json -s 7e-7 constituencies=input.json
topojson -o output6.json -s 7e-7 --id-property NAME constituencies=input1.json
topojson -o output7.json -s 7e-7 --id-property NAME --width 980 constituencies=input1.json
topojson -o output8.json -q 1e4 -s 7e-7 --id-property NAME --width 980 constituencies=input1.json
topojson -o output9.json -s 7e-7 --id-property NAME constituencies=input2.json
topojson -o output10.json --id-property NAME constituencies=input4.json

Where am I going wrong?


texts:
Using openlayers, I would like to create a temporal filter to query a GeoServer WFS service.

There is a WFS filter example here without the temporal (date) portion that i need:
http://openlayers.org/dev/examples/wfs-filter.html

My dates on GeoServer are stored in the following format: "2010-09-12 08:10:40+00"

I seek any links or code snippets that might help me?


texts: Displaying contour lines calculated from SRTM in Google Earth?

texts:
I'm relatively inexperienced in GIS, but would like to display contour lines calculated from SRTM data in Google Earth.

It doesn't matter whether they are actual linestrings or just an image overlay.

In a previous life, I worked with SRTM height files in Matlab, so I could probably botch a solution in Matlab or Scilab.

Is there a simple way to do this using free software?


texts: Installing PgPointCloud extension for PostgreSQL on Ubuntu Bionic

texts:
I'm stuck with the installation of the PgPointCloud extension for PostgreSQL. I found this article:
https://connect.boundlessgeo.com/docs/suite/4.6/dataadmin/pointcloud/install.html

But when I run the following:

apt-get install postgresql-10.4-pointcloud

I get:

Reading package lists... Done
Building dependency tree
Reading state information... Done
E: Unable to locate package postgresql-10.4-pointcloud
E: Couldn't find any package by glob 'postgresql-10.4-pointcloud'
E: Couldn't find any package by regex 'postgresql-10.4-pointcloud'

I am using Ubuntu 18.04 [Bionic] and postgreSQL 10.4.

Also, when I type:

 pip3 install glob

I get:

  Collecting glob
  Could not find a version that satisfies the requirement glob (from versions: )
  No matching distribution found for glob

Package regex I installed successfully, but I still get the same error that it can't be found when I try to install pgpointcloud extension.