texts: Determining if Layer is Point using Python Toolbox?

texts:
I have been working on a Python toolbox which contains 6 parameters. I am a little lost on all of the different data types available in ArcMap (it seems like there's hundreds of them!). However, one of the input parameters has a file extension of .shp (shapefile). I have always used the GPFeatureLayer data type and it has always worked for everything else I hope to do with my tool, although I don't completely understand why. I have only been doing my testing with point shapefiles for the time being.

What I would ultimately hope to do is check if the user has entered a shapefile that is not a point (or a collection of points). If it is not a point, throw a warning to the user, but still allow them to use the tool. From what I understand, you can check the shape type by utilizing the Describe() class. However, when I attempt to pass my parameter into Describe (in order to return a Describe object), ArcMap throws an error:

Here is how I define the parameter I am focusing on:

def getParameterInfo(self):
    """Define parameter definitions"""

    # First parameter = layer with customized features
    param0 = arcpy.Parameter(
        displayName = "Source layer with customized features",
        name = "source",
        datatype = "GPFeatureLayer",
        parameterType = "Required",
        direction = "Input")

And I was hoping to throw the warning in the updateParameters() method:

def updateParameters(self, parameters):
    if parameters[2].value == True:
        parameters[3].enabled = True
        parameters[4].enabled = True
        parameters[5].enabled = True

    else:
        parameters[3].enabled = False
        parameters[4].enabled = False
        parameters[5].enabled = False

    describe = arcpy.Describe(parameters[0])
    if describe.shapeType == 'Point':
        arcpy.AddWarning('Warning: This is a point shapefile')

Any ideas as to why I am getting this error?

Am I using the correct data type?

I have tried using other data types for this parameter, such as DEFeatureClass, GPPoint, and DEShapefile, but I received the following error if I tried adding my shapefile in as my parameter:

The error was the same regardless of which of the three I used.

My issue with the Describe() class has been cleared up by adding the .value extension in my updateParameters() method:

def updateParameters(self, parameters):
    describe = arcpy.Describe(parameters[0].value)

However, now when I attempt to use the setWarningMessage() function, I get no response from ArcMap. My updateParameters() consists of the following:

def updateParameters(self, parameters):
    describe = arcpy.Describe(parameters[0].value)
    #check = str(describe.shapeType)
    if describe.shapeType == 'Point':
        parameters[0].setWarningMessage('This is a point feature class')
    else:
        parameters[0].clearMessage()

    if parameters[2].value == True:
        parameters[3].enabled = True
        parameters[4].enabled = True
        parameters[5].enabled = True

    else:
        parameters[3].enabled = False
        parameters[4].enabled = False
        parameters[5].enabled = False

    return

When I add a shapefile that I know in fact is a point, the tool does not output a warning. It simply accepts it without any feedback to the user.

Does it matter if the Describe.shapeType function has an output of Unicode?

I tested this through the ArcMap Python window with the following sequence of commands:

However, even if I try to convert the Describe.shapeType to a string, I still don't get a warning message. I've done the following tests:

Rather than checking for: if describe.shapeType == 'Point', I tried to check for: if describe.shapeType == u'Point', which would force it to check for a Unicode data type.
Converting the shapeType into a string via str(decribe.shapeType).
Using the encode method to change the Unicode into a sequence of ASCII characters: describe.shapeType.encode('ascii', 'ignore').

Hence, I've hit a wall once again.


texts:
I want to import an elevation layer into GeoServer for later WMS access. This is the base layer I would like to import - http://en.wikipedia.org/wiki/File:Srtm_ramp2.world.21600x10800.jpg . What I'd like to do is reproject the layer to Mercator so it looks like this http://upload.wikimedia.org/wikipedia/commons/thumb/f/f4/Mercator_projection_SW.jpg/350px-Mercator_projection_SW.jpg , then somehow convert it to GeoTIFF so that I can import it into GeoServer.

How should I proceed?

I tried using gdal_warp to reproject my layer but it warps the layer completely leaving nothing behind. I had some success reprojecting the layer using G.Projector, but it doesn't work for bigger images and it doesn't georeference the image anyway (GeoServer reports no layers in GeoTIFF file).


texts: Creating network in ArcMap from an RGB image to find optimal path between locations?

texts:
I am hoping to find optimal routes from an Origin (car park) to destinations (Library, Blocks, Rec Center).

I have an RGB ortho-mosiac image taken from an UAV that I want to digitize paths, car parks, and infrastructure, and then find an optimal route from one to another.

How would I go about making sure the polygons and lines can be used to distinguish which path might be the best to choose from?

And therefore adding stops etc. using network analysis toolbar.


texts: Find coordinate of point on the boundary of a circle given a latitude/longitude and radius

texts:
I understand how to use the haversine formula to compute the distance between two points on the earth. However, I would like to answer a similar, albeit different question. How can I compute the latitude/longitude of any point on the boundary of a circle?

For instance, I have a center point at 31.782643, -102.291938 and a radius of 245km. I would like to find the latitude/longitude coordinate for any point on the edge of this circle. How can I do this?


texts: Extract attributes from line to point shp

texts: If anyone can help me it would be great. I am using ArcMap 10 and would like to extract the attributes from a line shp (Roads) to a point shp created from the line shp. The points fall on the intersections and in the mid point of the roads. I would like to populate the point attributes table with that of the roads attributes? Is there an easy way to do this as i have over a 100 000 records? For example a point at the intersection i would like to create a field (Road names) of the two roads intersecting.


texts: How do I use a generic List with my Arc Objects code?

texts:
/// <summary>
/// Gets a list of all geoms in a Feature Class
/// </summary>
/// <param name="FClass">The Feature Class</param>
/// <returns>a List of Geometries</returns>
private List<T> GetGeometries<T>(IFeatureClass FClass)
{
    List<T> List_Geometries = new List<T>();

    IQueryFilter aQF = new QueryFilter();
    aQF.WhereClause = "";
    ISelectionSet selectionSet = FClass.Select(aQF, esriSelectionType.esriSelectionTypeIDSet, esriSelectionOption.esriSelectionOptionNormal, null);
    IEnumIDs enumId = selectionSet.IDs;
    int nFeatures = selectionSet.Count;
    enumId.Reset();
    int id;
    for (int j = 0; j < nFeatures; j++)
    {
        id = enumId.Next();
        IGeometry pGeom = FClass.GetFeature(id).Shape;
        List_Geometries.Add((T)pGeom);
    }
    return List_Geometries;
}

I don't get any diagnostics errors on the above routine. However I do get an error when I attempt to call this method. How can I call this method from the following set of code?

private List<IPoint> GetPointsList(IFeatureRepository Bfr)
{
IFeatureClass Bfc = Bfr.GetFeatureClassByName("PointFeatureClass");
List<IPoint> listPoints = GetGeometries(Bfc);
return listPoints;
}

By the way, this is not the way to do it, but I'm just trying to give an example of what I'm trying to do.


texts: Calculating bearing of building using ArcGIS Desktop?

texts:
Is it possible to get the main bearing of buildings by floor plan?

I want to have a building (by ID) and it's main bearing (e.g. E/90° or W/270°). The idea is to change the polygon-feature class into a line-feature class. Then I can calculate the lengths of each side. I can also calculate the bearing of each side. Buildings with a more complex floor plan I planned to cumulate the side's lengths on the bearing.

The problem is that one building with four parallel sides get four different bearings because of the starting point for the bearing-calculation (see image). Furthermore the calculation doesn’t always start at the same end of a line.


texts: DEM of difference - estimate the error - simple method

texts:
I have two DEMs of the same area, measured at two different dates. Subtracting one from the other, I can create a DEM of difference (DoD) to find out about the volumetric change in m³.

DEM 1 has a resolution of 0.05 m with an uncertainty (RMSE error) of +/- 0.03 m.
DEM 2 has a resolution of 1.00 m with an uncertainty of +/- 0.10 m.

The volume that I calculated between both DEMs is 150.000 m³.

How can I include the uncertainty of both DEMs in the resulting volume? I would like to state that there was a volumetric change of 150.000 m³ +/- xy m³ or a volumetric change of 140.000 - 160.000 m³.

I have already been searching for error propagation methods, but I am not very into the topic and I often do not understand how to do it. Any ideas?


texts: QGIS 3.4 and loading raster from multiple variable NetCDF file

texts:
Using QGIS version 3.4.4 when loading a NetCDF file out.nc which contains multiple variables on a fresh Ubuntu 18.04 install I encounter the following error:

CRITICAL    Invalid Layer : Raster layer Provider is not valid (provider: gdal, URI: out.nc

The same file load correctly on another 18.04 Ubuntu environment with more package setup. Unfortunately I cannot find which package is making the difference. The libgdal package seems to be the same on both environments.


texts:
I am very new with OWS, specially WPS, and I am testing 52North WPS for a project. I have run the Sextante visibility process, which result format is a tiff image.

The execute response starts as follows:

<?xml version="1.0" encoding="UTF-8"?>
<ns:ExecuteResponse xmlns:ns="http://www.opengis.net/wps/1.0.0"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsExecute_response.xsd"
    serviceInstance="http://localhost:8081/wps/WebProcessingService?REQUEST=GetCapabilities&amp;SERVICE=WPS"
    xml:lang="en-US" service="WPS" version="1.0.0">
    <ns:Process ns:processVersion="1.0.0">
        <ns1:Identifier xmlns:ns1="http://www.opengis.net/ows/1.1">visibility</ns1:Identifier>
        <ns1:Title xmlns:ns1="http://www.opengis.net/ows/1.1">Visibilidad</ns1:Title>
    </ns:Process>
    <ns:Status creationTime="2013-04-17T12:13:16.298+02:00">
        <ns:ProcessSucceeded>Process successful</ns:ProcessSucceeded>
    </ns:Status>
    <ns:ProcessOutputs>
        <ns:Output>
            <ns1:Identifier xmlns:ns1="http://www.opengis.net/ows/1.1">RESULT</ns1:Identifier>
            <ns1:Title xmlns:ns1="http://www.opengis.net/ows/1.1">Visibilidad</ns1:Title>
            <ns:Data>
                <ns:ComplexData mimeType="image/tiff"><![CDATA[MM?*?????????????????????

To extract the image to a file (to store it in Geoserver) I get the CDATA, decode with Base64Decoder to a byte array and the write it to a file, like visibility.tiff.

But this image seems to be "corrupted", I can't preview it and Geoserver fails when I try to create the datastore with the error: "Could not list layers for this store, an error occurred retrieving them: El argumento "value" no debe ser nulo." The last sentence means the argument "value" must not be null.

Am I doing something wrong?


texts:
I opened a headerless flat binary file in ENVI.

I used menu->save file as->TIFF/GeoTIFF to save it as tiff image.

I have 100s of images to convert into tiff.

How can I do this batch processing in ENVI?

In other words, for batch processing in ENVI is there a GUI or command line?


texts: When the bearing of a line is calculated using the "Add Geometry Attribute (LINE_BEARING)" tool, is the bearing referenced to Grid North or True North?


texts: Historical maps as cylindrical projections

texts:
Objective: I want to overlay old, historical maps on top of Google Maps. According to Wikipedia, Google Maps uses a variant of the Mercator projection, a close relative of cylindrical projections.

Research: I have found various historical maps sources, particularly the website by David Rumsey. However, all these maps (as well as others I've found online) are conical projections.

Questions: Here are three scenarios to achieve my objective.

Are there historical maps available as cylindrical projections?
If not, then are there tools that can 'rubber-sheet' existing conical projections into cylindrical projections?
If neither is the case, then what is the best way to 'pin up' a conical projection onto a cylindrical one?


texts: I have digitized the internal of a building's floors using ArcMap 10.2.2 and need to georeference the shapefiles now. The only raster file I have, is a google image file which isn't accurate at all. Any suggestions on how can I georeference the shapefiles?


texts: Select and copy attributes features in ArcMap using Python add-in tool?

texts:
I'm trying to add a tool to an add-in toolbar to select features from an existing feature class and copy attributes to another feature class, from point to point.

TopoHeight_Field_Name='Comment'
def testing(mxd,Transfare_Anno_row):
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name=="sManhole":
           msg= "Your Layer is " + str(lyr)
           with arcpy.da.UpdateCursor(lyr,TopoHeight_Field_Name) as cursor:
              pythonaddins.MessageBox(msg, "My Layer")
              for row in cursor:
                   pythonaddins.MessageBox(msg, "My Layer")
                   row[0]=Transfare_Anno_row
                   cursor.updateRow(cursor)

class SelectByLine(object):
"""Implementation for SelectByLine_addin.tool (Tool)"""
def __init__(self):
    self.enabled = True
    self.shape = 3
def onMouseDownMap(self, x, y, button, shift):
    mxd = arcpy.mapping.MapDocument("CURRENT")
    pointGeom = arcpy.PointGeometry(arcpy.Point(x, y), mxd.activeDataFrame.spatialReference)
    searchdistance = getSearchDistanceInches(mxd.activeDataFrame.scale)
    lyr = arcpy.mapping.ListLayers(mxd)[0] # assumes you want to select features from 1st layer in TOC
    arcpy.SelectLayerByLocation_management(lyr, "INTERSECT", pointGeom, "%d INCHES" % searchdistance)
    with arcpy.da.SearchCursor(lyr,TextString_Field_Name) as cursor:
        for Anoorow in cursor:
            Anno_row=Anoorow[0]
            SourceFlag=1
    message = "Your mouse clicked:" + str(x) + ", " + str(y)+"Your Layer is" + str(lyr)+"Your Anno Filed is" + str(Anno_row)
    pythonaddins.MessageBox(message, "My Coordinates")
    if SourceFlag:
        testing(mxd,Anno_row)

When I Call the testing function it will send the first message box but can not send the second message box and is not updating the row.


texts: ModelBuilder adds numerous file extensions to outputs as it runs?

texts:
Using ArcMap 10.3 I am trying to create a model that will take all the shapefiles in a folder and convert them to KML files. I first create an inline variable for the shapefiles in the folder, then when I convert them to feature classes in a geodatabase (because the layer to kml tool needs layer files as input), my output ends up at %variable%.shp.shp.shp.shp... etc. Then it fails on the next tool, which is converting those feature classes to layer files. I have attached two screenshots.

Does anyone know why I might be getting multiple .shp extensions when converting shapefiles to feature classes?


texts:
Can someone help me with this error concerning installation of geopandas even after installing 'shapely'?

  Error:
           cwd: C:\Users\JASHVA~1\AppData\Local\Temp\pip-install-ccci6ax2\fiona\
      Complete output (1 lines):
      A GDAL API version must be specified. Provide a path to gdal-config using a GDAL_CONFIG environment variable or use a
  GDAL_VERSION environment variable.
      ---------------------------------------- ERROR: Command errored out with exit status 1: python setup.py egg_info Check th e logs for
  full command output.


texts:
QGIS 1.8.0-1 does not launch (nor its companions: Browser, Grass, etc.) neither on machine with Win XP SP3 nor with Win 7 (GRASS installed separately works well). Console displays a message on launch: "'findstr' is not recognized as an internal or external operable program or batch file".
I tried all solutions suggested here earlier with no result.
Could anybody advice I'd be grateful.
Ian.


texts: Is there an effective way to display integrate 50+ old 1:50,000 maps into public mapping tools such as Google Maps?

texts:
Some context:

These are Dutch colonial maps of Indonesia from about 1910-1930
At the time there were limited roads and often villages were in the middle of the jungle and have since been abandoned.
Under Indonesian law, the primary administrative concept is the 'administrative village' which does not necessarily correspond at all to any physical village, and can be several hundred square kilometres comprising and splitting up multiple geographical places. For navigational purposes the administrative villages are not very useful for this reason. Often they have further administrative subdivisions, but these are again more regions than specific places.
For historical analysis, land ownership, and place naming, it is useful to overlay historical maps onto both satellite imagery and onto modern street maps using a tool such as Google Maps, OpenStreetMap or whatever.

The question is, given that I have 54 such maps each covering around 20,000 square kilometres in total, which are scans of map sheets currently in TIFF format, is there a good way to make the maps publicly accessible using some web tool such as Google Maps so that it is not necessary to download the images, etc.?
Consider that:

The scans are each 80 MB tiffs; I guess this could be reduced using a jpeg to 20MB without too much trouble.
The scans are sheet maps with some white space at their border, also each with their own key and legend. While it would be possible to crop the white space and key off, given that not all maps are of the same date, etc., and they do not match up 100% perfectly with each other, it's probably best to have each map as a separate object overlaid onto an underlying base map, each with adjustable transparency, so that you can compare the historical mapping with the current satellite imagery and when you scroll between maps fade out one map.

Possibly Google Earth is a good tool, but I'm not sure if the maps can be hosted remotely and all 50 made persistent and accessible via a URL?


texts: Using loops within an equation for Python and inputting values into QGIS

texts:
I would like to incorporate the use of loops within my Python coding equation. I am looking to find the value for the radius (distance) for multiple earthquakes and their corresponding intensity values. I am getting the magnitude and depth data from the USGS Excel Spreadsheet file. At the moment my code looks like this (see attached) but I would like to remove the Input MMI Intensity Value code and use a loop that runs intensity values from 0 to 12, in 0.5 intervals, into the equation.
ml = float(input("Enter Magnitude Value: "))
h = float(input("Enter Depth Value: "))
Iso = float(input("Enter Intensity MMI Value: "))
import math

Io = 1.5 * (ml - math.log(h) + 1.4)
a = (Io - Iso) / (1.8)
d = (10 / 8) * (a)

d = d - 1
if d <= 0:
    print("Distance is zero")
else:
    d = math.sqrt(h * h + d * d)
    print("Radius (km): ", d)

Once I have the distance for each individual intensity value, I would like to input these distances into QGIS, along with other Earthquake data, to create an earthquake hazard map that has radii rings of varying intensity around an earthquake point.


texts: Why does this python OGR code for creating a layer and then adding a field not work?

texts:
Spot the difference!
So this works...

fp = "C:\\LiDAR\\MyLayer.shp"
layer = "MyLayer"
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource(fp)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = data_source.CreateLayer(layer, srs, ogr.wkbPoint)

field_name = ogr.FieldDefn("Ben", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)
layer = None

But this does not (Error 6: CreateField() not supported by this layer).  Anyone know why and know what I need to do to fix? Thanks.

def createLayer(fp,layer):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(fp)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    return data_source.CreateLayer(layer, srs, ogr.wkbPoint)

layer = createLayer("C:\\LiDAR\\MyLayer.shp","MyLayer")

field_name = ogr.FieldDefn("Ben", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)
layer = None


texts: Plotting Street Map for City - Where do I start?

texts:
Please excuse my ignorance of this and if I break any of the rules please be gentle as this is my first post here and I'm not GIS expert either. I need people to steer me in right direction so that I can get started.

I'm experienced programmer working with Python, C/C++, Java, and PHP to name few and have a task at hand to accomplish. We want to make an application that necessitates accurate mapping of the city. And so We need to establish accurate mapping. Collecting positions is not supposed to pose major problem as We can do it with mobile phone (correct me if am wrong here)

Now the problem is what is really needed to plot it as google map does? A street where one can zoom to see roads and buildings (only roads and buildings matter everything else is ignored)

Is there any tutorial or video that explains how I can achieve that? I looked at postGIS and looks interesting but am afraid of putting too much effort on wrong tools.

Sorry for long post, you can leave a comment to help me ask better.


texts: Algorithm(s) for Sensor Footprint at limbs

texts:
Is anyone aware of algorithm(s) that have been developed for computing remotely sensed image footprints that handle body limbs? The limb is the edge of the body / atmosphere where only the blackness of space is then visible.  Images with a visible limb are captured, possibly, when the sensor is highly oblique or at great distance from the body.

For the vast majority of remotely sensed data, the image footprint (see here for an example of Landsat footprints over texas) is roughly a square (for a framing camera) or a long rectangle (for a line scan camera).

An approach to generate a footprint for images without a visible limb can simply walk the image boundary and compute the body intersection at some node interval. For example, the upper left corner intersects the body at (x,y,z in some reference frame), the 100th pixel along the top edge intersection at (x1,y1,z1), etc. Those intersections define the footprint of the image on the body. In the case of a limb (or the terminator honestly), walking the edge is not possible because a body intersection does not exist for some subset of pixels (because we see space). A naive approach would be to create a footprint with a node at every pixel, but that could be > 70,000 geometries for some long line scan images.

Before I go building an algorithm, I wonder if one might exist.


texts: Using Spatial Join/Merge/Table View in ModelBuilder

texts: I have developed a model that spatially joins a series of line segments to line and point segments using the "One to One", "Closest" function and creating a "distanceTo" field in feet. The model then merges and creates a table view however my output seems to be duplicated multiple times. The model for the spatial joins is using the "in_memory\layer" function with the merged layer being a parameter. I would like the output to be just one record per unique id (Priority_ID).


texts:
I am new to PyQGIS.
I have a Geodatabase that contains feature classes.
I am trying to get those feature classes into PyQGIS but how do I loop a Geodatabase?


texts: CGIAR DEM + gdalwarp -- projection to Cartesian is not happening

texts:
I have downloaded a 30-by-30 tile from

http://srtm.csi.cgiar.org/srtmdata/

Then gdalinfo on the obtained single *.tif-file would yield

Driver: GTiff/GeoTIFF
Files: cut_n30e000.tif
       cut_n30e000.tif.aux.xml
Size is 36024, 36023
Coordinate System is `'
Origin = (-0.010416618510297,60.009584059976817)
Pixel Size = (0.000833333353512,-0.000833333353512)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.0104166,  60.0095841)
Lower Left  (  -0.0104166,  29.9904167)
Upper Right (  30.0095841,  60.0095841)
Lower Right (  30.0095841,  29.9904167)
Center      (  14.9995837,  45.0000004)
Band 1 Block=36024x1 Type=Int16, ColorInterp=Gray
  Min=-444.000 Max=4783.000
  Minimum=-444.000, Maximum=4783.000, Mean=348.628, StdDev=398.861
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=4783
    STATISTICS_MEAN=348.62818252965
    STATISTICS_MINIMUM=-444
    STATISTICS_STDDEV=398.86067773429

So far so good. So I would try to project this onto a 2D plane, as follows:
gdalwarp -t_srs '+proj=tmerc +lat_0=44.00 +lon_0=15.00 +units=m +no_defs' cut_n30e000.tif result_europe.tmerc.tif.
After this operation, result_europe.tmerc.tif is supposed to be represented
with respect to a Cartesian system of coordinates.
But I was surprised to find out that the lat/lon system persists:

Driver: GTiff/GeoTIFF
Files: result_europe.tmerc.tif
Size is 36024, 36023
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["unknown",
            SPHEROID["WGS84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",44],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-0.010416618510297,60.009584059976817)
Pixel Size = (0.000833333353512,-0.000833333353512)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -0.0104166,  60.0095841) ( 15d 0' 0.00"E, 44d 0' 1.94"N)
Lower Left  (  -0.0104166,  29.9904167) ( 15d 0' 0.00"E, 44d 0' 0.97"N)
Upper Right (  30.0095841,  60.0095841) ( 15d 0' 1.35"E, 44d 0' 1.94"N)
Lower Right (  30.0095841,  29.9904167) ( 15d 0' 1.35"E, 44d 0' 0.97"N)
Center      (  14.9995837,  45.0000004) ( 15d 0' 0.67"E, 44d 0' 1.46"N)
Band 1 Block=36024x1 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

Why the original lat/lon bounds are still there? I.e. what was wrong with my gdalwarp command above?


texts: How to recognise building type from polygons?

texts:
I have a number of building polygons which are extracted from the aerial images, along with their postal address. (The aerial images are not available in my case.) For example, the polygon of a square shape bungalow is represented as:

[(1,1), (2,1), (1,2), (2,2)]

There are several types of building such as house, flat, commercial and bungalow.

What I am currently doing is running shape recognition using Python. However, since most of the building polygons are square or rectangular, the recognition accuracy is very low. So I'm wondering if there is any open source packages that I can use? Or if there is any external data source that can lookup the building type with a given address?


texts: I need to create a fishnet grid (presence and absences) form species distribution areas (polygon shapes). But I need to define that presence (value 1) is correct when more than 30% of the grid is covered by the species distribution area. Any ideas of how to do this?


texts:
I am attempting to style Ordnance Survey Vectormap Local on QGIS 1.8.0 by using the styling sheets suplied by Ordnance Survey. The styling sheets are available from the OS website (http://www.ordnancesurvey.co.uk/business-and-government/help-and-support/products/vectormap-local.html). However when I load one of the Vector layers - take 'road centrelines' as an example and then load an SLD (OS VectorMap Local - Road_Centre_Lines.sld) in the style tab, it only displays two rules which don't style the layer. It's as though QGIS doesn't read the SLD correctly? Or am I not loading/applying the SLD files the right way? I was expecting the SLD to configure the symbology so that Motorways would be displayed as thick blue lines for example and minor roads to be displayed as small white lines.

I have 16 Vector layers in total and it would be a very long process to manually style all of them so I would rather avoid that approach at all cost. I would be grateful for any suggestions.


texts: Retrieving start and end point coordinates with ArcPy?

texts:
How do I retrieve start and end point coordinates with ArcPy for a polyline feature class?

I expect to pass a segment identifier to a subroutine and have it pass back start and end coordinates. The Field Calculator method doesn't work for me, because I need the value for other calculations that can't be performed within it. (I also prefer not to change the data to store these coordinates as attributes.) I am attempting to calculate address ranging for a "center out" addressing scheme. The address value depends on the distance to the "county center".


texts:
I have just drawn a set of bar charts over a layer in ArcGIS 10.2 for Desktop.

How can I manage rotate the charts keeping my dataframe with the same rotation?

The options offered in Symbology Display only lets you to display the chart as bar or column.

Any chance to set a rotation angle??


texts: Creating polygon in turf surrounding route generated by OSRM?

texts:
I am generating a route using the osrm rendering it by leaflet, and I was trying to create a  buffer automatically surrounding this route

As you can see in the previous image, there are parts of the route that are not inside the buffer generated by turf.

Does anyone know some router that can generate me more points of the route so the buffer can be more precise?


texts: Add a filter Filter to SQL Query based on geometry generated on-the-fly

texts:
I am creating a LINESTRING geometry type that is based on two point's geometries. To get the two different points, I am joining on a to and from id field.

I have an sql query like this:

SELECT geometry::STGeomFromText(
CONCAT('LINESTRING(',
(SELECT CONCAT(fstruc.Shape.STX, ' ', fstruc.Shape.STY)),
',',
(SELECT CONCAT(tstruc.Shape.STX, ' ', tstruc.Shape.STY)),
')')
, 103772) geom,

FROM dbo.SAN_TELEVISE vid

JOIN dbo.san_struc fstruc
ON fstruc.cid = vid.f_struc_cid

JOIN dbo.san_struc tstruc
ON tstruc.cid = vid.t_struc_cid

WHERE fstruc.cid IS NOT NULL
AND tstruc.cid IS NOT NULL
-- THIS DOESN'T WORK: AND geom.STIsValid() = 1

The query works, I am getting a pretty good result, but there are some bad data rows here though, and ArcMap is choking on some line strings that only have one point. So I want to apply a filter to the query to eliminate these bad rows before ArcMap gets them.

Something like geom.STIsValid ( )  = 1 would do it, but I can't get my sql query to run if I add this to the where clause. I get the error:
Cannot find either column "geom" or the user-defined function or aggregate "geom.STIsValid", or the name is ambiguous


texts: 2010 census population data by block groups for Florida

texts:
I would like to map population density by block group for the state of Florida, but I'm having a hard time trying to locate the correct data on the Census website.  I know I'll probably have to calculate the density myself, but I can't find just plain total population by census block groups.  I was going through the suggestions on this page, even though the question was about Census Blocks, but still not having luck:

Where to get 2010 Census Block data?

Do I have to download each county's census block group data separately or am I just missing something?  Their website is very cumbersome.

I have the shapefile of the block groups just can't locate the 2010 census population data pertaining by block group so I can join it with the shapefile in ArcGIS.


texts: Find end of arc given start/end angles, start point and radius

texts:
Using PostGIS, how do I work out the end point of an arc given the start point, start and end bearings (relative to true north) and radius of the arc? Eg. say I need draw an arc of 4000m radius from bearing 180 degrees to 90 degrees with the start (the 180 degree point) at POINT(10 20). How do I calculate the end point (other than approximately drawing the arc, point by point)?

This is only for short distances, so working on a plane should be OK. I can work out the distance between the start and end point using the length of chord formula and if I could just get the bearing between them as well then I could use ST_Project - but I don't know how to work out the bearing. (In the simple example above it should be 45 degrees.)

Solution:

It ended up being quite simple in the end. Thanks to whuber for explaining it very well. I wanted to put the PostGIS code here in case anyone else wants it:

centre_point := ST_Project(start_point, radius_m, normalise_angle(start_angle + pi() / 2 * CASE WHEN is_clockwise THEN 1 ELSE -1 END));
end_point := ST_Project(centre_point, radius_m, normalise_angle(end_angle + pi() / 2 * CASE WHEN is_clockwise THEN -1 ELSE 1 END));

normalise_angle is a simple function that ensures the angle is between 0 and 2PI (as required for ST_Project), start_point is a geography(Point) and radius_m is the radius in metres. Angles are in radians.


texts: Remove Color Attribute Table - GDAL - geoTiff

texts:
How can I remove the GDALRasterAttributeTable with GDAL? It interferes with styling in Geoserver.

GDALInfo output below:

Driver: GTiff/GeoTIFF
Files: gcev3-product1-cropland-extent-australia.tif
       gcev3-product1-cropland-extent-australia.tif.aux.xml
Size is 155316, 125062
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (112.509138433435837,-10.116018245998344)
Pixel Size = (0.000269494585236,-0.000269494585236)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 112.5091384, -10.1160182) (112d30'32.90"E, 10d 6'57.67"S)
Lower Left  ( 112.5091384, -43.8195501) (112d30'32.90"E, 43d49'10.38"S)
Upper Right ( 154.3659594, -10.1160182) (154d21'57.45"E, 10d 6'57.67"S)
Lower Right ( 154.3659594, -43.8195501) (154d21'57.45"E, 43d49'10.38"S)
Center      ( 133.4375489, -26.9677842) (133d26'15.18"E, 26d58' 4.02"S)
Band 1 Block=155316x1 Type=Byte, ColorInterp=Palette
  Description = Layer_1
  Min=0.000 Max=3.000
  Minimum=0.000, Maximum=3.000, Mean=1.411, StdDev=1.467
  Metadata:
    LAYER_TYPE=thematic
    STATISTICS_HISTOBINVALUES=-2147483648|463732937|403592142|-2147483648|
    STATISTICS_HISTOMAX=3
    STATISTICS_HISTOMIN=0
    STATISTICS_HISTONUMBINS=4
    STATISTICS_MAXIMUM=3
    STATISTICS_MEAN=1.4106150300956
    STATISTICS_MEDIAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_MODE=0
    STATISTICS_SKIPFACTORX=1
    STATISTICS_SKIPFACTORY=1
    STATISTICS_STDDEV=1.4672104301177
  Image Structure Metadata:
    NBITS=2
  Color Table (RGB with 4 entries)
    0: 0,0,0,255
    1: 255,255,0,255
    2: 0,255,255,255
    3: 46,48,71,255
<GDALRasterAttributeTable Row0Min="0" BinSize="1">
  <FieldDefn index="0">
    <Name>Histogram</Name>
    <Type>1</Type>
    <Usage>1</Usage>
  </FieldDefn>
  <FieldDefn index="1">
    <Name>Red</Name>
    <Type>0</Type>
    <Usage>6</Usage>
  </FieldDefn>
  <FieldDefn index="2">
    <Name>Green</Name>
    <Type>0</Type>
    <Usage>7</Usage>
  </FieldDefn>
  <FieldDefn index="3">
    <Name>Blue</Name>
    <Type>0</Type>
    <Usage>8</Usage>
  </FieldDefn>
  <FieldDefn index="4">
    <Name>Opacity</Name>
    <Type>0</Type>
    <Usage>9</Usage>
  </FieldDefn>
  <FieldDefn index="5">
    <Name>Class_Names</Name>
    <Type>2</Type>
    <Usage>2</Usage>
  </FieldDefn>
  <Row index="0">
    <F>9847120537</F>
    <F>0</F>
    <F>0</F>
    <F>0</F>
    <F>0</F>
    <F></F>
  </Row>
  <Row index="1">
    <F>463732937</F>
    <F>255</F>
    <F>255</F>
    <F>0</F>
    <F>255</F>
    <F>1. Croplands</F>
  </Row>
  <Row index="2">
    <F>403592142</F>
    <F>0</F>
    <F>255</F>
    <F>255</F>
    <F>255</F>
    <F>2. Pastures</F>
  </Row>
  <Row index="3">
    <F>8709683976</F>
    <F>46</F>
    <F>48</F>
    <F>71</F>
    <F>255</F>
    <F>3. Non croplands</F>
  </Row>
</GDALRasterAttributeTable>


texts: Zoom to Postcode QGIS Plugin - "Postcode not Found"

texts:
I'm trying to get the Zoom to Postcode plugin for QGIS to work on my machine (Mavericks, QGIS 2.2.0)

The plugin downloads fine, and when I first try to search for a postcode it alerts me that it needs to download the postcode data, which I approve.

The loading dialog box appears and seems to complete OK.

Has anyone else had this problem? Maybe I can download the file someplace else and manually include it.


texts: Generate random points in polygon. Algorithm question

texts:
i created a Random Point Generator in python using the algorithm described here by @whuber. The algorithm is implemented using Shapely, PostGIS and psycopg2. To boost performance pythons multiprocessing and threading libraries are used.  The source code can be found in my GitHub Repository. The code generates a total of 57.990.970 Points within 6388 Polygons in ~3517 sec on my Xeon-1231v3 with 16GB Ram.

Procedure SimpleRandomSample(P:Polygon, N:Integer) {
    U = Sorted list of N independent uniform values between 0 and 1
    Return SRS(P, BoundingBox(P), U)
}

Procedure SRS(P:Polygon, B:Rectangle, U:List) {
    N = Length(U)
    If (N == 0) {Return empty list}
    aP = Area(P)
    If (aP <= 0) {
        Error("Cannot sample degenerate polygons.")
        Return empty list
    }
    t = 2
    If (aP*t < Area(B)) {
        # Cut P into pieces
        If (Width(B) > Height(B)) {
            B1 = Left(B); B2 = Right(B)
        } Else {
            B1 = Bottom(B); B2 = Top(B)
        }
        P1 = Clip(P, B1); P2 = Clip(P, B2)
        K = Search(U, Area(P1) / aP)
        V = Concatenate( SRS(P1, B1, U[1::K]), SRS(P2, B2, U[K+1::N]) )
    } Else {
        # Sample P
        V = empty list
        maxIter = t*N + 5*Sqrt(t*N)
        While(Length(V) < N and maxIter > 0) {
            Decrement maxIter
            Q = RandomPoint(B)
            If (Q In P) {Append Q to V}
        }
       If (Length(V) < N) {
            Error("Too many iterations.")
       }
    }
    Return V
}

But concerning the algorithm I have an open question. It mentions a U = Sorted list of N independent uniform values between 0 and 1 which i could build using numpy.uniform(0, 1, numberOfValues). But then the algorithm will generate weird results, like one part of the polygon will not have any points in it. I work around this issue by weighting the points to be created by the recursive call using k = int(round(n * (p1.area / polygon.area))).

So my question is why this List of uniform values between 0 and 1 should be used in the first place.


texts: How can I label the coordinate position of the grid coordinate lines

texts:
I need to display the coordinates of the grid lines. I saw this article 'Fun with Graticules', which has much common with that I want to do.

But I still have a question. How did you create the new instance of the CoordinateFormatter class in GeoTools?

And also, how did you overload this setAttributes function? -

GridFeatureBuilder grid_builder =
       new GridFeatureBuilder(GridType) {
            @Override
            public void setAttributes(GridElement ge, Map<String, Object> map) {
                 some code...
            }
};


texts: ArcSDE PostgreSQL mixed environment error connecting to schemas

texts:
I have an enabled ArcSDE enterprise geodatabase with PostgreSQL 9.4 and ArcGIS 10.5.1. I am the database administrator and have only 2 other users who I work with use these SDE/DB connections. I personally use PostGIS but my coworkers use strictly ArcGIS tools. The database I set up has a sde/postgis mixed environment see here PostGIS and ArcSDE mixed environment

So I will probably get scolded for attempting to setup my DB design like I did, but I tried to do an unconventional DB/arcsde setup for easier viewing.

Because I will have hundreds of tables in this environment I didn't want a massive laundry list of every single table when I open a database connection. So instead of storing all the spatial tables in the sde schema and assigning roles and groups with different privileges, I attempted connect to specific PostgreSQL schemas instead of the standard sde schema. so every schema has name with that specific data in it, and a unique username and password.

so for example the base_layers schema contains only the base_layers data and when you sign in with the base_layers user name and password you only get the data in that schema not a hundreds of layers

so this all was working for me when it was just me connecting and using this database (which is hosted on my computer not a server).

QUESTION

now my other coworkers have started using these database connections and initially they worked fine for my other colleagues and then out of nowhere all my database connections stopped working (except when I connected with the sde user and to view all the tables). the error I would get is

what is even stranger is if I go to add query layer and sign in with one of the database schemas (that was not working when I tried to use it from add database connection) it worked and I am able to select any table I want to bring in...

no now nobody can connect specifically to those schemas, we can only use the sde account

UPDATE

from pauldzy response I ran this the query below and it returned a row in the sde_table_registry which points to a table whos name has changed. I then deleted this row and so far all the schema log in names that were not working have worked. Will update when I have tried all of them


texts: How to create hibernate layers at runtime?

texts:
We need to dynamically create hibernate layers from shapefiles. As far as we understand, the workflow is following:

Import shapefile to database(creating new table for each shapefile?)
Is it good approach?
Dynamically create hibernate entity from the table
Make the Entity visible for the  Hibernate mechanism
Create Hibernate layer in runtime

Are we moving in the right way? How would you handle our problem?

Is it good approach to creating new table for each shapefile?
If it so, how can resolve steps 2-4.


texts: Scaling map for certain layers using CARTO?

texts:
I would like layer setup in CARTO - scaling map for certain layers.

How do I change the zoom, for example, on map markers, so that they are not displayed at all scales?


texts: Calculating average raster to point distances within census tracts

texts:
I have some experience using ArcGIS, but mostly with the field calculator in the attribute table/spatial joins, so I apologize if I say something that doesn't make sense.

To provide a bit of background, I have a vector layer containing census tract data. I have another layer with geocoded supermarket locations. I'm trying to determine the number of supermarkets per census tract, which I can then use for various regression models. When I first tackled this issue I simply used raw counts for each census tract. Those counts, however, assume that residents from one census tract have absolutely no access to supermarkets in other census tracts. I've been struggling with a way to account for this in my regression model, so I decided to ditch the regression statistics for now and find a way to account for this issue in my count data.

After doing some searching, I found a method described by Liam Downey (2003) that would seem appropriate for what I'm trying to do. I'm not sure if we're allowed to directly quote academic articles, so I'll paraphrase his method in steps.

Convert census tract maps to raster maps (he used 25x25m grid cells).
Determine the number of facilities within a 1 km radius of the center of each grid cell.
Calculate the distance (m) from the center of each grid cell to the nearest facility.
Sum the values in (2) and (3) for each tract and divide by the number of cells in each tract.

You're then left with two new variables: the average number of facilities in each tract, and the average distance to the nearest facility for each tract.

This method sounds interesting, but beyond converting the census tract vector map to a raster map using polygon to raster (and truth be told, I'm not even sure I did THAT correctly), I'm not really sure how to proceed. What tools should I be using for this?


texts: Soft coding for loop to change file output of geoprocessing tool to create identical shapefiles for multiple years

texts:
I would like to perform the arcpy.CreateFeatureclass_management tool to copy a template shapefile several times to create the same shapefile for multiple years.
I would like the output files to be "neb_precip_2008.shp","neb_precip_2009.shp","neb_precip_2010.shp"... up until 2016. I understand it should be looped through, but I want to learn how to do this without a list. I'm still in the beginning stage of python learning.

try:
  immport arcpy
  from arcpy import env
  env.workspace = r"C:\Users\Documents"

  out_path = r"C:\Users\Documents"
  out_name = "neb_precip_2008.shp"
  geometry_type = "POINT"
  template = "practice_dataPrecip2008Readings.shp"
  has_m = "DISABLED"
  has_z = "DISABLED"

spatial_reference = arcpy.Describe("C:\Users\Documents\practice_dataPrecip2008Readings.shp").spatialReference

arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference)

except:
  print "Script failed to complete"
  print arcpy.GetMessages(2)


texts: CRS undefined when copying to new layer, QGIS 2.18.12

texts:
I'm using QGIS 2.18.12 and have this problem:

When I copy objects and paste them into a new layer, I get the message CRS undefined. (I'm using EPSG:28992, Amersfoort/RD New as CRS). In the long term release 2.14 version, I get the same message, but I get a window to change the CRS. In the 2.18.12 version, there is no window to change the CRS, actually, and this is pretty annoying, QGIS stops the whole action. In short, I can't paste objects anymore into new layers.

What's wrong?


texts:
I am working on a project on islands in inland waters in Europe.
In the first step I was sorting to large point layers which I got from a state office, trying to identify only points lying on such islands. I used a small Python-Script in QGIS to alter large OSM layers mapping inland waters for that, which worked quite good. But since OSM-data contain some errors, I would also like to try other approaches. Also, for comparison and some statistics I would like to identify all other islands, for which I did not get any points.

I thought, there might be a way to extract islands from the free EU-DEM dataset (see https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem), but I am not sure how to start on that.

So the question for me is: Is there a way to extract islands from DEMs like the EU-DEM?


texts: GeoPandas not projecting from 4326 to a UTM zone properly?

texts:
I am at my wit's end trying to figure out why my data is not being reprojected properly in a geodataframe.

Given this input data, saved in a variable called power_lines:

{'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'properties': {}, 'geometry': {'type': 'LineString', 'coordinates': [[34.059162, -18.676658], [34.193058, -18.677471]]}}]}

I create a geodataframe like this, specifying that the data is in ÈPSG:4326:

power_line_gdf = gpd.GeoDataFrame.from_features(power_lines["features"], crs="epsg:4326")

When I save this out to file and look at it, it looks like this:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "LineString", "coordinates": [ [ 34.059162, -18.676658 ], [ 34.193058, -18.677471 ] ] } }
]
}

I then reproject the data like this:

power_line_gdf = power_line_gdf.to_crs(epsg=3036)

And then save it back out to file again, and it looks like this:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3036" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "LineString", "coordinates": [ [ -4448151.722213513217866, 15265232.02820273861289 ], [ -4436899.918848425149918, 15281201.161440940573812 ] ] } }
]
}

It is somewhat obvious from the data itself, but this puts the data way off in another location, way off the west coast of Africa (it should be in Mozambique). Am I missing something blindingly obvious here?


texts: How to create (online / svg) map with rivers highlighted and without borders?

texts: I want to create a world map on a website without administrative items such as borders, labels and without roads, only the continent outlines, but with big rivers visible (not like in Google Maps, where rivers aren't visible if I zoom out to view the whole world). Some kind of background like a relief map would also be great. What's the best way to do this? I thought about using OSM data – but with which framework / library / editor?


texts:
My goal: I'd like to publish a secured WMS on ArcGIS Server and then use MapProxy for limiting the access to it. For each user, the same WMS layers shall be available, but a different geographical extent (i.e. each user shall see the data of a different municipality / province).

What I know: By using the MapProxy Authorization API (limited_to), a geographical extent can be set.

What kind of user management do I have to implement for the authorization by user?


texts: Is there a fast way to copy large raster datasets between File Geodatabases?

texts:
Why does it take so long to copy a large raster dataset from one file geodatabase to another?  I can copy the entire database in windows explorer in minutes but copying a single raster in ArcCatalog takes forever.
I always assumed it is not advisable to use the windows explorer to copy data between to File Geodatabases but is there a safe way to do this?  If so, how?


texts: R spTransform() in lidR - error: non quantizable value

texts:
I am trying to reproject a las file in lidR (doawnload file here) from epsg:2232 to epsg:32613, and get the following error message after about 3 seconds of processing:
las <- readLAS("LD31201353.las")
projection(las) <- sp::CRS("EPSG:2232")
spTransform(las, CRSobj = "+init=epsg:32613")
#> Error in fast_quantization(X, scalex, offsetx) :
#>  Non quantizable value outside the range of representable values of type 'int'

Has anyone encountered this, and might you have suggestions for a fix (or what I'm doing wrong)?
The las file natively does not have its 'WKT OGC COORDINATE SYSTEM' field or  projargs slot populated, but I get the error with and without assigning them.
las@header
#> File signature:           LASF
#> File source ID:           0
#> Global encoding:
#>  - GPS Time Type: Standard GPS Time
#>  - Synthetic Return Numbers: no
#>  - Well Know Text: CRS is WKT
#>  - Aggregate Model: false
#> Project ID - GUID:        00000001-0012-07e3-0000-000000000000
#> Version:                  1.4
#> System identifier:        Merrick
#> Generating software:      MARS
#> File creation d/y:        188/2018
#> header size:              375
#> Offset to point data:     377
#> Num. var. length record:  0
#> Point data format:        6
#> Point data record length: 30
#> Num. of point records:    1965664
#> Num. of points by return: 1634171 285040 42901 3469 82 1 0 0 0 0 0 0 0 0 0
#> Scale factor X Y Z:       0.001 0.001 0.001
#> Offset X Y Z:             3122231 1354500 11639
#> min X Y Z:                3121464 1353000 11008.84
#> max X Y Z:                3123000 1356000 11782.48
#> Variable length records:
#>    Variable length record 1 of 1
#>        Description: WKT Information
#>        WKT OGC COORDINATE SYSTEM: PROJCRS["NAD83 / Colorado Central (ftUS)",
#>     BASEGEOGCRS["NAD83",
#>    [...] (truncated)


texts:
I have installed PostgreSQL 9.2 and PostGIS 2.0.3 from Stackbuilder. My template_postgis_20 does not have a geometry table only spatial_ref_sys table.

Is that a problem and, if not, then how can i create table using this postgis template with Geometry?

I have both tried Windows 7 and 8.


texts: Importing shapefile to SQL database, attributes are all NULL

texts:
I used ogr2ogr to import a shapefile into a MySQL database using the following:

ogr2ogr MYSQL:"db,host=localhost,user=root" -nln "map" -a_srs "EPSG:4683" map.shp -overwrite -addfields -fieldTypeToString All

but when I checked the table, all the attributes said "NULL" instead of their original values. What happened?


texts:
I have a set of vector lines (which represent translation vectors for georeferencing).
Based on the layer CRS, how can I compute the coordinates of each line in an attribute table ?

i.e (Xmax-Xmin; Ymax-Ymin) for each line.

I have managed to get the line length with $length in the expression editor but I am struggling with other geometrical parameters.

Edit: I am looking to "store the formula" within the layer. The idea is to keep these values automatically updated with the geometry. (Auto)update should happen if a new line is added or modification of the current geometry is done.


texts: Mapnik PGRaster Plugin creating 8x8 tiles with visible boundaries

texts:
I have used the latest addition to Mapnik PGRaster plugin (https://github.com/mapnik/mapnik/tree/2.3.x) by Strk(https://github.com/strk/mapnik/tree/2.3.x-pgraster) and merged recently.

I am able to get the images properly projected on the Map. But the images are tiled into an 8x8 manner and I could see the boundaries of the images at certain zoom levels.

I have added a Raster with a black patch so that the tiling boundaries would be more clear.

I am using OSM, Mapnik(2.3.x branch), ModTile and Renderd.

Also, when I run Renderd, I am getting logs like:

  renderd[15011]: Rendering projected coordinates 14 4896 6056 -> -8061966.247298|5205055.878110 -8042398.368057|5224623.757351 to a 8 x 8 tile

I have added a layer via layer-shapefile xml as follows:

<Layer name="MS-GIS" status="on" srs="+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs">
    <StyleName>raster</StyleName>
    <Datasource>
      <Parameter name="type">pgraster</Parameter>
      <Parameter name="table">ms_test</Parameter>
      <Parameter name="raster_field">rast</Parameter>
      <Parameter name="dbname">gis</Parameter>
      <Parameter name="estimate_extent">false</Parameter>
      <Parameter name="extent">-20037508,-19929239,20037508,19929239</Parameter>
    </Datasource>
</Layer>

How can I get the tiling to be smooth without boundaries?


texts: PostGIS. Obtain the rotation between two equal Geometries

texts:
I've got two features with the same drawing (poliline) but one is rotated relative to the other. I would like no know in PostGIS the rotation between the two geometries to apply this rotation to a third geometry (poliline).

I know how to do the last step with the st_rotate function but i don't know how to do the first one.

Regards,


texts: Lost geometry when joining spatial and non-spatial tables

texts:
I have two tables I am joining: 1 spatial (Table 1) and one not spatial (Table 2).

Table 1 looks like this:

| ID | geometry |

Table 2:

| ID | ... |

When I perform the query:

  CREATE TABLE AS SELECT * FROM t1, t2 WHERE t1.ID = t2.ID;

I get a table with the geometry column included. But unfortunately it is no longer recognized as a spatial table (spatialite_gui at least). Also, this mean I can't do what I really want to do, which is dump the table as shapefile.

So my question is: Why is the new table not recognised as a spatial table with a geometry column?


texts: GRASS GIS not quiet when called from Python although quiet option is set

texts:
I am working with GRASS GIS 7.4 and I call it via Python and GRASS session.

If I run this command:

r.cost(input='costraster', output='costsurface', start_points='start', stop_points='stop',
       outdir='movementdirections', memory=ram, overwrite=True, quiet=True)

I get loads of console output although I have set quiet=True. Why is that so? Also it is in red color...


texts: Problem with first page of data driven pages

texts:
I'm working on creating a set of data driven pages of quarter sections in a township. QutrSecs is the layer, QTR is the Name field (e.g. NE, NW, etc.) and sort field is SECTIONS (e.g. 01 -36), so that pages will ordr by section 1, then section 2, and so on.

This works fine for all pages except the first page, where it displays the entire section. Name on that page is NE. I can't figure out why this would be the case, and more importantly, how to fix it. This should be something simple, but can't find it.


texts: Renaming several fields using arcpy.AlterField_management()?

texts:
After several attempts and researches (I haven't found examples for renaming several fields directly without using a list), I tried the following code to rename several fields. However, I get the following:

  ERROR 000664: Invalid input: The type of dataset is not supported.

arcpy.AlterField_management(in_table = "file.shp",
                       field = ['Old_title','Old_title_1','Old_title_2'],
                       new_field_name = ['New_title','New_title_1','New_title_2'])


texts: Using FME to draw polygons around roads using aerial photography

texts:
I am trying to use FME to draw polygons around the extents of road using aerial imagery.

I have a Road Centreline table that I have buffered 5m and put this through a RasterGCPClipper so only the road and a meter or so off the road is shown on the aerial photography (depending on road width). I have the output going through the RasterResampler and then to the RastertoPolygonCoercer.

After this transformer I am stuck since I know the output will be thousands of little polygons the size of the designated cell size I have specified in the RasterResampler.

What I want to do is group on a certain RGB range or something similar. I know this may be a little crude but I just want to see if I can do something like this, and some of the roads have tree coverage that will muck it up too.

I want to output this to MapInfo and I am using an ECW file for the imagery.


texts: Problem with intersects operator in Postgres

texts:
I have this table below

CREATE TABLE public.geometry_polyline_volume (
    id          serial not null primary key,
    distance    float not null,
    height      float not null,
    coordinates path not null );

I am trying to make this query:

SELECT id
    FROM public.geometry_polyline_volume
    WHERE box'((15,47),(16,48))' ?# coordinates;

I am getting this error:

ERROR:  operator does not exist: box ?# path
HINT:  No operator matches the given name and argument type(s). You might need to add explicit type casts.

How can I find if a path is intersecting with a box?


texts: Making tiny objects visible for their location when zooming out

texts:
I'm adding small features to QGIS that I need to be in a precise location (trail signs) like so:
1:500

However, when I zoom out, they become invisible:
1:10000

If I switch the sizing to millimetres (rather than meters at scale) then I have to move the sign as I zoom out to keep it from overlapping the trails, then the position is all out of whack when zoomed in.
So, how do I create a precise placement, but still make tiny objects visible for their location when I zoom out?


texts: Problem with Georeferencer not scaling TIFcorrectly

texts:
I am trying to use QGIS and the Georeferencer to digitize paper maps but I am unable to get the raster to rotate afte I pick my control points.  Below is a brief outline of the steps I used to transform my raster.

Added two control points and picked corresponding points in map canvas.
The transformation type is set to Helmert
Resampling method: Cubic
Compression: Deflate
Specified Output file name
Checked "Load in Qgis...", then Ok
Then clicked run button.

The attached image is the results in QGIS.

The image should be much smaller


texts:
When I try to use field calculator in ArcGIS with polygon data in a geodatabase, I get an error "Empty value for ObjectID=1". I can manually change that field but when I try to use field calculator, I get an error. If I try to calculate the field with different sorting, I get the same error, e.g. "Empty value for ObjectID=252845".

What could be the problem?


texts: Loading all rasters in a directory and naming the variables iteratively

texts:
Hi I'm working with prism data and I already have clipped (masked) rasters (*.bil files) of a specific area for each day of a year. I'd like to open all these rasters with rasterio, convert them into numpy arrays, get the average value for each array, then an 8 day average of daily arrays for the entire year.

The first step to doing this is to open up all the bil files iteratively, an keep the names of the files as the variable names, but I cannot seem to get this to work.

import glob, os
import fiona
import rasterio
from rasterio.mask import mask
import numpy

path = '/home/demios/Desktop/Prism/Temps'

notebook = dict()

for filename in [each for each in os.listdir('/home/demios/Desktop/Prism/Temps') if each.endswith('.bil')]:
    with open(path+'/'+filename, "r+") as filehandle:
        notebook[filename] = rasterio.open(filehandle)

Near as I can tell this is what my code is supposed start off looking like. It does not seem to do anything though.

open with rasterio having matching variable names as files and append
variable names to list "rasterload" while opening
using last "rasterload convert loaded rasters into numpy arrays with
iterative date matching variable names and make a list "arraydates"
of them.
get average of each "arraydates" into a tuple named "dailyavg" via
np.mean(variable)
Average dailyavg every 8 days (ie 1-8, 8-16, 16-24, 24-32) into a new
tuple 8 day averages.

Anyone able to provide a path forward?


texts: I have defined an OpenLayers.Select.Control - I need only hoveron and hoverout events. I have implemented it using callbacks that and it works. But I need to disable just on feature clicking functionality. How can I do that?


texts: ArcGIS Vehicle Routing Problem Service via REST

texts:
I am trying to get a Vehicle Routing Problem to work as a geoprocessing service on one of our servers.

I used the default_date with unix time stamp in miliseconds and set the times of orders and depot in the format of "10:00 AM UTC". Now I get a time constraint violation for the orders ("ViolatedConstraints": 32)

If I take out the line "MaxViolationTime1":0 from the orders in my REST url, a route is calculated with correct travel times and distances and I can see the results. Especially the following lines (attributes of output route) caught my eye:

...
"TotalCost": 42.5174562651664,
"RegularTimeCost": 42.5174562651664,
...
"TotalTime": 42.5174562651664,
"TotalOrderServiceTime": 20,
...
"TotalTravelTime": 22.517456265166402,
...
"StartTime": 1369180800000,
"EndTime": 1369183351047,
...
"TotalViolationTime": 30.35018770210445,
...

I can change the default date and the earliest/latest start time, but the route starts still on the current date, at 2 AM. Could it be that some date and time options are "hard coded" somewhere in the module that was used to create the geoprocessing service? I checked the module again and can't find any options to change default date and time settings.


texts: Why is geometry being dropped when outputting to a MapInfo Writer in FME?

texts:
In some instances geometry is being unexpectedly dropped when writing to a MapInfo writer. This occurs when FME has a conflict with the Format Attribute of 'mapinfo_type' which has been attached to a feature from a non-mapinfo source such as SQL Server Spatial. The conflict arises when the mapinfo_type is point or line but the feature it has been passed to is a polygon. Sorry I don't have enough reputation to post an image of a simple workbench so hope the URL works or the following description:

Description: It will happen if we use a SpatialRelator transformer to perform a simple 'point INTERSECTS polygon' exercise. If the base polygon features come from a SQL Server Spatial reader but the point candidates originate from a native MapInfo reader the 'mapinfo_type' format attribute of the points is transferred to the base polygon features (as this attribute did not already exist on the base feature).  The result when writing out to a MapInfo writer is that FME sees a mapinfo_type of point on a feature which is actually a polygon and drops the geometry of any polygons that intersected.

We've worked around it by the following:

Write to SQL Server instead,
or if native Mapinfo output is required do any one of the following:
Expose mapinfo_type on any offending MI input datasets then remove it just prior to writing to output,or,
add a Geometry_Coercer transformer to the MI input dataset and force it to polygons.

Any comments on better solutions?


texts:
I'm trying to calculate pair distances among the points of my dataset. I have point in EPSG:3857.

This is the code I ran:

I import the data, then:

coordinates(datPb) <- c("X","Y")

proj4string(datPb) <- CRS("+init=epsg:3857")

datPb_owin <- as.owin(list(xrange=c(min(datPb$X),max(datPb$X)),yrange=c(min(datPb$Y),max(datPb$Y))))

datPb.ppp <- ppp(x=datPb$X, y=datPb$Y,marks=datPb$Pb, window=datPb_owin)

pairdist(datPb.ppp)

The pairdist function return a matrices of distance but I know (from field measurements and comparison with QGIS on a subsample of points) that these distances are not correct. I though it might be an error in the understanding of the EPSG but when I save back my data as a shapefile and measure distance in QGIS again they are fine.

here is a subset of my data:

quadrat X   Y

28  271322.8329 6231012.4650

29  271313.2814 6231006.6400

30  271319.1768 6231013.2850

31  271326.7498 6231017.9710


texts: Create Text files from a loop with Python?

texts:
I need to create some files to hyperlink to from Arcmap. I have a folder of jpeg files (100.jpg, 101.jpg etc...). I need to loop through them and create a text file of the same name except replace the *.txt at the end with *.fsv (100.fsv, 101.fsv etc...). This text file of the same name (with the txt replaced with fsv) is a set of instructions on how to open that jpeg in a piece of software so I need to write some instruction in it for example 100.fsv would contain:

open=100.jpg
set to 100

everything in the txt file is static except for the name of the jpg it opens. Thanks for any help.


texts:
I am attempting to write a script to include an input from the user that represents a "Minimum Confidence Threshold". I want to use this value to only save & print records to the feature class that have a confidence level above (greater than) this threshold.
The data is from a .txt file.

My current code is still printing all 411 records. If I have a confidence level of 99 for example, it should only print out 23 record. I think there must be an error with the way I have arranged my statements or indentations.Any advice?

import arcpy

work = raw_input("Enter the full path of WildlandFires.mdb: ") # C:\Scripts\Lab 7 Data\WildlandFires.mdb
arcpy.env.workspace = work # Set the workspace to the geodatabase
arcpy.env.overwriteOutput = True

iFile = raw_input("Enter the full path of wildfire text file: ")# input text file C:\Scripts\Lab 7 Data\NorthAmericaWildfires_2007275.txt

output = raw_input("Enter the name of the output feature class: ")
threshold = raw_input("Specify the minimum confidence threshold: ")

f = open(iFile, 'r')# input text file C:\Scripts\Lab 7 Data\NorthAmericaWildfires_2007275.txt in "read" mode

lstFires = f.readlines() # Read the lines of fire data from the input file
f.close() # close the file

field = ["SHAPE@", "CONFIVALUE"] # represents the field names for each row in the feature class

# Process: Create Feature Class
arcpy.CreateFeatureclass_management(work, output, "POINT")

arcpy.AddField_management(output, "CONFIVALUE", "FLOAT")

cursor_new = arcpy.da.InsertCursor(output, field)# add the new points that you will create to the "NewFires" feature class

cntr = 0 # counter variable

for fire in lstFires:
        if 'Latitude' in fire: # Skip the header
                continue
        pnt2 = arcpy.Point() # Create a new Point object
        lstValues = fire.split(',') # create a list of values "lstValues" based on the type of delimiter
        str_int = lstValues [-3:]
        str_int = list(map(float, str_int)) #turns string to integers
        latitude = str_int[0] # Latitude
        longitude = str_int[1] # Longitude
        confid = str_int[2] # Confidence Value

        if (confid > threshold):
                pnt2.X = longitude # Assign the X and Y properties of the Point object
                pnt2.Y = latitude

        row = [pnt2, confid] # Create a new row for the feature class
        cursor_new.insertRow(row) # Insert the new row to the feature class
        cntr = cntr + 1 # update the counter
        print "Record # " + str(cntr) + " written to feature class"

del cursor_new # release the cursor lock on the feature class


texts: Using a geodjango PointField with geography=True, my distance calculations are wrong

texts:
Here's my model:

class ExchangeArea(models.Model):
    """
    An Exchange Service Area, calculated from a phone exchange and surrounding census tracts.
    """
    coordinates = models.PointField(null=True, geography=True)
    geom = models.MultiPolygonField(null=True, geography=True)
    objects = models.GeoManager()
    adopts = AdoptionManager()

    def __unicode__(self):
        return '(%f, %f)' % (self.coordinates.coords[1], self.coordinates.coords[0])

And here's my ipython terminal output:

for e in ExchangeArea.objects.distance(ea.coordinates)[:10]: print e.distance, e.coordinates.wkt

0.0 m POINT (-121.4749999999999801 36.5399999999999991)
34.3149213171 m POINT (-87.4350000000000023 32.2049999999999983)
34.1483824654 m POINT (-87.6800000000000068 31.6400000000000006)
1.30461680198 m POINT (-122.4000000000000057 37.4600000000000009)
27.2054594705 m POINT (-94.2750000000000057 35.9950000000000045)
34.2874590776 m POINT (-87.7500000000000000 30.3550000000000004)
29.441773469 m POINT (-92.1550000000000011 33.8650000000000020)
5.51086451056 m POINT (-116.1352147188304826 35.1775313564249501)
2.11899150541 m POINT (-119.6700000000000017 37.6499999999999986)
27.6825400749 m POINT (-93.7950000000000017 36.1649999999999991)

Using this code as a reference, the first non-zero distance should be:

3151063.676 m

Update

So using south I've worked around this issue, ideally retaining the same level of accuracy (with reduced speed) by switching back to a normal geometry field (removing geography=true from my model) and using spheroid=True on my distance queries. They are now in a more sensible order of magnitude and still treat the earth as an spheroid (which is why I was using geography=True in the first place). Still would like to know what the deal with geography fields is. Leaving this open until someone can explain.


texts: GEE unsupervised classification (clustering) of Landsat 8 cloud free composite - User memory limit exceeded

texts:
Can anyone explain why I am getting the memory limit error when trying to classify using the wekaKMeans algorithm over a cloud-free composite of Bangladesh? I'm trying to differentiate mangrove pixels.

I'm using a rectangle polygon that is the same size as the one used in the GEE tutorial (this one: https://developers.google.com/earth-engine/clustering) to try to avoid this as I assumed it would input the same amount of data into the algorithm.

// Specify feature collection by country
var region = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw')
  .filterMetadata('Country', 'equals', 'Bangladesh');
var landsat8_collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');

// Filter to the scenes that intersect your study region.
var landsat8_studyArea = landsat8_collection.filterBounds(region);
var landsat8_SA_2015 = landsat8_studyArea.filterDate('2015-01-01', '2015-12-31');

 // Get an image collection and center the display.
Map.addLayer(landsat8_SA_2015, { min: 0.05, max: 0.8, bands: ['B6', 'B5', 'B4'] }, 'Landsat Collection');
Map.centerObject(region, 7);

// Cloud mask function
// Specify the cloud likelihood threshold.
var cloud_thresh = 40;
// Create the cloud masking function.
var maskClouds = function(image){
   // Add the cloud likelihood band to the image.
   var cloudScore = ee.Algorithms.Landsat.simpleCloudScore(image);
   // Isolate the cloud likelihood band.
   var cloudLikelihood = cloudScore.select('cloud');
   // Compute a mask in which pixels below the threshold are 1.
   var cloudPixels = cloudLikelihood.lt(cloud_thresh);
   // Mask these pixels from the input image.
   // Return the masked input image.
   return image.updateMask(cloudPixels);
};

// Mask the clouds from all images in the image collection with the map function.
var landsat8_SA_2015NoClouds = landsat8_SA_2015.map(maskClouds);
// Center map.
Map.centerObject(region, 7);

// Median reducer and display cloud free image composite
// Reduce the collection to the median value per pixel.
var median_L8_2015 = landsat8_SA_2015NoClouds.median();
// Print the information of the reduced image.
print(median_L8_2015, 'median_L8_2015');
// Display reduced image in the map window.
Map.addLayer(median_L8_2015,
 {min: 0.05, max: 0.8, bands: ['B6', 'B5', 'B4']},
 'median composite of cloud free images');

// Unsupervised clustering
// Make the training dataset.
var train = ee.Geometry.Rectangle(88.2, 22.9, 90.8, 21.2);
var training = median_L8_2015.sample({
  region: train,
  scale: 30,
  numPixels: 5000
});

// Instantiate the clusterer and train it.
var clusterer = ee.Clusterer.wekaKMeans(30).train(training);

// Cluster the input using the trained clusterer.
var result = median_L8_2015.cluster(clusterer);

// Display the clusters with random colors.
Map.addLayer(result.randomVisualizer(), {}, 'clusters');

Here is a link to my script: https://code.earthengine.google.com/0dc528e63b6f340baa82bef70ce63c92


texts: Accessing a OGC WMS server and publishing it to ArcGIS FlexViewer

texts:
I am trying to publish an OGC WMS service to my Flex Viewer. Within the configuration files amongst a few other layers I have the following:

        <layer label="EARSS National Outage Map" type="wms" visible="false" visiblelayers="earss:national_outage_map" skipgetcapabilities="false" alpha="0.6"
                autorefresh="60"
                url="https://earss.extranet.ornl.gov/geoserver/ows?service=wms&version=1.1.1&request=GetCapabilities/obs?srs=EPSG:4326">

The service does not appear within the TOC of the viewer. I tried the trick of creating "a service with a service" and that did not work.


texts: This is the scenario: I have thousands of shapefiles organized in different folders. My purpose is to merge all the shapefiles with the same name into a new single shapefile. I know I can import all of them into QGIS and merge them manually with MMQGIS plugin, but I was thinking if it would be possible to automate the operation with a script, as doing it manually would be extremely time consuming...


texts: Error when path contain spaces - Arcpy Integrate

texts:
I'm using a simple script to run Integrate tool in a default gdb and feature class, where user can indicate the path where this gdb is located.

This script works well when path indicated by user don't have any white spaces, but gives an error when path contain spaces.

I know that I need to use path in a raw string format (r' ').
But, how to do this since path is a parameter inputed by user?

import arcpy
import os
from os import path
from arcpy import env
from arcpy import da
from arcpy import management

# user input - GDB path
gdb_dir = arcpy.GetParameterAsText(0)

gdb = "Analysis.gdb"
fdsND = "ND"
network = "roads"
fullnetwork = os.path.join(gdb_dir,gdb,fdsND,network)

desc = arcpy.Describe(fullnetwork)
SR = desc.spatialReference

arcpy.Integrate_management (fullnetwork, "1.0 Meters")


texts: How to georeference point cloud data that was no spatial information?

texts:
I'm working on a project right now that the final goal is to compare two DTMs of a glacier surface from different years to visualize the change over time of the glacial surface. The first DTM was already created and given to me. It was created from data from opentopography (aerial photogrammetry mission flown in 2005). I am tasked with creating the second DTM and was given raw point cloud data in .las format.

This is were I things get tricky. The survey team that took these scans did not use any ground control and no GPS points were taken so the point cloud is floating in space. So far I have tried these methods:

Create a .tif in CloudCompare on the Z axis and then create hillshade
in ArcMap to try to georeference to a hillshade created from the 2005
data.This is proving difficult because there are no points that are
the same in both images because the surface is covered in snow so all
the results have a lot of assumptions.
Create DEM in ArcMap using LAStools. This produces a DEM at an angle,
because the point cloud was collected using terrestrial LiDAR. I
can't figure out how to project it on the Z axis.
Create a DEM by using the LAS to Multipoint tool in ArcMap, creating
a TIN then a DEM. Same thing as above, results at an angle.

If anyone has any suggestions on a route to take that would be immensely helpful. I am running out of ideas. Ultimately I will use both DTMs in Geomorphic Change Detector to visualize the change in the surface between the two years.


texts: Set an entire layer's opacity instead of feature-wise without using multiple map divs?

texts:
The issue here is that when you have a vector layer with polygons on top of each other and set layer opacity to something under 100%, their opacities seem to stack (which is a problem if you've used transparency to render an underlying layer visible). This also affects coloring.

So as the title says, do OpenLayers 3/4 have a native way of setting layer opacity post-feature-rendering? If not, is there a workaround for this that doesn't involve rendering multiple map divs on top of one another and passing events such as drag between them?


texts: Fill in gaps between rooms in floor plan using ArcGIS Desktop?

texts:
I am currently compiling a GIS of floors of a building. However, the current shapefiles have gaps between the rooms (the thickness of the walls-- shown below). I want to do some analysis on these rooms that require them to share boundaries.

How can I do this using ArcGIS Desktop?

I also use ArcPy and R.

EDIT

To follow up - yes I do need to know which rooms they share boundaries (Walls) with/are adjacent to. I also did a quick (rough) sketch of just a few of the areas that need to be sharing boundaries. Ideally- all of these rooms would be sharing walls at certain points with others. The dark green indicates areas that would most definitely be sharing a common boundary while the red indicates areas where there are furnaces, water heaters, etc. and don't need to be sharing common boundaries.


texts:
I use Google Earth Pro (7.1) on a daily basis to assess the acreage of US agricultural fields using the GE ruler.

What is the accuracy we can expect from Google Earth Pro for surface calculation? (the fields don't always have simple shapes and are sometimes in remote areas)

I found a few threads that partially cover this topic. But I didn't find any clear answer.


texts: GeoServer: points' symbols with a WKT polygon geometry?

texts:
I need to define custom points symbols in a SLD file for GeoServer. It seems that some WKT strings are accepted and others are not. I need to draw a closed line with a fill color.

The following example is accepted (symbols are displayed) but I can't specifiy a fill because it is just a linestring :

<WellKnownName>wkt://MULTILINESTRING((0 0, 0 1), (0 1, 1 1), (1 1, 1 0), (1 0, 0 0))</WellKnownName>

I have tried the followings to be able to fill the shape:

<WellKnownName>wkt://POLYGON(0 0, 0 1, 1 1, 1 0, 0 0)</WellKnownName>
<WellKnownName>wkt://MULTIPOLYGON(0 0, 0 1, 1 1, 1 0, 0 0)</WellKnownName>

The symbols are not displayed and I don't know how to debug it...


texts: Merging overlapping features using Geopandas?

texts:
I have a Shapefile with overlapping features.
Is there an easy way to merge them with geopandas?

I found a way with shapely, but I would like to do it with geopandas directly.


texts: Where is the "input point vector layer attribute" in the "Count Points in Polygon" in QGIS 2.18.9?

texts:
In older QGIS versions the Count Points in Polygon vector analysis tool had an "input vector layer attribute to aggregate" to perform a statistical analysis (sum, mean, min, max) on the counted points in a polygon.

However, in the newer QGIS (using version 2.18.6) this tool has disappeared. Only the counts in polygon is present in the vector analysis and you can't select the attributes of the point layer anymore to perform a statistical analysis. The tool now only counts the points in the polygon without any analysis.

Does anyone has a solution where to find the statistical attribute analysis?


texts: How to calculate area in raster data in ArcGIS 9.3

texts: I run a supervised classification on a Spot image and now want to calculate the area of different classes without converting to shapefile, because converting to polygon creates too much polygons which is very impossible to calculate area from them?


texts:
I'm brand new to ArcObjects. I've got it installed, and now I want to try using it in a .NET MVC application. In the controller, I think I need to initialize the license before I can do anything with ArcObjects since this is a standalone web application. Is there a code snippet out there to do that? I've been reading various threads and I'm having trouble putting it together. I'm using 10.5

I know at some point I have to call this:

ESRI.ArcGIS.RuntimeManager.BindLicense(ProductCode.Desktop)

but there's also some other code about initialization (AoInitialize ?) that I'm not sure where to put or what to do with.


texts: QGIS Plugin installed but reported as invalid

texts: I'm just discovering QGIS so please excuse my ignorance! I have QGIS 2.4.0 installed, and am trying to install Ordnance Survey Mastermap data. I have installed the OS2OGR translator to handle the decompression of the .GZ files, but the plugin manager reports the plugin is "broken - missing metadata file." I have checked the plugin folder, and the metadata file is there - what am I missing?


texts: Automate Exporting of Rasters Using Renderer in ArcGIS

texts:
It is possible to export rasters in ArcMap by right clicking the layer-> Data-> Export:

You can then click 'Use Renderer' check box to ensure the image you export uses the current symbology/rendering settings for that raster layer.

There does not appear to be any geoprocessing tools that match this 'Use Renderer' functionality.

I need to export over 200 images with these settings. I'd rather not do this manually. Is anyone aware of any way to automate this process?


texts: Drawing shapes like rectangle, polygons on fly using ArcPy in ArcGIS Pro?

texts:
I am trying to write a Python script that when imported into a toolbox can be used to draw user defined shapes on fly by clicking mouse points on the screen inside a GIS program. The motto of drawing should accomplish tasks like selecting the area of interest from a larger raster file etc.

How can I achieve this in Python using ArcPy with ArcGIS Pro?


texts: How do i rescale feature classes of different scales for intersect analysis?

texts: before I intersect 2 feature classes of different scales (e.g. 1:200,000 soil map and 1:25,000 land use map), should I rescale them like how u would resample raster datasets to the same resolution before overlaying them together? If yes, how do I do so in ArcGIS? (I have ArcGIS10.1).


texts: QGIS - merge overlapping polygons based on attributes while retaining attributes of one of them

texts:
This question seems simple, but I cannot figure it out. The problem is - I have a layer of polygons with animal species habitats. Different species can occupy overlapping or even THE SAME area, but polygons for a single species cannot overlap.

I have two types of errors in the layer that I want to remove:

partially overlapping polygons for the same species must be merged into one
100% overlapping features of the same species must be deleted, while 100% overlapping features of another species must be retained (no "delete duplicated geometries")

dissolve and v.dissolve don't produce what I want. I tried more elaborated solutions, like splitting layer by species and working on individual species. While I can delete duplicates I cannot find a way to merge only overlapping polygons without manual editing.


texts: Is it possible to run actions directly from the map canvas?

texts:
Is it possible to run actions directly from the map canvas in Qgis?
If yes how can I do that?


texts: ArcGIS online: Publishing map service, error code 400

texts:
I update a webmap on ArcGIS Online on a daily basis from an mxd in ArcGIS desktop, from where I update 8 of the 10 layers. I have used the "overwrite an existing service" with success until recently. Now I get the following message:

  Packaging succeeded, but publishing failed. ERROR: code:400. Unable to
  delete item. This Service item has a related Service item, Bad syntax
  in request. Please consult the GP results window for enhanced error
  description or contact your administrator for the server logs.

None of the layers have relates, or can these be hidden? What has happened since I started getting this error notification is that I have moved/copied these 10 layers and the mxd to a onedrive location, could this affect the update?

If i build up a new .mxd from scratch I get a notice saying that the caching size is 31GB and will be very slow to upload. I did not get this message earlier, even with the same layers. I have never touched the caching settings.

Is there perhaps an easier way of updating a webmap daily? Is there a difference between adding (altered) feature layers and overwriting the map service?

I am using ArcGIS desktop 10.3.1


texts:
Actually I figured out some problems and corrected it. I declared the overlay just outside the init function. So now I got the 'featureselected' event working properly. But 'featureunselected' event is still not working.

//Vehicle Overlay
      var  vmarkers = new OpenLayers.Layer.Vector("Vehicle Overlay",
                {
                    eventListeners: {
                        'featureselected': function(evt) {
                            var feature = evt.feature;
                            var popup = new OpenLayers.Popup.FramedCloud("popup",
                                  OpenLayers.LonLat.fromString(feature.geometry.toShortString()),
                                    null,
                                    "/*pop up content goes here*/",
                                    null,
                                    true
                                    );
                            feature.popup = popup;
                            map.addPopup(popup);
                        },
                        'featureunselected': function(evt) {
                            alert("unselected");
                            var feature = evt.feature;
                            map.removePopup(feature.popup);
                            feature.popup.destroy();
                            feature.popup = null;
                        }
                    }});
    function init()
       {
        map.addLayers([gmap_streets, gmap_hybrid, vmarkers]);
        selectControl = new OpenLayers.Control.SelectFeature(vmarkers, {
                click: true,
                autoActivate: true
            });
            map.addControl(selectControl);
         }


texts: Setting relative paths for KMZs in Arc2Earth dialog?

texts:
Using Arc2Earth, I'm converting a large vector dataset to a KML with regions.

The resulting kmz works great on my computer but the files paths are all actual paths and not relative, so they don't work on other computers or the server.

I would really prefer not to edit all of the resulting KML files, so I'm hoping somebody could help me figure out how to set a relative paths in the Arc2Earth dialogue box.


texts: ArcMap crashes with Create Enterprise Database tool

texts:
I am trying to create an enterprise database. I have ArcGIS for Server enterprise edition/workgroup and ArcGIS Desktop 10.4.1. Also I installed SQLEXPRESS 2012 SP1 included from Esri in the ArcMap 10.4.1 installation media and I am connected to the instance in ArcCatalog using local user authentication. Everything is on the same machine

I am trying to use the Create Enterprise Database tool and I have followed all the steps in the documentation but before the tool is completed the dialog box says Geodatabase Created Successfully and ArcMap crashes.

When I look in the DATA folder the database has been created but it is not functional.

ArcMap also crashes if I try to connect to the db that is created with the Create Enterprise Database tool.

And then arcMap crashes


texts:
I am looking for a way to simulate a gravity model using a point-based layer.

All my points are assigned a z-value and the higher this value, the larger is their "sphere of influence". This influence is inversely proportional to the distance to the center.

It is a typical Huff model, each point is a local maximum and valleys between them indicate the limits of the zone of influence between them.

I tried several algorithms from Arcgis (IDW, cost allocation, polynomial interpolation) and QGIS (heatmap plugin), but I found nothing that could help me. I also found this thread, but it isn't very helpful for me.

As an alternative, I could also be satisfied by an way to generate Voronoi diagrams if there is a way to influence the size of each cell by the z-value of the corresponding point.


texts: Adding layer from cliprasterbymask using PyQGIS

texts:
In my plugin I would like to use the gdal tool "cliprasterbymask" and add the new clipped raster to my project.

According to the Documentation the output has the type "raster". So it returns a raster object.

Thus, I have tried the following:

 clip = processing.run('gdal:cliprasterbymasklayer',{'INPUT': wasserstand,'MASK': gebeude,
                   'NODATA': -1, 'ALPHA_BAND': False, 'CROP_TO_CUTLINE': True, 'KEEP_RESOLUTION': True,
                   'OPTIONS': "", 'DATA_TYPE': 5, 'OUTPUT': outFile})

 self.iface.addRasterLayer(clip)

Unfortunately, I get the following error message which I don't know what it means.

  TypeError: QgisInterface.addRasterLayer(): arguments did not match any overloaded call:   overload 1: argument 1 has unexpected type 'dict'   overload 2: argument 1 has unexpected type 'dict'

How do I get the outputraster from cliprasterbymasklayer and what is actually output?

Apparently, its not just a simple raster layer.


texts: Draw cirlcle on Equirectangular projection

texts:
I'm trying to draw a circle (actually a great circle) on an equirectangular map to archive something like this http://i.stack.imgur.com/7EliU.png (from here).

So I have the lat/lon coordinates of the center of the circle on the surface of the sphere and I already have the coordinates (the position of the sun) and would like to have a function describing the border of this circle on an equirectangular map by using these coordinates.

I think the answer mentioned above fits my needs, but sadly I don't understand it entirely. Is there an easy way to do this?


texts: Getting Properties of Next Record Using Calculate Field Tool in ModelBuilder?

texts:
I'm trying to calculate the distance between the centroids of adjacent records (adjacent in the attribute table, not necessarily in space) into an attribute field. I'm planning on using the truecentroid.x and truecentroid.y properties of one record and the next record in the attribute table with the Pythagorean theorem to calculate distance. but can't figure a way to get to the properties of the next record. I'm not a great programmer but I can usually hack something together.

Can someone suggest a starting point for me?


texts: lines of a SpatialLinesDataFrame corresponding to points

texts:
I have a spatial lines data frame which can be projected on to a map. The projection is defined as WGS84. I have a Spatial Points Data Frame which is also in the same location.

How can I isolate the streets/lines in the line data frame which contain, or are close to, a point in the points dataframe.

I am not able to see the geometry of the lines described in the Spatial Lines Dataframe.


texts: catchment area using SAGA and Arcgis 10.2

texts:
I would like to calculate catchment area. I used D8 algorithm and running in both arcgis 10.2 and SAGA GIS 2.1.2. However, the results from two software are different. The lowest value of catchment area in arcgis = 0; high value = 109325, while SAGA GIS lowest value= 625, highest value=47820. I think the result should be similar and maybe I do something wrong. Could you please explain what I did wrong? These are what i am doing,

the input: DEM25
In arcgis, to calculate catchment area: Dem25, then run tool flow direction, then run tool flow accumulation, then 25 (cell size) * 25 * flow accumulation

In SAGA GIS, terrain analysis -- hydrology --catchment area (parallel), elevation = dem25, sink routes=no set, weight = not set, catchment area = create, Step =1, linear flow = no tick, linear flow threshold =500. convergence = 1.100000001

Thanks so much for your help.


texts:
I am looking to run a query in OGR to write back all instances where there is only one record of an attribute. I was thinking that I would use the following:

ogrinfo -dialect sqlite -sql "SELECT COUNT(*) FROM roads AS roadcnt WHERE roadcnt = 1 GROUP BY road_name" roads.shp >report.txt

This query returns an error: near "road_name": syntax error


texts: geology problems: three-dimensional fault line - r.plane (GRASS) from polyline

texts: I should use the r.plane script (GRASS) but imputing a polyline (trace of my fault line) and not a point; is this possible? I could use batch process, but I should be charged a lot of points and then fonder raster.


texts:
I am trying out the ArcEditor 10.0 trial version and am having problems with the Parcel Editor.  I am attempting to "select" using SELECT PARCEL FEATURE tab, nothing gets selected.  I did the START EDITING under Parcel Editor, clicked SELECT PARCEL FEATURE, tried selecting the parcel features to no avail.  After I click the SEL PAR FEATURE, it does show that features symbol on the pointer/scroll arrow.

When I click START EDITING, only the SEL PAR FEATURE is clickable, the rest of the other PARCEL EDITOR buttons are unclickable.

I am using shp files I created while using 9.3 ArcView.  Would that have anything to do with it (which I can't see any reason why it should)?

Some guidance in resolving this problem is appreciated.


texts: Using ArcGIS Online for public participation?

texts:
We have an ArcGIS Online Organizational account setup with several maps facing the public.

I would like to create an application using one of our maps that has the added functionality of allowing the public (anyone, not just those with an ArcGIS Online account) to place a point/marker on the map and submit it with a comment.

The comments would need to be posted to a database (or emailed to me) with coordinates.

Has anyone developed such an application, or seen any templates available similar to this?

ArcGIS Online's Map Notes are close to what I'm looking for, but editing map notes is restricted to the author.


texts: How do I display hatching labels in a different route measure unit than my coordinate system?

texts: I have a route that crosses several states from west to east which is in the  USA_Contiguous_Lambert_Conformal_Conic coordinate system. The linear unit of this coordinate system is in meters. I would like to keep this coordinate system because it has good distance accuracy but I need to display my hatching in miles. I have set it up so that the hatches are displayed every mile (hatch interval of 1609.34 m). When I label these hatches my labels are of course in meters. What is the easiest way to display my labels in miles? Is there an easier way than building a text expression?


texts:
I need to download tiles of an entire country to use offline.
I don't need a server o whatever, I just need to download tiles organised in folders (x y and z).
Should I use mapnick? Im on OSX but I can also use Windows 8.1
Which tool permit to organise tiles in such form?


texts: How to Legally Use Google Imagery and Maps in ArcGIS Desktop and Web Applications?

texts:
With companies like Arc2Earth and Spatial Energy offering solutions to view Google imagery and maps inside ArcGIS for Desktop, there is clearly a legal way to develop a similar product yourself. My question is how can this be accomplished and what licensing is required?

We would like to (legally) view and print Google imagery and maps inside ArcGIS for Desktop, Flex based web mapping applications, and JavaScript web mapping applications. My dream is to somehow provide this as a map service or image service using ArcGIS for Server.


texts: Making gridcell take on value of largest region within gridcell with QGIS

texts:
I am new to QGIS and very excited about all its features. We want to match gridcells to a region, when a gridcell lies in more than one region we want to assign the largest region to the gridcell.

Up till now we have been working with "connect with attributes on basis of place" to assign the region name to the gridcell. We are not familiar with R or Python.

Is there anyway we can do this without using these programs?