texts: Identify circle segments in polyline feature class

texts:
I am using ArcGIS and I try to find all roundabouts in the OSM data of Germany. Although, there is a OSM road type depicted as 'roundabout', the type is not used for the OSM data. At least not for the data I downloaded...

How can I find all those line segments that belong to a roundabout?

Do I have to check all features for their geometry? I tried already to see if the geometry object returns a 'curve' element when converted to a JSON object, but it only contains the coordinates of the vertices and the spatial reference.


texts:
I am working with copies of feature Datasets that come from a versioned SDE.  The copies are in an un-versioned 'development' SDE.  I am attempting to do a field calculation on a single feature class in the development sde.  I am trying to activate an editing session in python to work around the versioning but the arcpy.da.Editor is not activating an editing session.

I am more or less copying the code from these two sources:

http://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-data-access/editor.htm
https://geonet.esri.com/thread/77534

For some reason it will not work, and returns the "Objects in this class cannot be updated outside of an edit session" error.

What am I doing wrong?

import arcpy

workspace = "C:\\Users\\syoklic1\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\devsde.sde"
fc = "devsde.SYOKLIC1.main"
arcpy.env.workspace = workspace
arcpy.MakeFeatureLayer_management(fc, "fc_layer")
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()
with arcpy.da.Editor(workspace) as fieldedit:
    arcpy.CalculateField_management(in_table= "fc_layer",
                          field="MiantenanceArea",
                          expression="1234",
                          expression_type="VB",
                          code_block="")

edit.stopOperation()
edit.stopEditing(True)

Additionally, it works if I do it inside of ArcMap.


texts: Calculating percentage area within each cell using ArcGIS for Desktop?

texts:
I would like to quantify the percentage area covered by different categorical values (A, B, C, D, E, F, G) within each cell. The cells are defined by layer_1 while the categorical values are in layer_2.

Is there a specific tool or tool combination that I can use in ArcGIS to achieve this?

Hopefully the image gives a good example of what I mean.

I am looking at the spatial join tool + pivot tool to see if this can be done.

I am running ArcGIS 10.3.1 for Desktop (Advanced) with extensions for Spatial Analyst, XTools Pro and Geostatistical Analyst.


texts:
I made a Query in Mapbasic like :

Dim query As String
query = "Select * from some_table"

How to run query to form table? Is it possible to save selected value in variable without forming table?


texts: PostGIS ST_AsGeoJson makes geometry invalid

texts:
For the following SQL, I create a GeoJSON from a valid geometry, which makes it invalid:

select st_asgeojson(geometry, 6), st_isvalid(geometry), st_isvalid(st_geomfromgeojson(st_asgeojson(geometry, 6)))
from bd_region
where id = 310100;

the output:

{some GeoJSON}, true, false

the output JSON is here:
https://jsonblob.com/c62204c1-03f1-11e8-8467-f5eea275e9a5

After a search I found that may due to precision reason. How can I fix it?

  [
    [
      121.412516,
      31.244891
    ],
    [
      121.41456,
      31.243229
    ],
    [
      121.41456,
      31.243229
    ],
    [
      121.412516,
      31.244891
    ]
  ],


texts: Is a pdf geospatial if the Adobe Reader Geospatial Location Tool doesn't write anything?

texts:
Does a PDF contain geographical coordinates if the Adobe Reader Geospatial Location Tool doesn't write anything?

I mean that I am moving the mouse over the pdf and no latitude and longitude values will show.

Does this mean there are no coordinates stored in the map?


texts: MATLAB mapping toolbox: write multiple feature mapstruct

texts:
I have a few hundred road segments, described as lines, which I am trying to write to a single shape file using MATLAB's mapping toolbox. However I can only manage to get one feature per shapefile, which isn't very useful. I have developed a work around which calls GDAL from within MATLAB using system(), but the company I work for can be sensitive about us using software that they don't specifically support, so if I can find an all MATLAB option that would be ideal.

I assume that either it's not possible, or more likely that it's really simple and I'm missing something obvious. The shapewrite() function requires a 'mapstruct' structure item, but the documentation on how to write such a thing is not very clear, I've only succeded in writing single featured mapstructs. The shapewrite documentation only shows how to write a multi-feature shapefile from a mapstruct that was just read using shaperead(), but not how to write your own multi-feature mapstruct structure.

Here's what I have so far...

% Lines is a structure containing all of the individual features
LineIDs = fieldnames(Lines);
First = 1;
for LineIndex = 1:numel(LineIDs)
    LineID = LineIDs{LineIndex};
    Line = Lines.(LineID);
    % So now Line is a structure for that individual line, it contains
    % such fields as the X and Y vertices, the BoundingBox, the Geometry
    % type (line in this case), and other attributes. Line is my single
    % feature 'mapstruct'. It can be written to a shapefile using shapewrite.
    LineFileName = sprintf('TempFiles\%s.shp', LineID);
    shapewrite(Line, LineFileName)
    % Now use GDAL to merge this file with all the others into one large file.
    % These next lines are the ones I would like to replace.
    if First
        First = 0;
        MergedFileName = 'MergedFile'
        MergedFileNameF = 'MergedFile.shp'
        % If the file already exists, remove it.
        if exist(MergedFileNameF, 'file') == 2
            delete(MergedFileNameF)
        end
        % Write the file to be merged to.
        [status, ~] = system(sprintf('ogr2ogr %s %s', MergedFileNameF, LineFileName));
    else
        % Merge with the file.
        [status, ~] = system(sprintf('ogr2ogr -update -append %s %s -nln %s', ...
                      MergedFileNameF, TempFile, MergedFileName));
    end
    if status == 0
        % Merge was succesful
        % Delete temp files, etc
    end
end

My single feature mapstruct structures look like this:

>> Line
ans =

          Class: 'A Road'
             ID: 28841
       Geometry: 'Line'
           Name: 'Willowbank Road'
              X: [3.9323e+05 3.9335e+05 NaN]
              Y: [8.0548e+05 8.0547e+05 NaN]
    BoundingBox: [2x2 double]
      DftNumber: 'A93'
       FEATCODE: 25729
             CP: 80026

And they are stored within Lines, which looks like this...

>> Lines
ans =

    A: [1x1 struct]
    B: [1x1 struct]
    C: [1x1 struct]
    etc etc


texts: Avoid time out error Nominatim Geopy OpenStreetMap

texts:
How to avoid the GeocoderTimedOut('Service timed out') error in GeoPy using the OSM Nominatim geocoder? I am not interested in using an exception and skip some entries. I already use a time delay of 1.1 seconds in my script using time.sleep(1.1). Total database is 10,000 entries.
Nominatim Usage Policy


texts: OpenStreetMap converting to country -> administrative unit -> city database

texts:
I'm stuck on my task and I don't know what to do. The job is:

extract data from OpenStreetMap
create database model containing following schema and data

Schema looks like this:
+---------+     +------------+      +----------+
| country |     | adm_unit   |      | city     |
+---------+     +------------+      +----------+
| id      |--\  | id         |--\   | id       |
| name    |   \-| country_id |   \--| adm_id   |
+---------+     | name       |      | name     |
                +------------+      | coord_x  |  //coordinates - longitude and latitude
                                    | coord_y  |
                                    +----------+

So far I managed to:

download map itself (only for Poland to test it)
import it to database using osm2pgsql tool
perform some selects and realize, I don't know how to do it.

The main problem is: I can't see any relations between different items of map. Taking Poland as example: I can find Warsaw (capitol city), as well as "masovian vovoidenship" performing few selects, but there is no relation between these two. There are lot of empty columns in my db instead. Also none of column names suggests that it should hold any relation to "parent".
I don't know how to find those relations.
Is it something with imports that i'm missing?
Should I change anything in my default.style (I haven't changed anything inside it)?
Or maybe there's some kind of other approach to solve my problem?
Maybe there are free databases that contain data I need and are less painful to work with?
Geonames is no-go, their database is full of errors.
I'm really green when it comes to gis, I'm PHP programmer, that's completely different world for me. Unfortunately, even if I want to explore it, my deadlines won't allow that.


texts: How to create a raster from a CAD scan for ARCMap?

texts:
I have scans of CAD streets and such and want to get a raster for georeferencing into ARCMap.
What are the easy steps to make a targa, geotiff, etc. of a known size, say 2 miles along one side as a test size scenario or any other raster format that can be imported into a layer for subsequent geo reference? Not wedded to any particular form, size, technique, format; just best practice and practical steps. Thanks for any help.


texts: Geotransform data missing in L2A products on S3

texts:
I am trying to use Sentinel-2 L2A tiles from the sentinel-s2-l2a bucket. Unfortunately I noticed that the geotransform information (SRS, origin, scale, and skew) are missing on the tile files.

An example image is the following:

tiles/32/T/LS/2017/6/22/0/R10m/B03.jp2

Why do the L2A tiles not have the Geotransform specified?


texts: Spatial join in GeoPandas: Line within multiple polygons

texts:
I'm having a GeoDataFrame of lines and a GeoDataFrame of polygons. The polygons have an attribute with the altitude of that polygon.

For each line, I want to find in what polygon it is located. As output I would like to have the GeoDataFrame of lines, with the altitude attribute of the accessory polygon added to each row (each row is a line).

Below code gives the desired output for most lines.

output = gpd.sjoin(lines, polygons, how='left', op='within')

I say this works for most lines, as it does not work for lines which are located in multiple polygons, shown blue in attached image. For these lines, it gives a NaN value for the attributes. However, I would like it to find the attributes of both polygons, and eventually I keep the lowest altitude of the matched polygons.

Hopefully I made myself clear. Could anyone explain how I can do this?

Below I explain what I already tried, but this is a bit messy and maybe not the right direction:

As solution, I tried to iterate over each row of the polygon GeoDataFrame (from lowest altitude to highest), and then do the spatial join, only with lines which are not yet matched. But this didn't work for me, as the spatial join didn't work on single rows of the GeoDataFrame (crs issue).

polygons= polygons.sort_values(by=['altitude'])
output = None

for index, row in polygons.iterrows():
    if output is not None:
        toappend = gpd.sjoin(newlines, row, how='left', op='within')
        newlines = toappend[toappend['altitude'].isna()]
        output = output.append(toappend.dropna())
    else:
        output = gpd.sjoin(lines, row, how='left', op='within')
        newlines = output[output['altitude'].isna()]

AttributeError: 'Series' object has no attribute 'crs'


texts: Offsetting series of points perpendicular to line with another line as reference distance using ArcGIS Desktop?

texts:
I need to find out how to move/create a series of points which are already on a line perpendicular left or right to intersect with a second line which runs parallel to the first.

The distance between each of the output points does not need to be the same as the original, but the up/down relationship needs to stay in tact.

An example of this is that if I had a grid line running exactly north/south with points evenly spaced along it, I would need the new output points to be exactly east or west of the original line but sitting right on the second line

This may sound pretty complicated the way I am explaining it but maybe the picture attached may make it clearer.

When offsetting I need to be able to specify a maximum offset distance and would like the new points feature class to contain the same information as the original.

I would also like to be able to do it on a mass scale

I was thinking along the lines of:

Input Feature Point Class
Input Feature Line Class: Bearing reference lines
Input Feature Line Class 2: Intersect lines for new points
Output Feature Point Class: with original attributes

I am using ArcGIS Desktop 10.1.

Example


texts: Leaflet circleMarker's not working - but marker works fine

texts:
Using Leaflet, I'm adding around 15,000 circleMarker's to a geoJson layer. The circleMarker's, however, are not showing on my map when the layer is added, despite specifying the relevant styling options. Weirdly, the markers work correctly when I use .marker and a div icon, rather than .circleMarker. The layer is definitely being added to the map, and the circleMarker's are definitely being added to the layer.

What could be causing this issue? Heres my code:

function launchInteractiveMap ( attacksNodemap ) {

    const mapboxAccessToken = 'randomToken';

    map = L.map('interactive-map', { attributionControl: false } ).setView([51.505, -0.09], 2);

    var myRenderer = L.canvas({ padding: 0.5 });

    var darkStyle = L.tileLayer('https://api.mapbox.com/styles/v1/lbamber/ck13kjh340v741cmwjsgjqop9/tiles/{z}/{x}/{y}?access_token=' + mapboxAccessToken, {tileSize: 512,zoomOffset: -1}).addTo(map);

    // circleMarkers TRYING TO BE ADDED

    var nodeLayer = L.geoJson();

    attacksNodemap.forEach(function (item, index) {

        L.circleMarker([item.latitude, item.longitude], {
            renderer: myRenderer,
            color: '#FF0000',
            radius: 10
            }).addTo(nodeLayer);

    });

    // ADDING LAYER TO MAP
    map.addLayer(nodeLayer);
    map.flyTo([39.097845, -99.050912], 3);

}

launchInteractiveMap( attacksData );

I found this question which indicates that a crossover of d3 css could be to blame, however after removing all the links for d3 that I was using for other charts, the issue still persists.

I am also using other geojson layers with my map, however after removing these the issue was not resolved.


texts:
I would like to have a borderline around a subset South-South(SS) part of the map (Map 1) (More like highlighting the subset part).

I tried to clip a vector of the subset (map 2) to the entire map, but that did not work.


texts:
I would like to achieve the same output as for GDAL using the GeoTools library:

    CoordinateReferenceSystem crs = CRS.decode("EPSG:32632");
    String wkt = crs.toWKT();
    System.out.println(wkt);

GeoTools output:

PROJCS["WGS 84 / UTM zone 32N",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]],
  PARAMETER["central_meridian", 9.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["scale_factor", 0.9996],
  PARAMETER["false_easting", 500000.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","32632"]]

GDAL output:

PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32632"]]

There are some minor differences, including missing authority in some cases, indents, etc. For example,

GDAL:

    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],

GeoTools:

  UNIT["m", 1.0],

Is is possible to tune GeoTools such that it outputs the same WKT as GDAL?

What are the other ways to print the same WKT output to GDAL but not using native GDAL libraries?


texts: How I can create a raster and have the same overlay as another raster in QGIS 2.4?

texts:
I have made a graphical model in QGIS 2.4. I have to create a raster equation to rasterize a polygon shapefile. I want to rasterize my shapefile and I want to ensure that this new raster will perfectly overlay the position of another raster. I have used the tool "Rasterize" but it does not overlap my other image.

In ArcGIS we have some parameters that allow to specify to snap the pixel raster to another but I can't find the equivalent in QGIS. Can somebody help me?


texts: Proj/PyProj Convert between ocentric and ographic latitudes

texts:
I am trying to use Pyproj to perform the conversion between ocentric and ographic latitudes. I can do this using the cct command line tool:
echo 0 45 | cct -z0 +proj=geoc +a=3396190 +b=3376200
  0.0000000000   44.6617680466        0.0000           inf

based on the docs here.
When attempting to do this with Pyproj, I would expect to be able to define the projects and then apply the transformation. Something akin to:
from pyproj import Proj, transform
ocentric = Proj(proj='geoc', a=3396190, b=3376200)
ographic = Proj(proj='latlon', a=3396190, b=3376200)

ographic(0,45)

I can get the ographic projection to instantiate without an issue, but the ocentric projection fails with an invalid CRS error: CRSError: Input is not a CRS: +proj=geoc +a=3396190 +b=3376200  +type=crs
Has anyone been able to define a geocentric projection on a custom ellipsoid? If so, how?
Note: I also tried:
from pyproj import Proj, transform
ocentric = Proj(proj='latlon', geoc=1, a=3396190, b=3376200)
ographic = Proj(proj='latlon', a=3396190, b=3376200)

transform(ocentric, ographic, 0, 45)

and I get an error that the input is not a transformation. I believe that geoc is being silently ignored in this example.
Update
I was able to define a pipeline to make the conversion work as follows. Thanks @snowman2.
trans = pyproj.transformer.Transformer.from_pipeline('+proj=pipeline +a=3396190 +b=3376200 +step +proj=geoc')
lon, lat = trans.transform(0, np.radians(45), errcheck=True)
np.degrees(lat) = #44.6617680466192

Note the weird np.radians call. When I run this transformation using 45 (as in degrees, which should be the default because radians=False is the transform default) I see: ProjError: transform error: latitude or longitude exceeded limits.


texts: Using linear referencing in ArcGIS to plot line data

texts:
I have just started using the linear referencing tools in ArcGIS 10.2 but I'm having a bit of trouble. I'm trying to plot schemes as line data on to our road network layer which is also line data using start and end points (almost like routes).

This is what I have so far :

I've successfully created a route layer from our road network line data. This layer is split into around a 1000 sections ranging from 50m to 2500m with each having there own SECTION_ID.
I have a table which has numerous schemes in each with the following attributes : SCHEME_NAME, START_SECTION_ID, START_CHAINAGE (How far into that section in meters), END_SECTION_ID, END_CHAINAGE_ID (How far into that section in meters).

I'm using the Make Route Event Layer tool to create the schemes on the road network layer using the SECTION_ID as the route identifier field. The problem is I can only seem to plot the schemes which start and end in the same SECTION_ID as there is no option to provide an additional end route identifier field.

I'm not sure I'm going about this the right way so can anyone give me any advice?


texts: Intersect multiple binary rasters based on certain criteria in ArcMap

texts: I have multiple binary raster image (0.1). For example if I have n images, how can I extract the pixels that has the value of 1 and frequently occurred " the same pixels" in 70% of the n images (for example if I have 10 images, I need to extract the pixels that has the value of 1 in and frequently exist in at least 7 images.


texts:
I was wondering if there is a way to export in DXF from QGIS (version 3.4.4) this type of line (The one with the "E" as a pattern).

Here is how the line was created

Currently, when I export this type of line, the symbols are not kept and the result is two superposed polylines: one which is the green polyline and the other which is a black polyline with the same thickness as the size of my symbols "E".


texts:
I've been working with PostGIS for a while now, and I'm able to run queries over massive datasets reasonably easy and fast. However, I do have an issue now that I don't know how to solve. I did find some workarounds, although they mean lower performance, higher execution times, and more complex queries.

The objective is straightforward: to get all the points within a volume.
But, the volume is not a simple extrusion of a polygon, it is instead of the difference of a higher volume minus some smaller volumes here and there.

My first attempt: to tackle the problem was inspired by this question, although it didn't work:

Create the master polygon, using a third axis (Z) as the base for the volume:

ST_Polygon(ST_GeomFromText('LINESTRING(Longs, Lats, Zs...)'), 4326);

Extrude the polygon using the height of the volume, getting a POLYHEDRALSURFACE as a resulting geometry:

ST_Extrude(Previous_polygon, x=0, y=0, z=height);

Having a table with the information regarding each master volume and its small volumes to be removed from. I created a new geometry removing all these small volumes from the master volume:

ST_3DDifference(Master volume, combination of smaller volumes to be removed);

Query the table with the positions and joining the resulting master volume on

3dpoints.geometry &&& Master_Volume.geometry.

As I said, it didn't work because the points were also found outside the master volume. Not sure why! I guess I'll try to come back to that option should the following solution be too much of a hassle.

The second attempt I tried was:

Query the table with the positions straightaway, merging the volume as a 2D polygon, and including the filter on the Z-axis, something like:

ST_Contains(Master_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Master_Polygon.MinZ AND Master_Polygon.MaxZ

Of course, this included all the point within the smaller volumes I didn't want to consider, so then I added the condition to remove those points:

ST_Contains(Master_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Master_Polygon.MinZ AND Master_Polygon.MaxZ AND NOT (ST_Contains(Smaller_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Smaller_Polygon.MinZ AND Smaller_Polygon.MaxZ)

And it does work, although only when the smaller volume (polygon) is only one. Should the lower volumes be more than one, then the condition is not enough as it would create a different set of points based on each condition.  This gives all the points within the master volume but not in the lower volume A would be returned or all the points within the master volume but not in the lower volume B.

Question:
Is there a way to get all the points in the master volume by removing all those points with a condition which are in any of the smaller volumes?

I apologize for the long post, but its due to the complicated explanation of the issue.
I did find a solution, and it's by splitting the master volume into smaller volumes and merging it after the query. But it is not quite right, especially when I will have to update the geometries and make the master volume to change accordingly automatically.


texts: Count number of points within multipolygon shapefile using OGR

texts:
I have two shapefiles: 1) consists of points (hundreds of entries) and 2) multi polygons (tens of thousands of entries).

I simply need to count the total number of points contained within the polygon set (presumably using ogrinfo) and write out a shapefile containing these points.

Must use OGR/GDAL to do this. I feel that there is an ogr2ogr sql statement that allows me to select the points without performing a full on clip, which is expensive in terms of computation time (clipping the points by the polygon extent takes hours on hours and I cannot tell if the operation works). Or perhaps the best way would be a spatial join?

There is a similar question here:
https://stackoverflow.com/questions/29769991/python-ogr-select-and-write-points-which-intersect-with-polygon-layer/33399208#33399208


texts: ArcGIS 10 ToolValidator - Empty feature class and feature class within declared feature dataset

texts:
I'm new to ToolValidator for a script tool and I'd like to make two checks which then adds an error if either rule is violated.
I need to make sure parameter 2 and 3 (feature classes) are within parameter 1 (feature dataset).  The tool will fail (roughly 3/4 of the way through) if the feature classes are not contained in the feature dataset.

Also, the tool needs to display an error is the feature class is empty.

I know Tool Validator is the way to go.   I don't have much and am willing to read stuff.  Also, this doesn't work, but I've managed to piece it together with some other posts.

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    import arcpy
    self.params = arcpy.GetParameterInfo()

def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    if self.params[3].value==True:
        self.params[3].enabled==True
        fc=self.params[3]
        rows = arcpy.SearchCursor(fc)
    if rows.next():
        self.params[3].enabled==True
    else:
        self.params[3].enabled==False
    return

def updateMessages(self):
"""Modify the messages created by internal validation for each tool
parameter.  This method is called after internal validation."""
return
    if self.params[3].hasBeenValidated:
        self.params[3].clearMessage()
    else:
        self.params[3].setErrorMessage('Your feature class is empty. 300 Tool will not fail.')


texts:
I am trying to install pgrouting on ubuntu 14.04 like this:

Add pgRouting launchpad repository

sudo apt-add-repository -y ppa:ubuntugis/ppa
sudo apt-add-repository -y ppa:georepublic/pgrouting
sudo apt-get update

Install pgRouting package (for Ubuntu 14.04)

sudo apt-get install postgresql-9.3-pgrouting

But  "apt-get update" fails with:

W: Failed to fetch http://ppa.launchpad.net/ubuntugis/ppa/ubuntu/dists/trusty/main/binary-amd64/Packages  404  Not Found

W: Failed to fetch http://ppa.launchpad.net/ubuntugis/ppa/ubuntu/dists/trusty/main/binary-i386/Packages  404  Not Found

E: Some index files failed to download. They have been ignored, or old ones used instead.

How should I fix this?


texts: Correct way to determine the type of Leaflet layer?

texts:
What is the correct (or recommended) approach to identifying a layer "type" in Leaflet/Javascript?

I am able to successfully use layer instanceof, but that only works when you know the type at design time(which I expect is usually the case).

map.eachLayer(function(layer) {
    console.log('_leaflet_id=' + layer._leaflet_id + ' is layer type= '+ getLayerTypeName(layer));

function getLayerTypeName(layer)
{
    if (layer instanceof L.Marker){
        return 'Marker';
    }
    else if (layer instanceof L.Tooltip){
        return 'Tooltip';
    }
    else if (layer instanceof L.esri.DynamicMapLayer){
        return 'DynamicMapLayer';
    }
    else if (layer instanceof L.esri.BasemapLayer){
        return 'BasemapLayer';
    }
    else if (layer instanceof L.Layer){
        return 'Layer';
    }
    else
    {
        return 'Unknown';
    }
}

Is there a different approach I should be using?

https://jsfiddle.net/CumminsJP/tb28qwpo/


texts: How to know where views are being generated from in CartoDB?

texts: I see in my dashboard that I can hover over a map and see how much it is being viewed. Many services, like mapbox and even youtube show you where your views are being generated from or perhaps it just tracks links (eg. facebook, twitter, or a particular blog or website). Is this feature available?


texts: Diameter Symbol not appearing in Layer Labels

texts:
I'm trying to label a pipelines layer with Diameter information. I'm using following expression:
u"\u2300" + [Diameter_cm] + "cm"

But Unicode (u2300) appears as a box rather than appearing as a proper Diameter symbol.
Can anyone help?


texts:
I have bbox strategy for one source of data. Code looks like this:

bbox: function newBboxFeatureSource(url, typename) {
    return new ol.source.Vector({
        loader: function (extent) {
            let u = `${url}&TYPENAME=${typename}&bbox=${extent.join(",")}`;

            $.ajax(u).then((response) => {
                this.addFeatures(
                    geoJsonFormat.readFeatures(response)
                );
            });
        },
        strategy: ol.loadingstrategy.bbox
    });
},

I works fine but... When I pan/move the map then this loader is calling again and add another features which fit to new box. But there is a lot of duplicates then because some of new features are just the same as old.
So I wanted first clear all features using this.clear() before add new features but when I add this command then loader is running all the time and I have "infinitive loop". Do you know why? How can I disable loading new features after calling this.clear()?


texts: PostgreSQL PostGIS: append table to another table

texts:
I have a large polygon table  (30+million rows) and a smaller polygon table (45,000 rows). Both are PostGIS enabled.
The smaller table only contains polygons that do not exist in the large table. Both tables have columns named poly_id, ref_no, status and a geometry column (geom). The large table has a few extra columns that aren't important.

My query in pgAdmin is:

INSERT INTO schema.largetable(geom, poly_id, ref_no, status)
SELECT geom, poly_id, ref_no, status
FROM schema.smalltable;

However I get the error:

ERROR: duplicate key value violates unique constraint "largetable_pkey"
DETAIL: Key (objectid)=(1) already exists
SQL state: 23505

the small table doesn't have an 'objectid' column and even if both table has duplicate id values, I don't need the id to carry across from the small table. Can't I just append the small table to the large with new id numbers created for the new additional rows?

Here are the columns in large table:
objectid (pk, integer), poly_id (integer), ref_no(charvar 9), insert_ (charvar 20), update_ (charvar 20), status (charvar 1), shape_length (double), shape_area (double), wkb_geometry (geometry)

small table:
poly_id (pk, integer), ref_no(charvar 9), insert_ (charvar 20), update_ (charvar 20), status (charvar 1), wkb_geometry (geometry)

My 'create table' script for the large table is:

CREATE TABLE schema.largetable
(
    objectid integer NOT NULL DEFAULT     nextval('schema.largetable_objectid_seq'::regclass),
    poly_id integer,
    ref_no character varying(9) COLLATE pg_catalog."default",
    insert_ character varying(20) COLLATE pg_catalog."default",
    update_ character varying(20) COLLATE pg_catalog."default",
    status character varying(1) COLLATE pg_catalog."default",
    shape_length double precision,
    shape_area double precision,
    wkb_geometry geometry(MultiPolygon,27700),
    CONSTRAINT largetable_pkey PRIMARY KEY (objectid)
)
WITH (
    OIDS = FALSE
)
TABLESPACE pg_default;

ALTER TABLE schema.largetable
    OWNER to user_admin;


texts: If rectangle corner points have same nearest neighbor, does whole interior?

texts:
If all four corners of a spherical (but actually WGS84/elliptical) rectangle are closest to a given point p in a set S, does that necessarily mean all points inside the rectangle are also closest to p?

By spherical rectangle, I mean a polygon with coordinates:

{lon1, lat1}, {lon1, lat2}, {lon2, lat2}, {lon2, lat1}, {lon1,lat1}

(with the last point just to close off the polygon if necessary) for two longitudes lon1 and lon2 and two latitudes lat1 and lat2.

These are lines of constant latitude and longitude, not geodesics.

I'm trying to coax Mathematica into creating a geographical Voronoi map (https://github.com/barrycarter/bcapps/blob/master/REDDIT/bc-metro.m). This isn't too hard to do if I assume the Earth is spherical, but, I've decided to do this using the Earth's true shape, or at least WGS84.

Mathematica does have a WGS84 accurate GeoDistance function, but it's expensive to evaluate, so I want to use it as few times as possible.

My plan is to break up my region into latitude/longitude rectangles, and, if the 4 corner points are all closest to the same point in my set, assume the entire rectangular region is also closest to that point.

I'm pretty sure I could prove this on a perfect sphere, but I'm worried that it might not be true on an ellipsoid.

I created http://test.barrycarter.info/gmap8.php using this technique (https://github.com/barrycarter/bcapps/blob/master/MAPS/bc-closest-gmap.pl but the code is currently commented out since I tried to use qhull instead), but I have no idea if it's 100% accurate.


texts:
How many scenes are there for ASTER per year ?

From what I understand is that the temporal is every 16 days. So, some month it will have 2 and some months it will have one.


texts: Getting closest node point (within 2km) from current location using pgrouting?

texts:
I'm very new to the GIS, and grasping the pgrouting function to see the possibles develop this feature.

Basically I wanted to know all the next end points distance within 2KM from the current location.

As in the image, I wanted to know all the D nodes points from the P node point.


texts:
When I try to execute geoprocessing tools (slope, viewshed...) I don't get any results. When I open the results window I see an hourglass icon like in the picture attached near the process.
I'm using ArcGIS Desktop 10.6.1


texts: Does GPS signal quality depend on the chip manufacturer? Does it depend on weather?

texts:
(Sorry if GIS is not the right place to post this question to).

Recently me and my school friends had a GPS software application project to pursue at the university.

A device very similar to this one was used by the application.

One day it refused to show us our current location - it was continuously showing a point in the Atlantic Ocean instead of land [we were outside of any buildings].

I did a double check by opening a default map application on an HTC Desire with Android 2.1 - it showed no current location as well.

Questions:

Does GPS signal quality depend on the chip manufacturer?
Does GPS signal quality depend on weather?
Is it possible to have a commercial hardware-software system which never loses the GPS signal and always shows the correct current position of a user?


texts: Hyperlink from attribute table gets send to browser missing the last character, }

texts:
I am using QGIS 3.10.2-A. I'm working with city data and parcels. I populated my attribute table with a lengthy hyperlink that ends in }. The string is about 455 characters, and it shows in the attribute table complete. However, when clicked, it gets sent to the browser (Firefox, Chrome or Internet Explorer) missing the last }. It renders the hyperlink useless.

I imported the hyperlink to QGIS from a CSV form (no problem)
The details side panel shows it complete When the hyperlink gets clicked,    prompts to an error page I tried cntrl+c/v and still misses the last character.
I can manually add the } to the URL and the URL works fine then.
Other hyperlinks work fine as long as they don't end in }
No actions from action-menu were added or needed with working    hyperlinks.
hyperlink shown in attribute table

https://www.zillow.com/tx/?searchQueryState={%22pagination%22:{},%22usersSearchTerm%22:%22texas%22,%22mapBounds%22:{%22west%22:-98.7578398,%22east%22:-98.7547018,%22south%22:29.79359607,%22north%22:29.79556007},%22regionSelection%22:[{%22regionId%22:54,%22regionType%22:2}],%22isMapVisible%22:true,%22mapZoom%22:18,%22filterState%22:{%22sortSelection%22:{%22value%22:%22globalrelevanceex%22},%22isAllHomes%22:{%22value%22:true}},%22isListVisible%22:true}

Hyperlink sent to browser

https://www.zillow.com/tx/?searchQueryState={%22pagination%22:{},%22usersSearchTerm%22:%22texas%22,%22mapBounds%22:{%22west%22:-98.7578398,%22east%22:-98.7547018,%22south%22:29.79359607,%22north%22:29.79556007},%22regionSelection%22:[{%22regionId%22:54,%22regionType%22:2}],%22isMapVisible%22:true,%22mapZoom%22:18,%22filterState%22:{%22sortSelection%22:{%22value%22:%22globalrelevanceex%22},%22isAllHomes%22:{%22value%22:true}},%22isListVisible%22:true


texts:
How do I ensure the same selection is not made on param[2] and param[3] both of which are populated from the same field using tool validation?

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
     """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()
    self.fcfield = (None, None)

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parmater
    has been changed."""
    if self.params[0].value and self.params[1].value:
      fc, col = str(self.params[0].value), str(self.params[1].value)
      if self.fcfield != (fc, col):
        self.fcfield = (fc, col)

        # Get the unique values of the field 'col' in the feature class 'fc'
        unique_values = [str(val) for val in
                            sorted(
                                set(
                                    row.getValue(col)
                                        for row in arcpy.SearchCursor(fc, None, None, col)
                                )
                            )
                        ]

        # Assign the unique_values list to parameters 2 and 3
        self.params[2].filter.list = unique_values
        self.params[3].filter.list = unique_values

    # Set the default values of parameters 2 and 3 to the first item in the list
    if self.params[2].value not in self.params[2].filter.list:
      self.params[2].value = self.params[2].filter.list[0]
    if self.params[3].value not in self.params[3].filter.list:
      self.params[3].value = self.params[3].filter.list[0]

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""


texts: How to Calculate Field in joined table using ArcPy?

texts:
I have to calculate a field based on another field in a joined table, but an error occurs:

"Parameters are not valid.
ERROR 000728: Field !Table_A.field_A! does not exist within table"

My code:

import arcpy
from arcpy import env
env.workspace = "C:/TMP/BASELINE.gdb"
Table_A = "dist"
Table_B= "matrix"
arcpy.AddJoin_management(Table_A, "id_grid", Table_B, "id_matrix")
arcpy.CalculateField_management(Table_A, "!" + Table_A + ".field_A!", "!" + Table_B + ".field_B!","PYTHON")

Both tables are in the BASELINE.gdb. The Join occurs without any problem.

I forgot to say that, of course, "field_A" exists in "Table_A" and "field_B" exists in "Table_B".

Last but not least, if I run the same calculation in ArcMap, using Field Calculator, no error occurs.


texts: Running my Python script from outside the QGIS Python Console

texts:
I have a Python script that runs in the QGIS Python console but I would like to run it outside of the console, so that a map can be created without having to be within QGIS. The script runs perfectly inside in the console but I just need to find a way how to use it from the outside, so I can then eventually use task scheduler to repeat this process at set time intervals. What are the steps I need to take in order for my script to run outside? I vaguely know that I need import modules/libraries that I have used into Python but I am unsure on how to do this and what other steps I need to take. I have looked at this hand book but I can't seem to find a way what to do - http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html. I have included an example of a small section of my script but included all imports:
import math
import json
import os
import glob
from qgis.PyQt import QtGui

#Add Base Map
urlWithParams = 'type=xyz&url=https://a.tile.openstreetmap.org/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0&crs=EPSG3857'
rlayer2 = QgsRasterLayer(urlWithParams, 'OpenStreetMap', 'wms')

if rlayer2.isValid():
        QgsProject.instance().addMapLayer(rlayer2)
else:
        print('invalid layer')

#Add Vector Layers
vlayer = iface.addVectorLayer(latest_file, "Earthquake Event", "ogr")
if not vlayer:
    print("Layer failed to load!")

# Creating the distance buffers
layer = iface.activeLayer()
feats = [ feat for feat in layer.getFeatures() ]
epsg = layer.crs().postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&    field=y:real&field=point_id:integer""&index=yes"

    #Using zonal statistics to work out affected population sum column
    #######
    # there is no need to look for the raster layer again since you already saved it to a variable before(line 53)
    prefix = rlayer.name()
    params = {'INPUT_RASTER': rlayer.name(), 'RASTER_BAND': 1, 'INPUT_VECTOR': f'MMI {intensity} Buffer', 'COLUMN_PREFIX': prefix+'_', 'STATS': 1}
    processing.run("qgis:zonalstatistics", params)

#This adds a map item to the Print Layout
map = QgsLayoutItemMap(layout)
map.setRect(20, 20, 20, 20)

#Choosing export folder and exporting layout
data_dir = os.path.join(os.path.expanduser('~'), 'Path_to_output')

output_file = feat['id'] + ' Event.png'
output_path = os.path.join(data_dir, output_file)

exporter = QgsLayoutExporter(layout)
exporter.exportToImage(
    output_path, QgsLayoutExporter.ImageExportSettings())

print("Finished")


texts: How to create the attribute table for a raster that was produced by the Slope Tool (ArcGIS10)?

texts:
I use the Slope Tool  (Spatial Analyst Tool – Surface) of ArcGIS10 by ASTER GDEM (that was exported to UTM projection).
The program produce new raster (with slope value), but the resulting raster doesn’t have the attribute table (it has it in inactive form). But it’s possible to create a map from this raster by slope value, it’s possible to represent each pixel of the secondary raster in its unique value, but impossible to open the attribute table.

The same problem has the Aspect Tool (Spatial Analyst Tool – Surface), maybe some others.
Also I can’t convert this raster to a polygon, because Raster to Polygon Tool doesn’t see it as a raster.

How I can make the attribute table of the resulting raster  active?


texts: How to convert length in Decimal degrees to meters programmatically?

texts:
I'm new to qgis and I want to find the length of some roads in a map. So I used the following code,

from shapely.wkb import loads

def getLength():

layer = qgis.utils.iface.mapCanvas().currentLayer()
total = 0
for feature in layer.selectedFeatures():
    geom = feature.geometry()
    wkb = geom.asWkb()
    line = loads(wkb)
    total = total + line.length
return total

print getLength()

It shows the length in degrees I think. Is there any way to convert the result I got to meters, without changing the CRS of project. I want to do it programatically(without adding new field using field calculator).


texts: Problem with nested query in PhP and PostGIS

texts:
I have a simple problem regarding a nested query. I would like to select all buildings in a given circular buffer and create a kml file of those selected objects. The following query does the job very well and returns the results when I run it in PgAdmin.

SELECT total_pop, ST_AsKML(ST_Transform(geom,26918)) AS kmlgm FROM h_buildings WHERE ST_DWithin(ST_Transform(geom, 26918), (SELECT ST_Transform(ST_GeomFromText('POINT(1116858.6062789 7086290.6680572)', 900913), 26918)), 1000)

the "POINT(1116858.6062789 7086290.6680572)" is the center of the circle, and the value 1000 is the radius.

As you can see the query has another SELECT inside it (nested).

The following PHP code is able to generate kml file based on the results of query. I have already tested it with a simple SELECT query and it works fine. However, when I use the above-metioned query it does not work, and it does not give an error.

Does this mean that PHP does not support nested query?!

I'm a beginner, Can anybody please help me to find where the problem is and how it can be solved?

PHP code:

<?php

$dbconn = pg_connect("dbname='mydb' user='postgres' password='****' host='localhost'");
if (!$dbconn) {
    echo "Not connected : " . pg_error();
    exit;
}
$srid = '4326';

$sql = 'SELECT name, ST_AsKML(ST_Transform(geom,26918)) AS kmlgm FROM Buildings WHERE ST_DWithin(ST_Transform(geom, 26918), (SELECT ST_Transform(ST_GeomFromText('POINT(1116858.6062789 7086290.6680572)', 900913), 26918)), 1000)';

 # Try query or error
$query_result = pg_query($dbconn, $sql);
if (!$query_result) {
    echo "An SQL error occured.\n";
    echo pg_last_error($dbconn);
    exit;
}

$kml ="<?xml version='1.0' encoding='UTF-8'?><kml xmlns='http://earth.google.com/kml/2.1'><Document>";

for ($i = 0; $i < pg_numrows($query_result); $i++) {

    $townname = pg_result($query_result, $i, 0);
    $townkml =  pg_result($query_result, $i, 1);

    $kml .= "<Placemark><name>Town: ".$townname."</name><description>".$townname."</description>".$townkml."</Placemark>\n";
}

$kml .= "</Document></kml>";

header('Content-Type: application/vnd.google-earth.kml+xml');
header("Content-Disposition: attachment; filename=sample75.kml");

echo $kml;

pg_close($dbconn);
?>


texts: Select by attribute with field equal to element of vector in ArcPy?

texts:
I am trying to select by an attribute with a condition that field equals to an element of a vector. I use a for loop but could not get the right expression for this condition. Here are my codes.

random37 = [30449, 46617, 17309]
pointlayer = "network_nodeSelected"

for random in random37:

randomname = "n" + str(random)
arcpy.MakeFeatureLayer_management(pointlayer,pointlayer)
arcpy.SelectLayerByAttribute_management(pointlayer,'NEW_SELECTION','"ARC_" = random')
arcpy.CopyFeatures_management(pointlayer, randomname)

print randomname

The error message is presented as follows.

Traceback (most recent call last):

File "", line 6, in
    arcpy.SelectLayerByAttribute_management(pointlayer,'NEW_SELECTION','"ARC_" = point')

File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\management.py", line 7742, in SelectLayerByAttribute
    raise e

ExecuteError: ERROR 000358: Invalid expression
Failed to execute (SelectLayerByAttribute).

As you can see, my problem is the expression '"ARC_" = random'. But I do not know how to get the right expression.


texts:
I would like to get the output of the following command :

gdalinfo --formats > formats.txt

In Linux, it works

In Windows, it works but the output is not completed.

Do you know how to get the complete output in Windows ?

EDIT :
I know it's not complete because it finished in the middle of the line
End of my file :

GSC -raster- (rov): GSC Geogrid
FAST -raster- (rov): EOSAT FAST Format
BT -raster- (rw+v): VTP .bt (Binary Terrain) 1.3 Format
LA

I try different OS, windows 2016, 2019, 10 and differents versions of GDAL, the last one i try is 3.0.0.

EDIT 2 :
I run some test and found out the problem but don't know how to handle it.
So I took the gdal from here :
Gdal with mapserver

From here bin\gdal\apps I copy gdalinfo.exe and put in in bin.

I run gdalinfo --formats > formats.txt and it works

But after that I add a plugin ECW From Here

I extract the msi, I copy paste the gdalplugins and libecwj2.dll to the previous bin folder

And now the commands is not completed. I don't know why.


texts:
I have some regions in an AutoCAD dxf file that needs to be imported in QGIS as polygons.
I tried to import these regions from dxf but nothing happens also using the Dxf2SHP converter.
What can I do?


texts: Geoprocessing describe geometry object in ArcMap Python Window

texts:
I am attempting to identify features in multiple feature layers by their attributes and perform a union on the geometry object.  Essentially pull specific features from different layers and combine them into a new feature. A decomposed version of the function works in the console directly, but for some reason when I run it with a custom Python function imported into the Python Window (ArcMap 10.6.1) and run it there, it does not work.  I've identified the culprit (in the post title) but don't know how to fix it.

def union_multiple(fp_feat, field_name, target_val, fname):
  arcpy.env.workspace = 'path/to/default.gdb'
  for i in range(len(fp_feat)):
    feat = fp_feat[i]
    field = field_name[i]
    target = target_val[i]
    with arcpy.da.SearchCursor(feat, [field, 'SHAPE@']) as cursor:
      for row in cursor:
        if target_val in row:
          if 'shp' not in locals():
            shp = row[1]
          else:
            shp = shp.union(row[1])
  arcpy.CopyFeatures_management(shp, fname_out)

Then in the console I import the function and run it:
from zachs_arc_scripts import union_multiple
fp_feat = ['feat1', 'feat2', 'feat3']
field_name = ['util_type', 'util_type', 'util_type']
target_val = ['pline', 'pline', 'pline']
fname = 'plines_union'

union_multiple(fp_feat, field_name, target_val, fname)

Herein lies the problem - it runs and plines_union is created, but it's only containing one of the feats.  BUT it's not a problem with my code - I've tried.  It's this:
If I add a print(row[1]) to debug I get two different objects between directly running parts of this script decomposed in the Python Window and if I run using an imported function:
Directly in Python Window:
row[1] = <Polyline object at xxx>
From Function in console
row[1] = <geoprocessing describe geometry object object at xxx>
So it's NOT a syntax error, it's something with running via a function that returns the geoprocessing describe geometry object object instead of the Polyline object.  FYI I am utilizing Table of Content feature layers, NOT full file paths.
ArcMap 10.6.1


texts: Google maps marker is not getting displayed

texts:
I am trying to learn google map API. I see a map but I am unable to get a marker. Here is my code

Activity

public class HomeActivity extends Activity {

    @Override
    protected void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        setContentView(R.layout.home_page);
        GoogleMap map = ((MapFragment) getFragmentManager().findFragmentById(R.id.map)).getMap();
        map.addMarker(new MarkerOptions()
                .position(new LatLng(37.7750, 122.4183))
                .title("San Francisco"));
    }

}

home_page.xml

<?xml version="1.0" encoding="utf-8"?>
<fragment xmlns:android="http://schemas.android.com/apk/res/android"
          android:id="@+id/map"
          android:layout_width="match_parent"
          android:layout_height="match_parent"
          android:name="com.google.android.gms.maps.MapFragment"/>

What am I missing out?


texts: Spacing between ArcMap class name and description in the legend

texts: I am wanting to reduce the gap between the class name (crop type) and the description (ac) with in the legend. It appears to be based on the longest class name. is there a way to have it a set distance from each class name instead of following the longest class name. I know I can convert it to graphics and then move them but that is cumbersome and time consuming,, just want to see if there is an easier way.


texts:
Goal: I am attempting to create a fetch (wind exposure) raster for a portion of the St. Lawrence River in order to look at habitat suitability for a number of fish species.

Tool used: http://www.umesc.usgs.gov/management/dss/wind_fetch_wave_models_2012update.html

Tool is zipped near the bottom of page and has a pdf w/ instructions below it.

Link to image of raster input: http://imgur.com/a/gf8Kg#5

After making that album I changed output location to a regular folder.

Wind Direction I am keeping simple to start:
230,100

My result: http://imgur.com/a/7iJst#0

Everything is valued as 0???

Can anyone let me know why I am not getting pretty images like the ones in the PDF?

Additional Information:

-The location is the Thousand Islands at the intersection of Lake Ontario and the St. Lawrence River.


texts: QgsComposition loadFromTemplate Method Causing Segmentation Fault

texts:
I'm attempting to write some stand-alone python code that opens a qgis project and print composer template in order to generate pdf figures. I have this running on my Windows 7 machine, but I need to move it onto our Ubuntu Linux server.

The code crashes when trying to load a template into the print composer, printing the following error to the console:

Segmentation fault (core dumped)

This is the code I'm trying to run using the python interpreter:

from qgis.
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtGui import *
import sys
import os

app = QgsApplication(sys.argv, False)
app.initQgis()

os.chdir('/path/to/output')

renderer = QgsMapRenderer()
project_path = "/path/to/project/project.qgs"
template_path = "/path/to/template/template.qpt"

project = QgsProject.instance()
project.read(QFileInfo(project_path))

root = project.layerTreeRoot()

layers = []
for child in root.children():
    if child.isVisible():
        layers.append(child.layerId())

renderer.setLayerSet(layers)

template_file = file(template_path)
template_content = template_file.read()
template_file.close()
document = QDomDocument()
document.setContent(template_content)

c = QgsComposition(renderer)
c.loadFromTemplate(document)

c.refreshItems()

map_item = c.getComposerItemById('map')
map_item.setNewScale(14400000000,True)

#export to PDF and also save a QGS file
c.exportAsPDF('Fig' + str(1) + '.pdf')
project.write(QFileInfo('Fig' + str(1) + '.qgs'))

#clear project instance
project.clear()
app.exitQgis()
app.exit()
sys.exit()

The line that causes the problem is:

c.loadFromTemplate(document)

The qpt file was created in Windows, I'm not sure if that would be a problem. The qgs file was also created in Windows and it works fine.

EDIT:
I figured out that its the ComposerLabel tag in the qpt file that is causing the problem. Also, I'm using QGS 2.18.


texts: Drawing polygon in ArcGIS Pro not scaling up

texts: I am drawing some street polygons on ArcGIS Pro, when I zoom out while creating a feature, the line is not scaling up. My coordinates are set to UTM, the same with my new feature. Any hint on to fix the scale?


texts: Filling field using row value in ModelBuilder?

texts:
I'm making this model that modifies the table that I show:

MODEL.PNG + TABLE.SHP + TABLE.PNG

(I'm sending it this way because the network I'm using have limited access and this is my only option)

The next step is selecting the values from "FECHA_2" that are lower than the value in "FECHA_MIN" for later calculate another field

One of the ideas I have is to fill all the Null values in the field "FECHA_MIN" with the only value I have, but I don't know how to do it in ModelBuilder, so I'm looking for a way to solve this or another idea to make the selection.


texts: How to copy a legend in QGIS map layout?

texts:
I am creating a map where I have many legend items. Initially, I divided them into two columns but now I want to separate them. When they were together I have customised font size, style etc.

I can create a new legend entry however I am looking for a shortcut where I can copy(duplicate) the existing legends and after pasting I can do the alteration. In items I could see the option of "copy item" but I do not know how to paste the legend for further modification.


texts: Why does CalculateField set values to  when output extent is in a different coordinate system?

texts:
I'm quite new to ArcGIS, and now I have to work on a script that goes through a list of feature classes to:

add a field to their attribute table and
fill that field with a number that I retrieve from somewhere else

When I execute the following statements in the ArcMap Python window or in a minimal standalone script, it works just fine:

table = "JS4A"
arcpy.AddField_management(table, "BK_Code", "LONG")
arcpy.CalculateField_management(table, "BK_Code", "1", "PYTHON_9.3", "")

I can then open the attribute table of that feature class in ArcMap and see that this field has been added and populated with all 1s, just as expected.

Now the weird thing is, when I run the exact same code in the complex production script, the field is added as well, but when I open the attribute table, all the values are <Null>.

I managed to track the problem down to the following line of code in my production script, that comes before the part with AddField and CalculateField:

env.extent = desc.extent

where desc is the return value of an arcpy.Describe() call to some other, unrelated shapefile. All CalculateField_management calls before that line will work fine, all calls after that line produce <Null>. When I replace the above line with a line that sets the output extent to a hard-coded value, it works.

I don't understand this at all. Is this an ArcGIS bug?

Unfortunately I can't omit the above line because it is required in the production script.


texts:
Can I clip a part of a basemap on ArcMap 10.4?

I'd like to know if it is possible to clip a basemap, specifically OpenStreetMap.


texts: Host a WMS service from MapProxy SQLite cache?

texts:
I'm trying to start a WMS server that would not have access to internet and for a source, I have the SQLite layers cache from MapProxy .

I read the manual and saw that I could probably work with a mbtiles cache as a source, but not really SQLite layers databases. It seems to me that some configuration of caches should allow to host the cache itself as a source but I have been unable to get it working. Am I missing something obvious? Is it even possible? Should I be using some other software?


texts: Is ther open source tool for calculating Getis-O-ord (ARCGIS) and able to customize in own platform

texts: i would like to know is there any open source tool for calculating Getis-O-ord(ArcGIS)? If so can i customize for own application or can i write my script within that application like qgis.


texts: Get data from a raster using two coordinate arrays with Python

texts:
I'm trying to find a function in Python, that given two arrays with geographic coordinates, returns a list or array of raster values in those positions. that is to say:

import numpy as np
"""geog coordinates on x axis """
xc = [878571.041, 878583.54109659, 878596.04119319]
""" geographical coordinates on Y axis """
Yc = [187189.14, 187176.64015077, 187164.14030153]
"""list to array"""
X = np.array (Xc)
Y = np.array (Yc)
"""raster values at those coordinates"""
z = function (X, Y)
or
z = function (Xc, Yc)

The output I expect is an array with raster values in the x, y positions. Is there a function in GDAL, Python or any other library that allows me to do this? Or must I necessarily use a for cycle to obtain these values.


texts: Is there an online application to interactively see Openstreetmap tags and keys?

texts: In order to find out how OSM Elements are categorized I was wondering whether there is an online app which shows which elements are displayed and what their respective keys/tags are by a mouse click in the map. I have checked out the OSM Inspector, but this doesnt show all elements in the OSM Map.


texts:
I tried to open a Microstation DGN version 8 file with QGIS/ogrinfo, but GDAL/OGR supports Version 7 only.  So I want to write my own parser by looking at the format of DGN.

Is the Microstation DGN version 8 format specification available publicly ?


texts: edit vertices of intersections from multiple lines

texts:
I would like to edit the vertex where two lines intersect, but doing so makes me manually move each vertex of each line. I tried to merge the lines but it did not work.

Is there an easier way to move the point of intersection without having to move vertices from both lines independently?


texts:
This is the same question as here:
Cartodb - SQL query of public visualization of a private table

We have a private table and would like to use setSQL method to query the sublayer as we did with private data. I have created a basic template per documentation but what's next? Do I have to use this template instead of visualization? I can't get how the concept of the named maps and all that template stuff is related to the actual application code in Javascript.
I have seen that you can have placeholder for SQL in the template but how is this placeholder filled from the application side in this case? We want to query by one or multiple fields. I just generated the code dynamically based on the filter object collecting values from UI elements (field1 and field2 are text and field3 is numeric with min and max values):

var filterFunction = function() {

        ////default query which is TRUE for all rows
        newQuery = "SELECT * FROM mytable WHERE 1=1";

        ////collect filter parameters and append them to WHERE clause in SQL
        if (filter.field1!="undefined"){
            newQuery = newQuery + " AND field1='"+filter.field1+"'";
        };
        if (filter.field2!="undefined"){
            newQuery = newQuery + " AND field2='"+filter.field2+"'";
        };
        if (filter.field3.min!="undefined"){
            newQuery = newQuery + " AND field3>="+filter.field3.min;
        };
        if (filter.field3.max!="undefined"){
            newQuery = newQuery + " AND field3<="+filter.field3.max;
        };

        ////send query to the layer
        sublayers[0].setSQL(newQuery);
        return true;
    };

So, how do we convert this to work with named map? I'm lost.


texts: Moving points to the largest polygon within another polygon

texts:
I need to move a large number of points (street addresses) into the largest polygon (a building, usually a house) within another polygon (taxlot).  The idea here is to put an address point into a house, not a garage or other outbuilding.  Its easy enough to use snap to put it into a building but it needs to be the right one and size seems to be the only option.  Shape_area has been calculated for each building.

How would I go about doing this?  The solutions I've tried wind up putting some points into smaller buildings.

License type = Advanced


texts: Nested Geometric Functions in SLD in GeoServer

texts:
I want to use nested geometric functions in SLD of GeoServer. For example I want to set some offset for certain type of polygons in my data. Then output of the offset geometric function needs to be given as input to another geometric function isometric to give required extrusion. I tried with following code.

    the_geom
    0
    7

  the_geom
7

It gives no validation error but in layer preview, it is throwing the following error:

        java.lang.RuntimeException: Unable to find function isometric Unable to find function isometric

Is it possible to use nested function, or I am doing something wrong here?

Thanks in Advance.


texts: Specify start/stop within color ramp in QGIS

texts:
I'm trying to color values along buffered two-point lines with a color ramp in QGIS. The color ramp should represent the height of the lines, that varies linearly from start to end point. Now I cannot figure out how to make that happen. If I choose the bounds of the color ramp to be the minimum and maximum value of the heights of all features, what would need to happen is that the start value on the color ramp is start_height/(max_height-min_height) and the end value on the color ramp is end_height/(max_height-min_height).
Example: say all heights vary between 0km and 40km. Then for a feature that goes from 13km to 1km, the start point no the gradient should be 32% and the end on the gradient should be 2%.
However, I can't specify the start point and end point on the color ramp.

A workaround would be to use a two color gradient and have begin and end color data defined, but then I would be stuck with a linear interpolation between two colors, which is not what I want.


texts: How to cut a SHP based on another SHP using R

texts: I have a Shapefile of points with certain information. I have another Shapefile with polygons corresponding to the states of the country I am working with (in this case Mexico). Since I only want to work specifically with 3 states, I wanted to know if there is a way to get rid of the points that dont fall into these 3 states.


texts: I'm using ArcGIS Pro 2.4.0 with Windows 10 to edit a few pixels in a raster. I'm creating a map of areas connected to the ocean and there are a few corners of narrow inlets that I need to correct to ensure they connect. The Pixel Editor in the Imagery toolbar is disabled (see picture). I've tried changing the raster format from TIFF to GRID, tried setting the Source Type to "Thematic" or "Elevation", and tried converting from integer to float values. None of that worked, and restarting the program did not work either. Any tricks for getting Pixel Editor to work?


texts: Using If precondition is false, continue iterating next input in ModelBuilder?

texts:
I've built a model in modelbuilder but cannot find an answer to this problem.  Please help if you can.

Firstly I have an 'iterate feature class' tool drawing inputs out of a folder and then selecting attributes from the input table using SQL.  I am then using the "get count" tool (which counts the rows selected) to test if any attributes were selected, this is a precondition to continue the model.  However, if this precondition is 'false' then the model stops.  I'd like for it to just continue with the next iteration of inputs.  Is this possible, and if so, any insight would be great.  I was thinking I might need if-then-else logic but not sure how to use it in this setting.  Please advise.  Here is an image of the model.


texts:
I made some more complex leaflets before, but here I'm stucked. You can find my geojson file here:
https://gist.github.com/Supermarmotte/16e0466c9d17f4ab341b6fdb805bb849
It is displayed correctly and is stored beneath the html file.

    Grrr


texts:
I'm using ArcMap with the Basic license.
I have turnout gates (points) on or very close to irrigation facilities (lines). The irrigation facilities have a field for type (canal, pipeline, culvert, etc.).
What join tool would I use to only add the facility type attribute to my turnout gates feature class attribute table?


texts:
I would like to know if POSTGIS and MySql with Spatial Extension use true geometry to find intersections/overlaps between geospatial objects. But i don't find it precisely in the documentation. I know that doing an ST_Intersects is more precise than only using the && command.

I want to know if in this example the round geospatial objects would intersect/overlap: https://dl.dropboxusercontent.com/u/13064345/example.jpg
The squares are the minimum bounding boxes, they overlap. But i am not sure of the ouput in POSTGIS and MySQl. Hope you can help me.


texts: Importing GEE python coordinate data into QGIS


texts:
I have all my points laid out on my map and i need to connect these points to create polygons. The points are pretty hectic and close together, but I need them created a certain way, so i was wondering how I can do that. I've tried Points2polygon and while it does create a polygon its not what im looking for.

I just think it would be easier to manually connect the dots because I know what they're supposed to look like.


texts: Performing a Distance Matrix based on specific attributes

texts:
I have two point layers with a field called names (A1, A2, A3....) in both layers. I want to perform a distance matrix for those points, but just only for those with the same attribute 'name'. If  I run the tool distance matrix it takes too much time due to the high amount of data.
Edit: I did some work around by creating connecting lines and then calculate the length, but maybe there is a smarter solution.


texts: One to Many Relationship in QGIS with PostGIS

texts:
Does anyone know of any tool that can do a One-to-Many relationship in QGIS?

I have my data in a PostGIS database.  I would like to be able to query spatial layers that are related to spatial or non spatial tables and vice versa.  I believe there was a way in ArcGIS 9.x that allowed to do something similar.


texts: Polygon to Point Grid with the same information

texts:
I want a grid of points (with defined spaces between them) inside a Polygon. All the points should contain the same information as the original Polygon did.
I tried the centroid (only one point -> not good) also the random point (no information of the original Polygon -> also not good). Is there a simple way to do this in Qgis?
To show my thought-progress I added my poorly made picture of paint. The black circle represents the Polygon and the red dots should be my Point-grid.


texts: How do I determine programmatically what values to add to the proj4js constructor for a CRS or SRS such that I can iteratively render proper overlays

texts:
I think this is primarily a GIS based question as far as for example:

var crs = L.CRS.proj4js('EPSG:2154',
                    '+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
                    transformation);

So if I read in a wms tile layer through the leaflet api to create an overlay on a wfs base map how do I decide how to adjust those lat/lon values? And what is ellps, towgs, and units?


texts: I know this is simple but am unable to find this online.  Using ArcMap 10.5 I want to drop the first digit of the three in a numeric column.. e.g., (023) where I want it (23).  Which query can I use to make that happen?


texts:
I am working with a point shapefile, for which I have created a new field NEW_ID. This represents a unique code for each point. This is more convenient to me than using the OBJECTID.

I want to fill the new field with the values coming from the following list:
myvalues=range(1000000,1002332)

And I am trying to get this working, but the field is still <Null> after the snippet has run:

myFeatureClass=r'C:\temp\testGDB\test_shp'
cur = arcpy.UpdateCursor(myFeatureClass)

for row in cur:
    row.setValue('NEW_ID', row.getValue(myvalues[row]))
    cur.updateRow(row)

I suspect I cannot use row as an iterator. So, what should I use instead?

ArcGIS version: 10.3.1


texts: Get elevation of GeoTIFF using GDAL bindings in Java

texts:
I'd like to read an elevation value from a GeoTIFF that I have. It's already working, but I'd like to interpolate the values for better results.

I currently do so like this:

Dataset dataset = gdal.Open(pathToGeoTiffFile, 0);
String projection = dataset.GetProjection();
double[] transform = dataset.GetGeoTransform();

SpatialReference src = new SpatialReference();
src.SetWellKnownGeogCS("WGS84");
SpatialReference dst = new SpatialReference(projection);
CoordinateTransformation coordinateTransformation = new CoordinateTransformation(src, dst);

double[] xy = coordinateTransformation.TransformPoint(lng, lat);

float[] result = new float[1];
int[] band = {1};

int x = (int)(((xy[0] - transform[0]) / transform[1]) - 0.5);
int y = (int)(((xy[1] - transform[3]) / transform[5]) - 0.5);

dataset.ReadRaster(x, y, 1, 1, 1, 1, 6, result, band);

double elevation = result[0];

This basically works, but there is no interpolation. So basically a coordinate gets mapped to a pixel (no matter if it's close to the next pixel or not) and then I take that value.
But I'd prefer to interpolate the value. So if I am on. f.e the bottom of a pixel and the next pixel underneath has a different elevation, the value that I return should be interpolated.

I could of course calculate that on my own, but I was wondering if there is a function that I am missing?


texts: How to keep NULL values in attributes until the user choses from the value map?

texts: When I set up a attribute domain by using value map in QGIS then create a feature which in turn pops up the attribute pop up window - the value mapping for the attribute inserts the first value from the mapping into the attribute by default.  I would like the value mapping to leave the attribute as null UNLESS an analyst chooses to click on the attribute in the attribute pop up window to reveal the value mapping pull down AND THEN choose a value for the attribute.   Does anyone have an idea how to do that?


texts: AttributeError: __exit__ when trying to use arcpy.UpdateCursor

texts:
I'm trying to populate a field from values of an existing field in a feature class using arcpy.UpdateCursor.

I'm getting this error:     AttributeError: __exit__

Here's the chunk of code I'm using:

with arcpy.UpdateCursor(polyline, ["SEWER_TYPE"]) as cursor:
    for row in cursor:
        if row.getValue("Layer") == "FUTURE_STORM":
            row.setValue("SEWER_TYPE", "STORM")
            cursor.updateRow(row)
del row, cursor


texts: Geoserver import does not recognize projection

texts: Using GeoServer 2.8 and the Data importer plugin I can import shapefiles generated in ArcMap with WGS84 geographic coordinates (EPSG:4326) and UTM+WGS84 projection (epsg:32724). As a reminder, ArcMAP will not add the epsg code to the prj file. That is, it appears that GeoServer can identify the EPSG just based on ESRIs WKT projection file. However, if I try UTM+SIRGAS2000 (epsg:31984) the importer does not recognize the SRS and I have to manually select it. Now, if I substitute the prj file with a more complete one (which includes the EPSG code), geoserver will use it. Is there a way to make GeoServer accept the prj generated by ArcMAP and correctly identify the SIRGAS2000 projection?


texts: Losing hyperlink between ArcMap and Google Earth

texts: Most of my colleagues use Google Earth to check our infrastructure. A few years before I started here, the hyperlinks from our GIS to Google Earth lost connection. The path is still in the "plan_path", but now the users on Google Earth have to copy/paste the link to the file explorer. My predecessor couldn't figure out, and I can't either, so I decided to ask here. I attached a couple screenshots.


texts: QGIS Atlas - Making legend and size assistant automatically update based on which area selected

texts:
I'm hoping to use the atlas feature to generate a batch of maps for London postcode districts which will show supply of private indicated by size of dots within the PCD. I was wondering if it was possible to make it so the size of the dots would automatically adjust on each map generated by the atlas based on the supply within each district.

For example, in an area such as Mayfair, where number of units per development is very low, the sizes would need to be scaled down (0-10,10-20, etc..) differently to Battersea, where there are a lot of developments with larger number of units (0-100,100-200,etc..)


texts: POSTGIS: function st_collect(geography) does not exist

texts:
I am trying to use some Geography functions, but I am getting an error that (all) functions are not existing. I have checked other Stackexchange posts (e.g. this one), but the solutions did not work for me.

Since no function is available, I assume that something is wrong with the geography tables. I am using POSTGIS 2.0.1, therefore the geography functions should be supported.

database# select postgis_full_version();
POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL
1.9.1, released 2012/05/15" LIBXML="2.7.3" LIBJSON="UNKNOWN" TOPOLOGY RASTER

I can see the geography tables in the database and the database owner also owns the geography tables.

 public   | geography_columns | view     | developer
 public   | geometry_columns  | view     | developer
 public   | raster_columns    | view     | developer
 public   | raster_overviews  | view     | developer

What could cause the error below?

function st_collect(geography) does not exist
LINE 1: SELECT ST_xMin(ST_Collect(location)) AS id, ST_xMax(ST_Colle...
                       ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.


texts: Rounding values of a raster dataset in a GetFeatureInfo response

texts:
Is there a way to round decimals of floating values for a GetFeatureInfo request from a raster (geotiff) in Geoserver? I would like to round up values to 2 decimals instead of 15. Is this something I have to do in the SLD? I generated a SLD from QGIS which now looks something like this:

            <sld:FeatureTypeStyle>
                  <sld:Rule>
                    <sld:RasterSymbolizer>
                        <sld:Geometry>
                        <ogc:PropertyName>grid</ogc:PropertyName>
                        </sld:Geometry>
                        <sld:ColorMap type="intervals">
                            <sld:ColorMapEntry color="#4575B4" label="4,0 - 4,5" opacity="1.0" quantity="4.5"/>
                        </sld:ColorMap>
                    </sld:RasterSymbolizer>
                </sld:Rule>
            </sld:FeatureTypeStyle>

I saw a similar question asked here: Removing Trailing Zeros in Geoserver SLD but I couldn't figure out how and where to implement this part:

<Function name="numberFormat">
   <Literal>#.##</Literal>
</Function>

Usually I do rounding in PostGIS for vector data, but I am not sure how to do it with GeoTiffs.


texts: Get the maximum distance between two consecutive points in a linestring

texts:
I want to check if a track is valid by measuring the distance between each consecutive points of my track and discarding the tracks which have two consecutive points really far away. My tracks are stored in a linestring. Do you know if its possible to do that?

input: LINESTRING( P1,P2,P3 .. Pn ) size n

output: Max Distance between each consecutive points


texts: Projected Coordinate System for COGO descriptions?

texts:
I am COGOing a legal description from 1950.  When I COGO the legal description it does not match the boundaries from the data in our database.  I suspect it is the coordinate systems not matching.  The projection my company uses is NAD_1983_StatePlane_Utah_Central_FIPS_4302_Feet.

Is there a way to determine what the coordinate system would be from a document from the 1950's?


texts: What happens if some polygons are empty in reduce() in Google Earth Engine

texts: I am drawing the NDVI time-series trend for two types polygon in my target area using Sentinel-2. I don't know why I miss one day's NDVI mean. Is it because some polygons are unmasked after cloudrmove? The second picture describe the image in that day.(12.22.2017)


texts: Which training courses in Geodetics can you recommend?

texts:
Can anyone recommend any training courses in geodesy which would cover the following objectives:

Case studies on coordinate issues.
Overview of geodetic datums, coordinate reference systems and transformations.
Use of Global navigation satellite systems ie. GPS.
Methods in projection methods and methods in measuring the Earths surface.
Vertical datums, geoidal models, vertical CRS and their transformations.
Orthorectification, georeferencing and other data referencing methods.


texts:
I have white polygons on same layer and I want to fill the space around the polygons (i.e. the whole background) with a single black polygon. The result would be a black & white map.

Is there anyway to do it by using rasters instead?


texts: QGIS Composer Legend not showing all map items

texts:
I have 4 focus areas on a map... They are categorized by name, but the symbology is the same for all 4. When I add them to the print composer only the first 3 focus areas show up in the legend...

I'd like all 4 to be there as I want to have the focus areas named in the legend.

See attached screenshot.


texts: Working with NaN/null and different extents in QGIS Raster calculator?

texts:
When it comes to NULL value and different extents of rasters I find it extremely cumbersome to use QGIS raster calculator. As far as I am aware there is no way to operate with NaN/Null value. So no matter what you do NaN/NULL in any input become NULL in output.

As far as I am aware you can't do anything like:

(raster1 is NULL)*raster2+(raster1 is not NULL)*raster1

I have similar issue with extents of raster. I have one raster of large extent and second covering only small part of it. I need to update null/zero values of large extent raster to values in small extent raster. If I do calculation below everything off of small extent becomes NULL/Nan.

(largeraster = 0)*smallraster + (largeraster != 0)*largeraster

Is there any way to be able to operate with NULL value in raster calculator? Is there any way to not trash all values in areas where all rasters do not overlap?


texts: I traced sidewalks lines from parcel data provided to me, but now my boss wants the lines shifted closer to the street. Is there a way to do this in bulk, and keep the attributes? I have 1150 lines and I'd rather not have to do each one individually. I know I'll probably need to touch up the corners again, which I am fine with. Some need to be shifted to the left and some to right. I know you can use the Copy Parallel tool, but I can't seem to keep the same attributes.


texts: Satellite moving target tracking simulation

texts:
I am trying to create a Simulation of an earth-pointing satellite.
The satellite should track a moving target on the ground via an on-board camera.

How can I calculate the goal attitude of the satellite?

I assume I have every data needed.
XY Point in image,
Camera data (FOV...)
Rotation and transformation of the satellite...

I am especially interested in detailed literature for this topic.
Till now I only found papers explaining control algorithms, but
no explanation of getting the goal attitude.

Another idea would be, to control the satellite's attitude by the target XY pixel offset
from the middle of the picture. I am not sure if this is possible. I couldn't find any
literature about it.


texts: qgis python + layout/template + exporting in .png + filtering layers

texts:
I need to do a lot of pictures at different points from vector layer on each. And the picture (.jpg / .png) must be signed and so on - then the layout is best suited. Now I do a picture on a map posted on the points (and even get .png) how to create a layout with the data or work from the layout and apply it using a filter python?

My code:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
import qgis.utils
layer = iface.activeLayer()
layer.setSubsetString('')
iter = layer.getFeatures()
for feature in iter:
    name_dot = feature['name']
    saving = "C:/project_qgis/"+"%s"%name_dot+".png"
    print saving
    layer.setSubsetString('"name"=\'%s\'' % name_dot)
    mapRenderer = iface.mapCanvas().mapRenderer()
    c = QgsComposition(mapRenderer)
    c.setPlotStyle(QgsComposition.Print)
    x, y = 0, 0
    w = c.paperWidth()
    h = c.paperHeight()
    composerMap = QgsComposerMap(c, x ,y, w, h)
    c.addItem(composerMap)
    dpi = c.printResolution()
    dpmm = dpi / 25.4
    width = int(dpmm *w)
    height = int(dpmm *h)
    image = QImage(QSize(width, height), QImage.Format_ARGB32)
    image.setDotsPerMeterX(dpmm * 1000)
    image.setDotsPerMeterY(dpmm * 1000)
    image.fill(0)
    imagePainter = QPainter(image)
    sourceArea = QRectF(0, 0, w, h)
    targetArea = QRectF(0, 0, width, height)
    c.render(imagePainter, targetArea, sourceArea)
    imagePainter.end()
    image.save("%s" %saving, "png")

In fact, the biggest problem is the size of the image, scale, degree grid and the label on the finished picture that can be done during layout...


texts: Import map server data: unable to edit data in ArcMap

texts:
I am trying to import the data from the below map server but have been having some trouble. I am trying to work with the Census 2015 (EA) (0) layer on the map server. I am able to use the option View In: ArcMap, which allows me to download a .lyr file that I can open in ArcMap. This file appears to have the map that I would like to use. However, I am unable to save this layer to a shapefile that I am able to edit. How do I go about getting a version of the map from the map server from a .lyr to a .shp that I can edit?

MapServer: https://www.integems-gis.solutions/server/rest/services/ms/sle_admin_boundaries/MapServer


texts: How to assign the attributes of the largest overlapping polygon area to each overlapping gridcell?

texts:
I am new to R, and I am trying to use it to calculate the areas of polygons that overlap with each gridcell, in order to assign the attributes of the largest overlapping polygon area to the corresponding cell and build a dataframe, in which each gridcell represents a row, and attributes represent columns.

.

As shown in the figure above, I have a polygon shapefile representing indegenous lands and an overlayed gridcell (polygons).

Is there anyway I could do this?
I have seen this question/answer: How to measure the area of polygons within a fishnet but for some reason the intersect function in QGIS (which is the GIS software I am using) does not generate a new layer.

Just in case, here you could download my spatial data: https://www.dropbox.com/sh/hdbpcq80uzd1ben/AADdG5Lg-BBbvFwA3QRLycSRa?dl=0


texts:
I was used to download SRTM/DEM data from the following link: https://dds.cr.usgs.gov/srtm/version2_1/SRTM3/
But I'm not being able to download anymore, I don't know the exact causes of that, it was really practical since the images were all segmented by regions.
Does anyone have any idea why the link is down? Is there anywhere else where data can be acquired easily just clicking on links and downloading it?


texts: Is there a way to set the default coordinate precision to 7 for GeoJSON files in QGIS? This would help me save time when I go to save files.


texts:
I want to create and display a graph on Main Window which I designed in QtDesigner. These plot graph values are taken from table widget, it is not problem. Main problem is here, what I use to show the graph on Main Window. When I click the "Çiz" button I want to give the plot values and show the graphic. My Main Window is here:

"Hacim" and "Yüzey Alanı (m2)" headers are the columns for values to graph. Above the "Çiz" Button, there is a gap to show the graph.

EDIT

I look every topics about the issue but I can not handle it. I try to use it but it gives me the error. The code is:

    self.layoutw = QtWidgets.QVBoxLayout()
    self._dlg.widget.setLayout(self.layoutw)
    self.addPlotWidget()

def addPlotWidget(self):
    x = np.random.normal(loc=0.0, scale=2, size=10)

    widget = pg.PlotWidget(title="Some pasaasa")
    widget.plotItem.plot(x)
    self.layoutw.addWidget(widget)

self._dlg.widget is here QWidget that I added from QtDesigner.

When it runs, gives me a plot like that:

When I move the cursor on the plot or something else - I can nor realize it clearly - it gives me this error or other one:

:


texts: How to query breadwinning mothers in IPUMS USA microdata?

texts:
Overall goal: Map rates of breadwinning mothers throughout South Carolina.

I haven't been able to find exactly what I need in American Factfinder. I am looking for something that comes close to what this Pew Research article did.

Since they mentioned IPUMS microdata, I figured I would try it but I am having a hard time figuring out how to query for breadwinning mothers.

For now I am using IPUMS SDA to query and analyze the data online, here.

For the rows, I selected my geography, all PUMAs in SC, so cpuma0010(915-924) has been working so far.

However, I would like the columns to show something close to the following:

Median income of married women who earn more than their husbands
Married women who are the sole provider for their family (with or without children)
Median income for single mothers

Which codes do I need to select to get this to work. I have been querying over and over in SDA and still haven't found quite what I need. I plan to then join this data up to PUMA areas and map it.

I tried asking on their own forum, and it is not working, their forum is timing out.


texts: Creating linestring layer in epsg:3857 from geojson data in 4326

texts:
In OL3, I am confused about creating linestring layer in epsg:3857 from geojson data in 4326.
What should I do?


texts: arcpy.da.Walk not adhering to datatype filter "Locator"

texts:
I have a file gdb containing feature classes and address locators. I am listing and storing the names of the datasets using

arcpy.env.workspace = "C:\input.gdb"
(root, fds, ftrs) = arcpy.da.Walk().next()

The feature class names appear in ftrs as expected, but the locators do not. By filtering on FeatureClass

(root, fds, ftrs) = arcpy.da.Walk(datatype='FeatureClass').next()

the feature class names are populated once again, as expected. Filtering on datatype=Locator returns an empty list. I set the filter to Table as a check, and the list returned empty as expected.

I would like to know why the address locators in the file gdb are not being listed when arcpy.da.Walk is used to specifically filter out Locator?


texts: Area and dimensions of ee.Geometry.Rectangle

texts:
I am using Google Earth Engine to analyze small areas and defining an square around each point. I need some help understanding how to get the correct measurements. I wanted to know if there is way to

Make a rectangle or square passing specific dimensions in meters. In the code below, 1E-4 is not the same as 21 m.
From a rectangle, calculate width and length in meters. I used Rectangle.Perimeter and Rectangle.Area which give me the correct results but Rectangle.Length returns zero. I also thought of looking at the coordinates of the rectangle to calculate the width and length but the coordinates seem to be the same (shown in figure below)

Here is my code. I know GEE returns values in meters which are very close to what I expect to see.

var makeRectangle = function(point, xDistance, yDistance, proj){
  var geometry =
    ee.Geometry.Rectangle(
      [point[0] - xDistance,
       point[1] - yDistance,
       point[0] + xDistance,
       point[1] + yDistance], proj, false);
  return geometry
}

var p = [-80.2081553905,25.8888877858]
var xDistance = 1E-4
var yDistance = 1E-4
var rect = makeRectangle(p, xDistance, yDistance,'EPSG:4326')
Map.addLayer(rect)
print(rect.area(1E-3))
print(rect.perimeter(1E-3))
print(rect.length())
print(rect.coordinates())

Map.setCenter(p[0], p[1], 21);
Map.setOptions('satellite');
Map.addLayer(ee.Geometry.Point(p), {color: '#FFFF00'});


texts: ArcGIS Desktop 10.3 - definition query to exclude integers starting with certain values

texts:
I have an integer field in a polygon shapefile that I want to definition query out everything number that beings with 18.  These numbers aren't sequential and are of any length, ie I want to exclude 18, 181, 18111, etc but not the numbers in between that don't start with 18.

After some research and going over this SE thread I'm not sure this is possible.  Any ideas?


texts: storing date range in PostgreSQL/PostGIS for time parameter in WMS/WFS

texts:
Using GeoServer I want to deliver some objects/features with WFS and WMS. The objects and their attributes have a date range during which they are valid. Say, we have street lamps and each street lamp has an ID, POINT (geom), start_date, end_date.

The WMS and WFS request can use the time parameter to filter the resulting objects (see WFS GetFeature with time, WMS GetMap with time). Using the time filter, only the street lamps that a within start_date and end_date would be returned by GeoServer .

However, I now saw that for time the respective time column must be stored with a date type. That is, in geoserver one cannot configure start or end date columns for a WMS or WFS layer. It is only possible to configure a single date column (same for MapServer).

What's the best way to store my street lamps in PostgreSQL/PostGIS? For a lamp that was created on Jan 01, 2014 and destroyed on Jan 20, 2014 would I need to have 20 data sets in my database:

ID GEOM DATE
01 xy   2014-01-01
01 xy   2014-01-02
...
01 xy   2014-01-19
01 xy   2014-01-20

Is there a way to store date ranges in PostgreSQL/PostGIS that can be used by GeoServer? Or is this impossible?


texts: PyQgs: how to get a list of all textlabels in my Layout

texts:
I have a print layout with a bunch of items, two text labels among them. For the pluging I'm wrting right now I need a list with the two text labels. My initial idea was
labels = QgsProject.instance().layoutManager().layouts().QgsLayoutItemLabel().currentText()

However when I try to execute my plugin I get the following python error:

  AttributeError: 'list' object has no attribute 'QgsLayoutItemLabel'

Any ideas how to get a list with the text labels?


texts: How to set the color of a feature depending on field attributes with PyQGIS?

texts:
from PyQt4.QtCore import *
from PyQt4.QtGui import *
def Symbology():
    canvas = qgis.utils.iface.mapCanvas()
    layer = qgis.utils.iface.activeLayer()
    for feat in layer.getFeatures():
        Value = feat['field']
        print Value
        if Value != '':
            featColor = {
                  Value :  ('red',Value ),
                    }
        elif Value == '':
             featColor = {
                  '' :  ('green','')
                    }
        ranges = []
        for rType, (color, label) in featColor.items():
            print rType,(color, label)
            symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
            symbol.setColor(QColor(color))
            category = QgsRendererCategoryV2(rType, symbol, label)
            ranges.append(category)

    # create the renderer and assign it to a layer
    expression = 'field' # field name
    renderer = QgsCategorizedSymbolRendererV2(expression, categories)
    layer.setRendererV2(renderer)
    canvas.zoomToSelected()
    layer.removeSelection()

Here i have to set based on the feature value  but when i am executing this script it is setting for only one value. I f i want to set colors based on the values how can i do that?

The value contains for example Value = ['x','y',NULL]


texts: Transfering complex map symbology from ArcGIS Desktop using AutoCad Map object classes?

texts: I'm trying to find the best way to transfer maps form ArcGIS to AutoCad Map. In ArcGIS map model is  separated from it's representation. View is created by applying symbology and representation rules to  underlying model. As far as I know this approach could be achieved in AutoCad Map using object definition classes. So I looking for a way to create complex topographic map symbology(like in ArcGIS) using AutoCad Map object classes. Official AutoDesk help describes this very poor. Is it really possible? Could someone provide links to some articles, how-tows about this issue?


texts: I'm a new user of QGIS and OpenLayers. I had used it some time ago, then I had to format my PC and reinstall all programs. After this operation, I've found I cannot use OpenLayers because the following window appears.


texts: Choosing interpolation method for grid data?

texts:
I have a grid data of evenly spaced points that gives the number of people with a facebook account within a certain radius (approximately equal to the mid-distance between to points). From that I'd like to draw a smoothed map of my variable density of facebook users.

I read that available interpolation methods for that were :

natural neighbours
kriging, or IDW in the simpler way
rectangular interpolation

1) Kriging seems to be a good idea in any case, still I wonder if rectangular interpolation could be more appropriate and far less complex since I have perfectly evenly spaced values ?

2) Are there any interpolation models that takes also as input the radius around each data point ? From what I understand all of the above methods use only X-Y coordinates and the value of the interest feature, not any radius around X-Y in which the value is the same.

From @Spacedman's comment, I see that I misunderstood my dataset. What I actually have is z values for a set of discs defined by a set of points (X,Y) and a unique radius r*. Which means, at that point I don't have a true grid to interpolate on.

An easy workaround is to assume that the density at each point (let's call them f(X,Y) by contrast with z(X,Y,r*)) is uniform in my disc, so that f(X,Y) = a in disc A. Since I know that z(A) = a*pi*r^2, I know that a = z(A)/(pi*r^2). I know have values for a set of points, which are still consistent with the values of the disc they belong to.

What if the discs overlap ? I can't assume anymore that f(X,Y) = a around a whole disc since I may have different values for a same (X,Y). Would taking the mean be outrageous ?