texts: ArcGIS Operations Dashboard - Display metadata in rich text

texts:
Is there a way to display with a rich text the content of the dashboard metadata ? Maybe using Arcade ?
Specifically I would like to display the description in the item's metadata somewhere in the dashboard (so editing it in the metadata automaticall updates the display on the dashboard) :
Metadata:

Dashboard:


texts: How to improve the horizontal and absolute accuracy of UAV layer

texts:
I am working with two digital elevation model data sets: LiDar and UAV DEMs. Both flights were conducted in the same location along a coastal shoreline showing the intertidal area with lots of bare rock and sand. The UAV layer is older and has poor accuracy compared to the LiDar file which is a newly conducted flight. The UAV digital elevation model is inaccurate by ~5-10m. The photos below show the logarithmic error of the UAV DEM.
SW-NE cross-section over the figure above, showing a logarithmic increase in bias from West to East:

NW-SE cross section, showing not as strong but still considerable gradient:

Is there a tool in QGIS, Global Mapper, or Pix4D that can modify this UAV layer to improve the accuracy?
I have tried the 'rubber sheeting' method in both QGIS 3.1 and Global Mapper, but did not have any luck. My assumption about the UAV DEM is there were issues during the pre-processing stage, perhaps with the flight overlap % or the qualities of the GCPs used to process it. I don't have access to the older files so cannot go back and add anything to the project to improve the model, unfortunately.
I am new to remote sensing, and have yet to work with python, CAD, or postgis.


texts:
I am brand new to postgis. I just learned that geometries model simple cartesian objects and that geographies model shapes on the surface of the earth. I have a bunch of data in columns X and Y of a table -- where X and Y are state plane coordinates in a medium-sized US city.

It seems  like you use ST_Point() to get the coordinates into a point and then ST_SetSRID to get them into a certain coordinate system and then ST_Tranform to get them into a lat/long format. But these are points on the surface of the earth -- so I really want a want a geography not a geometry, right?

Am I describing the correct method to do this? Do you do ST_Transform on geometries only?


texts: OSM tile server -- too many highway labels

texts:
I'm running a local OSM tile server.  I've combined three pbfs (three neighboring states) using osmconvert, then loaded them into postgis database:
osm2pgsql --slim -d gis --hstore --multi-geometry --tag-transform-script openstreetmap-carto/openstreetmap-carto.lua --style openstreetmap-carto/openstreetmap-carto.style
combined_areas.pbf

then generated the style.xml carto project.mml > style.xml I did not make any changes to this style file.  The style looks the same as the standard style on openstreetmaps.org (which is what I want) except that there are way way too many road labels, especially at lower zoom levels. How can I fix this?


texts: Problem to convert geometry to coordinates of geometry(MultiLineString) data


texts: Create additional contact details in GeoServer WMS capabilities document

texts:
I'd like to have separate contact details in the GeoServer WMS capabilities for the provider of the service and the provider of the contained data.

When setting up a GeoServer WMS, it adds the general contact details at the  beginning of the Capabilities document:

<WMS_Capabilities>
    <Service>
    ...
    <OnlineResource/>
    <ContactInformation>
        <ContactPersonPrimary>
            <ContactPerson/>
            <ContactOrganization/>
        </ContactPersonPrimary>
        <ContactPosition/>
        <ContactAddress>
            <AddressType/>
            <Address/>
            <City/>
            <StateOrProvince/>
            <PostCode/>
            <Country/>
        </ContactAddress>
        <ContactVoiceTelephone/>
        <ContactFacsimileTelephone/>
        <ContactElectronicMailAddress/>
    </ContactInformation>
    ...
    </Service>
    ...
</WMS_Capabilities>

Is it possible to add another contact information to the capabilities document?

When enabling a workspace, the provided details overwrite the general information but do not get added separately.


texts: How to change the projection and resolution of a file using gdal?

texts:
I have a file with these information:

lat_min -60.0
lat_max 80.0
lon_min -180.0
lon_max 179.9
lon_sampling 0.1 degree
lat_sampling 0.1 degree
samples = 3600
lines   = 1401

I want to re-project the file to

I tried

 gdalwarp -of "ENVI" -t_srs `EPSG:3410`  -ts 1383 586 but give error
 gdalwarp -of "ENVI" -t_srs "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs" -ts 1383 586 -ot Float32  file.envi file_1383586.envi
ERROR 1: No PROJ.4 translation for source SRS, coordinate transformation initialization has failed


texts: Generalising Rasters using .ovr pyramid data?

texts:
I have to produce some lower spatial resolution rasters from a large number of very big geotiffs (c.800GB).

It occured to me that since I've produced the pyramid overviews using GDAL already  I could do this by using the data in the *.ovr files without incurring much additional processing.

Does anyone know of an elegant way of doing this using GDAL or Python?


texts:
I really am at a loss with trying to figure out how to parse only walking routes using OSM2PO.

I've modified the config file to add foot to the routes which can be walked along, as follows:

wtr.tag.highway.motorway =       1,  11, 120, car
wtr.tag.highway.motorway_link =  1,  12, 30,  car
wtr.tag.highway.trunk =          1,  13, 90,  car
wtr.tag.highway.trunk_link =     1,  14, 30,  car
wtr.tag.highway.primary =        1,  15, 70,  car|foot
wtr.tag.highway.primary_link =   1,  16, 30,  car|foot
wtr.tag.highway.secondary =      1,  21, 60,  car|foot
wtr.tag.highway.secondary_link = 1,  22, 30,  car|foot
wtr.tag.highway.tertiary =       1,  31, 40,  car|bike|foot
wtr.tag.highway.tertiary_link =  1,  32, 20,  car|bike|foot
wtr.tag.highway.residential =    1,  41, 40,  car|bike|foot
wtr.tag.highway.road =           1,  42, 50,  car|bike|foot
wtr.tag.highway.unclassified =   1,  43, 50,  car|bike|foot
wtr.tag.highway.service =        1,  51, 5,   car|bike
wtr.tag.highway.living_street =  1,  63, 7,   car|bike|foot
wtr.tag.highway.pedestrian =     1,  62, 5,   car|bike|foot
wtr.tag.highway.track =          1,  71, 10,  car|bike|foot
wtr.tag.highway.path =           1,  72, 10,  car|bike|foot
wtr.tag.highway.cycleway =       1,  81, 15,  bike
wtr.tag.highway.footway =        1,  91, 5,   car|foot
wtr.tag.highway.steps =          1,  92, 5,   car|foot
wtr.tag.route.ferry =            2,   1, 10,  ferry
wtr.tag.route.shuttle_train =    2,   2, 50,  rail|car
wtr.tag.railway.rail =           3,   3, 50,  rail

I've also changed

# only convert ways with one of these flags

wtr.finalMask=foot

As opposed to being

wtr.finalMask=car

The problem is that Osm2po is still not converting routes which I've specified it to convert by defining those types of route as "foot". I'm getting an output which contains only residential streets for some reason.

Is there something else which needs changing in the configuration file that I'm missing?


texts: given sets of vehicle passing on each edge, how to get most used vehicle to avoid transition?

texts:
I have two tables: one, the resulting table from converting osm.xml via osm2po so I can use pgRouting, and another one that contains information on what edges each public vehicle passes through (|vehicle_code|edge_ids|). Note: one edge_id might contain many public vehicle codes.

So after executing pgRouting shortest path query, I want to know which vehicles are the best to use to navigate the shortest path to avoid too many transitions.

I'm kinda new to SQL; please help. Thanks!


texts:
I want to fill  a new field in a shapefile using arcpy. I have differents values so I am using for iterator to run all rows.
But when I run the code all rows are filling with the latest row value.
The code is:

#-----------------------------
camp= "GROUP"

a= parameters[0].valueAsText

arcpy.AddField_management(a, camp, "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
#-----------------------------

cursor = arcpy.SearchCursor(a)

for row in cursor:
    n = row.getValue("CODE2012")
    arcpy.CalculateField_management(a, camp, n,"PYTHON_9.3",code_block="")`

I obtain this attribute table:


texts: Getting error when clipping shapefile to shapefile using ogr2ogr?

texts:
I want to clip a shapefile to a shapefile. The input shapefile is a 1 m vector grid generated in QGIS. The tile I'm clipping to is just a simple polygon (the geometry is contained within the 1 m vector grid).

I put into the command line:

  ogr2ogr -overwrite -t_srs EPSG:27700 -clipsrc
  ./Crewe_boundary_tile_grid_ID__0.shp
  ./Crewe_boundary_1m_grid_clipped_ID__0.shp
  ./Crewe_boundary_1m_grid.shp

However I get the error message:

  ERROR 1: Attempt to write non-polygon (LINESTRING) geometry to POLYGON
  type shapefile. ERROR 1: Unable to write feature 500 from layer
  Crewe_boundary_1m_grid.

  ERROR 1: Terminating translation prematurely after failed translation
  of layer Crewe_boundary_1m_grid (use -skipfailures to skip errors)

Yet if I try and run this for a different shapefile, i.e. a simple rectangular polygon, it works. I have also used exactly the same command on a 1000 m grid and that works. Also, if I do 'ogrinfo' then each of the input files reads as a polygon.

How do I fix this?


texts:
I'm new to GIS, trying to generate random points within a polygon in QGIS 2.18.7 through the GUI.

I go to:

Vector/Research Tools/Random points inside polygons (fixed)
This menu prompts me to select my input layer. But no layers are available from the pull down menu. How do I get my project layers to be available through this pull down menu?


texts: Using Google Maps images in Arcgis Online pop-up windows

texts:
I typically configure pop-ups in ArcGIS to show the image in the pop-up, by using dedicated URLs to the image which typically ends with .jpg or .png etc.

But, I am having trouble getting ArcGIS to recognize the URL to images hosted on Google Maps; a URL such as (https://www.google.com/maps/uv?hl=en&pb=!1s0x94b2354a240140d:0xfab1aeb455c4b700!2m22!2m2!1i80!2i80!3m1!2i20!16m16!1b1!2m2!1m1!1e1!2m2!1m1!1e3!2m2!1m1!1e5!2m2!1m1!1e4!2m2!1m1!1e6!3m1!7e115!4shttps://picasaweb.google.com/lh/).

I do not have any coding experience.

Any ideas?


texts: Calculating Centroids with Zonal Geometry in ArcGIS

texts:
I have a shapefile containing polygons. Each polygon is defined by its postal code. I want to calculate the centroid of each polygon using Zonal Geometry.

Sadly, for some polygons there are no centroids being calculated. It is the same when I use Zonal Statistics as Table. Some polygons seem to be ignored. There are no error messages.

It would help if someone could point me in the right direction because right now I do not even know what to fix. I have ArcGIS 10.1.

Thanks for any help!


texts: Is it possible to set a Tile Map Service only by copying the folder structure to a hosted server?

texts:
I used gdal2tiles to generate a directory with tiles of a bigger raster. Now, I'm able to see these local tiles on my browser, using the generated web pages (e.g. OpenLayers.html).

Thus, I copied the folder on a hosted server. The goal is to access the tiles with QGIS using the following an URL like http://mysite.com/gis/{z}/{x}/{y}.png. But I get no results.

Is it mandatory to configure the XML files (described here: Tile Map Service Specification)? Or am I missing something simpler?


texts:
I'm quite new to QGIS and trying to figure out how to analyze the following problem.

I have three layers:

settlements point layer with 253 entries
waste treatment facilities point layer with 24 entries
roads polyline layer

My problem is to find the distance between each settlement from layer (1) to each treatment facility in layer (2), however it has to be subjected to roads (layer 3), namely instead of getting aerial distance matrix I am looking for to the matrix that describes the actual driving distance.

I'm kindly asking your help with suggestion of how to perform this analysis with QGIS 1.7 (note that I don't know how to use PostGIS layers, if relevant)

what plugin can or should I use?
Should the polylines be connected to the points in both layers?


texts: Updated instructions for setting up TileCache on IIS 7.5?

texts:
This article seems great...but it is hard [for me] to follow when setting up for IIS 7.5:

http://viswaug.wordpress.com/2008/02/03/setting-up-tilecache-on-iis/

Any updated instructions available with latest versions of TileCache and IIS?


texts: How to apply Melt Area Detection Index to MODIS imagery?

texts:
I am working with some MODIS imagery for a project (I'm relatively new to remote sensing) using MOD09GA, and hoping to apply a melt area detection index (MADI) to images to calculate melt day length for different years. This was done in a paper by McCabe et al in 2010, where they present the equation for this. It simply divides band 1 by band 7, but I have no idea how to do this in ArcMap or ERDAS Imagine.

Looking for tips on how I would go about this. Do I need to create a composite bands image, or use raster calculator?


texts: QGIS3 for Mac: KyngChaos vs Homebrew vs MacPorts?

texts:
How have been your experiences between the three ways of installing QGIS3 on a Mac? I've used KyngChaos (which right now I think is the QGIS official distribution for Mac) and the Homebrew's one.

Some of what I have learned throughout the process:

If you want to keep the number of Python's installations on your computer low (if for example you are using Conda for Jupyter, GDAL for R, etc.), it's better to install everything using Homebrew. In my case, before installing everything using Homebrew, I ended up having official Python3 for KyngChaos QGIS3, Conda's Python3 for Jupyter, Homebrew's Python 3 & 2 for GDAL - GEOS - PROJ4. Now everything is Homebrew's Python2 & 3 (installing Jupyter using pip).
The main issue with QGIS3 from Homebrew is that it takes a really long time to compile and install (about 40 minutes), which makes it a little annoying to wait all that time when you upgrade your packages.

I still haven't installed QGIS3 using MacPorts, so I'm wondering how it is compared with the other two options.


texts: Creating map documents based on similar feature class selection

texts:
I wondering if something like this is possible. I would like to create a PDF output using a pre-established map template. The tricky part is I want to create maps in batch by adding (what will usually be 2)feature classes as layers to a data frame in the map document. These feature classes have similar names, the only difference is one of the names is preceded by "buffer_". The feature classes will then be symbolized. I am very new to arcpy mapping and seems like this may be possible however, I am stumped as to how to do this. The limited code I have so far (which doesn't run properly) looks as follows:

import arcpy, os, fnmatch

from arcpy import env

mapTemplate = arcpy.mapping.MapDocument(r"template location")
try:
df = arcpy.mapping.ListDataFrames(mapTemplate)[1]
print df
arcpy.env.workspace = r"file geodatabase"
fc = arcpy.ListFeatureClasses()
for fcs in fc:
    name = fcs.getValue().replace("Buffer_","")
    for i in os.listdir("file geodatabase"):
        if fnmatch.fnmatch(i, name):
            arcpy.MakeFeatureLayer_management(fcs, "fc_layer")
            layer = arcpy.mapping.Layer("fc_layer")
            arcpy.mapping.AddLayer(df, layer + fcs, "AUTO_ARRANGE")
            mapTemplate.save()

except Exception, e:
import traceback, sys
tb = sys.exc_info()[2]
print "Line %i" % tb.tb_lineno
print e.message

I know this is no where close but I am having trouble even using the .getValue statement for the feature classes. Thanks for any suggestions or help you can provide.


texts: Overlap map and frame in QGIS print composer

texts:
Often in old maps you see where the cartographer has "cheated" by allowing part of the map to extend over the frame. You can see an example at the top (or the bottom of the inset) of this map:

Is there any way to do this in QGIS print composer?


texts: PostGIS. Find intersection of a polygon from 2 points

texts:
I have 2 polygons and 2 points stored in postgres tables and Im trying to find the front and back of one of the the polygons (green dots) based on a line from the 2 red dots.
I'm trying to find the green dot that is intersected by the line between the 2 red dots.
select
    st_astext(ST_Intersection(p.poly_geog,st_makeline(st_makepoint( c1.longitude, c1.latitude),
                            st_makepoint( c2.longitude,c2.latitude))))
from
    polygonfence p
join index i on
    p.index_id = i.id
join coordinate c1 on
    i.start_coordinate_id = c1.id
join coordinate c2 on
    i.centre_coordinate_id = c2.id


texts: How to find Chainage of points using QGIS

texts: I have two types of Road shape files. Road (Polyline) and important points (waypoints) along the road. Both of these shape files are prepared using GPS. I want to find the chainage of the points. Is it possible using QGIS ?


texts: When trying to join Excel datafile to shapefile in ARC GIS 10.2.2 I keep getting the following message:  "Failed to connect to database.  An underlying database error occurred.  Class not registerd."  Any ideas how to resolve?


texts: True-North Bearing between two GPS coordinates

texts:
I am trying to find the correct equation to find the bearing off of true-north between two GPS coordinates. I have seen a lot of online calculators and a lot of formulas on this GIS stack exchange, but many of them give me different results! What is the most correct method?

It seems the most popular one is this method: Calculate bearing between two decimal GPS coordinates (Arduino / C++)

Another method is this one: Calculate bearing between two decimal GPS coordinates , but, I am getting completely different answers when using this method (off by like tens of degrees), and I can't even get a matching result on the example from that post: "I get 250.20613449 as expected."

On all of my attempts, I am consistently using the degrees notation for lat and lon (I'm not using the hours, minutes, second notation). Also, on all of my attempts, I swap the start/stop position of both coordinates and recalculate just to make sure I'm not solving relatively to the wrong start position.

Online calculators seem to give me a variance as well. Sometimes they match the first method, sometimes they don't.

Can someone help nip this in the bud, once and for all?

Here's an example prob we can focus on:

Start: 42.766698 -75.746404
Stop: 42.761405 -75.757653

If I plot these points on a map and guesstimate the bearing by eye:


texts: QGIS 3.10 projection error using BNG27700

texts:
When opening a workspace with CRS BNG 27700 I am now getting an error that the following transformation is missing - OSTN15_NTv2_OSGBtoETRS

If I use the next transformation in the list it puts my WMS layers with CRS 3857 slightly off where it should be, up to a few metres but enough to know it is not correct.

I have tried adding the Suggested file OSTN15_NTv2_OSGBtoETRS into the folder named 'proj' it suggests. Note this folder was missing and was created by myself at the suggested location.

When I do this I get the following error when loading in my WMS in CRS 3857 ie. Google roads for example

as you can see this is far from ideal as it only shows this single strip of the layer I am wanting to view. Has anyone else got this same issue when using BNG 27700 in QGIS 3.10, running on Windows 10?


texts: My georeferencer plugin window doesn't open

texts: Ever since I downloaded QGIS 2.2 I can't get the georeferencer plugin to open. I've removed all older versions of QGIS & reinstalled & have just installed 2.8.6 but am still having the same problem. I select Raster -> Georeferencer-> Georeferencer but no window opens. Has anyone else had this problem & is there a solution for it?


texts: Creating buffers that automatically update when new feature added using ArcGIS Online

texts:
Is it possible to create a new feature in ArcGIS Online and the buffer automatically update?
I need to create a point in ArcGIS Online and produce the point together with its buffer.


texts: How to remove long linear strips from polygons?

texts:
I have a polygon layer representing impervious land cover.  This was taken from a land cover raster dataset.  The impervious data represents all impervious surface that is not buildings or roads.  I am trying to identify parking lots within the impervious features, but I am getting all sorts of omission and commission errors.
What I am finding is that I will identify parking lot features in the impervious layer, but they will be connected to other impervious (non parking lots) by long narrow strips of impervious (could be sidewalks or unofficial roads).  I'm thinking if I can somehow remove these strips from my impervious features, I could get more accurate parking lots.

As you can see in the image, the selected polygon (in cyan) covers the parking lot at an airport, but is connect by impervious to the runway.  I just want the parking lot.
Can anyone recommend a way to remove the strips from the polygons?  I can't think of way using Vector tools in ArcGIS 10.0.  Is there a way if I convert to raster?


texts: (pyqgis) Is there a way to disable the progress bar in the processing tool

texts:
When using processing tool in pyqgis and run, it always creates a progress bar. Is there a way to disable it?

processing.runalg('qgis:fixeddistancebuffer', _layer, distance_tolerance, 10, False, None)

I want to disable it because I get this error when I run my script in travis:

  wrapped C/C++ object of type QProgressBar has been deleted See log for
  more details

  RuntimeError: wrapped C/C++ object of type QProgressBar has been
  deleted


texts:
I'm trying to automate my production of maps with GRASS.

What I want to do can be described in fairly generic terms: I have two layers (one vector, one raster), which will be identical in every map. I have one layer (a raster heat-map) that will come from a different file for every map. I can write a shell or Python script (Mac OS X) to loop some GRASS command over all the heat-map files in a directory, but I am not sufficiently competent with GRASS to have figured out how to handle that end of the process.

Could someone give me an example of GRASS command-line code that will take a few files as inputs and produce a single image file as output?


texts: How can I specify the ranges in the gdaldem's color configuration file?

texts:
I want to render an IMG file according to a specific renderer, and export to a jpeg File. I have asked a question about this, and got a satisfactory answer.

The issue now is that I need to write a color configuration file which matches the symbology that is required. My Symbology is supposed to be like this:

Based on this answer, I have found that you can write the colors config file somewhat like this:

40  140 0   3
35  204 9   9
30  255 67  11
25  255 165 0
20  255 246 0
15  54  254 0
10  0   255 209
5   1   127 255
0   4   0   213
-4  133 1   243

My problem is that I don't understand if the given values for the pixels are the lower bound, or the upper bound. I don't know how to specify the ranges according to the symbology that I need.

Where can I Find documentation for Gdaldem's color configuration file? Or alternatively, how can I specify the ranges in the configuration file?


texts: Convert specific RGB (GeoTIFF) color to transparent

texts:
I have an RGB GeoTIFF file and I want to convert a specific RGB colour to transparent. I would like to use GDAL for doing this, but external Python scripts are also an option (though I haven't learned the art of Python scripting yet).

My thoughts are:

Add a 4th band to my GeoTIFF file and fill with 255.
Insert 0 into the 4th band where RGB colour == target colour

The first part I think I can solve using

gdalwarp -dstalpha source.tif target.tif

The second though I don't know how to solve. matt wilkie suggest using a Python script he wrote, but I think that is done separately for each band and not for a specific colour.

How can I assign a value to the 4th band of my TIFF when the RGB colour is equal to a specific RGB colour?


texts: Performing raster calculations along line feature or within polygon using QGIS?

texts:
I need to add a value to a dem file along a polyline feature or within a polygon.  How can I add 100 to the values of the raster cells intersected by my line feature or within my polygon feature.

I would like to know how to do this in QGIS, if possible.


texts: Selecting certain area of polygon according to size in ArcMap?

texts:
I have a set of polygons with a certain size (3500 sq meters each for examples). I only need 3000 sq meters of every polygon.
Is there an automatic or semi-automatic way to create sub-polygons out of those polygons with the desired size ?
Also, if indeed it's possible, is it possible to select from what side of the polygon you want to create the sub-polygon ?


texts: Squaring up photo into Google Earth using spatial adjustment

texts:
I have a photo such as:

that I want to add into a Google  Earth  map. I know where the image must be: North: 44°16'48.00"N, South: 35°15'0.00"N, East: 4°42'0.00"E, West: 10°1'12.00"O and I have obtained this:

I think I can adjust the photo into the Google Earth  photos with software such as QGIS but how do I do that?


texts:
OS=Ubuntu 20.04
I have imported OSM data into a postgis database.
I'm now trying to render tiles with generate_tiles.py which is in the repo mapnik-stylesheets. However, this doesn't work as mapnik tells me that all fonts in my xml file are not installed, which is definitely not true (renderd uses mapnik as backend and works with said xml file).
Where does mapnik look for fonts? What can I adjust to get this working?


texts:
I have my own OSM file, which is created by some application which I written.

You can see some of nodes in my OSM file bottom of the post.

I want to fill +/-0.0001 latitude, +/-0.0001 longitude around my point.

I've tried these to style points, however, sadly, height and width attiributes seem to not working :/

What may be the other solution(s) to fill +/-0.0001 lat-long around my point?

layer-amenity-points.xml.inc

<Style name="points">
    <Rule>
      &maxscale_zoom18;
      <Filter>[amenity]='mycustomamenity'</Filter>
      <PointSymbolizer height="52" width="52" file="&symbols;/overlay_img.png" placement="interior"/>
    </Rule>
    <Rule>
      &maxscale_zoom17;
      <Filter>[amenity]='mycustomamenity'</Filter>
      <PointSymbolizer height="26" width="26" file="&symbols;/overlay_img.png" placement="interior"/>
    </Rule>
</Style>

<Layer name="amenity-points" status="on" srs="&osm2pgsql_projection;">
    <StyleName>points</StyleName>
    <Datasource>
      <Parameter name="table">
      (select amenity
      from &prefix;_point
      where amenity is not null
      ) as points</Parameter>
      &datasource-settings;
    </Datasource>
</Layer>

custom.osm

<node id="10000001" lat="0.476527778159" lon="0.999861111111" visible="true">
    <tag k="amenity" v="mycustomamenity"/>
    <tag k="myvalue" v="573"/>
</node>
<node id="10000002" lat="0.476805555937" lon="0.999861111111" visible="true">
    <tag k="amenity" v="mycustomamenity"/>
    <tag k="myvalue" v="576"/>
</node>
<node id="10000003" lat="0.477083333715" lon="0.999861111111" visible="true">
    <tag k="amenity" v="mycustomamenity"/>
    <tag k="myvalue" v="538"/>
</node>

Update

With MarkerSymbolizer in the answer below, I couldn't find how to draw squares. There is just circles can be created by MarkersSymbolizer. I want to draw squares.


texts: How to limit the number of items to enter metadata on when importing from GeoServer into GeoNetwork

texts: We are a library looking at using GeoServer and GeoNetwork as the basis of a GIS repository/catalog. When we enter layers in GeoServer and then harvest using GeoNetwork, we end up with a large number of layers for which we have to enter metadata in GeoNetwork. University of Virgina have a tool that they use to group layers, and then publish the record for the group to GeoNetwork, but I was wondering if other people had come up with solutions that they would be willing to share.


texts: subtract two DEMs and extract values that are greater than 1 and set zeros to null or no data

texts:
I have two DEMs. my plan is to subtract 1 DEM from the other and extract only the values that are greater than 1 and set values that are less than 1 to no data or null

I can do this in 3 steps but I would like this to be done in 1 query in the raster calculator

1. dem2 - dem1 #output = subtract_DEM
2. subtract_DEM >= 1 #output = subtract_DEM_grt1
3. SetNull(0, subtract_DEM_grt1)

how can I compute this into one expression using the arcgis raster calculator


texts: Merge landsat8 bands and convert to tiles using gdal

texts:
I have downloaded some landsat8 TIFF files from Amazon and there are 11 bands. Then I use for BAND in {4,3,2}; do
  gdalwarp -t_srs EPSG:3857 LC08_L1TP_114027_20171116_20171116_01_RT_B$BAND.TIF farm/$BAND-projected.tif;
done and then
convert -combine {4,3,2}-projected.tif RGB.tif.
The color doesn't seem right. I want to convert the RGB.tif to tiles to put it on a leaftlet map, so I use gdal2tiles.py -t tile -e -a 0 -z 0-12 RGB.tif img but the error is shown:gdal2tiles.py: error: There is no georeference - neither affine transformation (worldfile) nor GCPs. You can generate only 'raster' profile tiles.
What should I do now to merge the TIFFs and put it on a map as tiles?


texts: Why is the curl syntax unable to create datastore?

texts:
I am implementing the following code to upload an existing shapefile as per the documentation,

curl -u admi PUT 'Content-type: application/zip' \
      --data-binary @/var/www/geo/shapefile/zip/Poa_horridula.zip \
      http://localhost:8080/geoserver/rest/workspaces/acme/datastores/Poa_horridula/Poa_horridula.zip

But I am getting following error,

  curl: (6) Couldn't resolve host 'Content-type'

    No such datastore: acme,Poa_horridula

The document says that the curl code will create the datastore then why am I getting this error?

UPDATE: I changed the code to upload a zip file but I am still not able to upload any shapefile.

Where am I going wrong?

UPDATE 2: Ok I am almost there

So I changed the command to,

curl -u admin:geoserver -v -XPUT -H 'Content-type: application/zip' --data-binary @/var/www/geo/shapefile/zip/Poa_horridula.zip http://localhost:8080/geoserver/rest/workspaces/acme/datastores/Poa/file.shp

The only problem is that I am able to add the shapefile but the layer is not being published. Is that normal? Do I have to manually publish it?


texts: Add columns if layer intersects more than one other layer in PostGIS

texts:
I have been spatially joining land use zones to parcels successfully with the following query:

CREATE TABLE new_parcel AS
SELECT DISTINCT ON(t.geom) t.*, m.basezone, m.overlay
FROM parcel AS t
LEFT JOIN zoning AS m
  ON st_intersects(t.geom, m.geom)
  AND st_area(st_intersection(t.geom, m.geom)) / st_area(t.geom) > 0.5

Currently, more than 50% of the parcel's area must be within a single zone for a join to occur.

However, it turns out that some of the parcels overlap 3 or more zones and are not getting joined to the zoning table because less than 50% of their area overlaps.  For example, the parcels of dark purple color overlap more than two zones:

How would one modify the query to add extra columns in the resulting table which would include data from all zones overlapped?

The existing code produces a table that has

basezone | overlay

How could one produce

basezone_2|overlay_2|basezone_3|overlay_3|...

and so on for all zones overlapped?


texts:
I have a problem with import data with fme and arcsde.
I have a loop of import data in the same arcsde table.

I have a problem with index :

  Error while executing stream insert operation for table
  'VLORA1A.TEST'.:SDE_ERROR_CODE(-51) Underlying DBMS error.
  SDE_EXT_ERROR(29861) SDE_ERROR_MSG1(ORA-29861: l'index de domaine est
  à l'état LOADING/FAILED/UNUSABLE ) SDE_EXT_ERROR_MSG2()

I drop the index and i register the table with sdelayer.

I think FME alter the spatial index.

Have you a solution ?

Thanks


texts: Match a list of shapefile names to list of raster names then do something

texts:
I have a list of shapefiles and a list of rasters that I would like to match the names then do some data extraction for later analysis. The beginning of each raster and shapefile begins with "yyyy-mm-dd_UNL_Wheat_NIN" and I would like to use this as the match criteria. I have the code below to show what I would like to do. Right now I only have it able to iterate over a list of orthos with one shapefile. How would I make it so there is a conditional loop prior to do something loop?
shape_path = glob.glob('/Volumes/ML - External/2017-2018 season/Shapefiles/'+'*.shp')

/Volumes/ML - External/2017-2018 season/Shapefiles/2018-05-05_UNL_Wheat_NIN.shp
/Volumes/ML - External/2017-2018 season/Shapefiles/2018-05-09_UNL_Wheat_NIN.shp

raster_path = glob.glob('/Volumes/ML - External/2017-2018 season/Orthos/'+'*.tif')

/Volumes/ML - External/2017-2018 season/Orthos/2018-05-05_UNL_Wheat_NIN_Flight 1.tif
/Volumes/ML - External/2017-2018 season/Orthos/2018-05-05_UNL_Wheat_NIN_Flight 2.tif
/Volumes/ML - External/2017-2018 season/Orthos/2018-05-09_UNL_Wheat_NIN_Flight 1.tif
/Volumes/ML - External/2017-2018 season/Orthos/2018-05-09_UNL_Wheat_NIN_Flight 2.tif

Here is what I would like to loop through and extract to an array for the shapefile that matches the raster.
# Iterate through ortho names
for ortho in ortho_path:
    # Print name as a check
    print(ortho)

    with rasterio.open(ortho) as img:
        # Code to automatically find date in yyyy-mm-dd in file name
        match = re.search(r'\d{4}-\d{2}-\d{2}', img.name)
        img_date = datetime.strptime(match.group(), '%Y-%m-%d').strftime('%Y-%m-%d')

        # assign location name
        #img_loc = img.name.split("/")
        img_loc = ""
        # Get features and barcode name from shape file
        for feature in shape:
            shape_cord = [feature["geometry"]]
            barcode = feature["properties"]['Barcode']
            out_image, out_transform = mask(img, shape_cord, crop=True)
            GLI = (2*np.mean(out_image[1,:,:]) - np.mean(out_image[0,:,:]) - np.mean(out_image[2,:,:]))/(
                2*np.mean(out_image[1,:,:]) + np.mean(out_image[0,:,:]) + np.mean(out_image[2,:,:]))

            # Appending values to array
            flightdate.append(img_date)
            flightloc.append(img_loc)
            brcd.append(barcode)
            gli.append(GLI)


texts: How to combine use gradient symbols for layer properties excluding 0?

texts:
I want to visualize my data by using graduated symbols. Now, there are some objects with value "0", they should not be printed on the map (which I solved by using rule-based symbols). But, anyway, "0" is included when automatically generating the (equal interval) classes. This is important because my data are values of more or less 1800 to 2000.
Also three different sizes are applied, but I don't think that will be important for this question.

How do I exclude "0", when assigning gradual colours?


texts: What is the status of any ongoing efforts that will allow QGIS to display .las files directly in the main 2D window?

texts:
After some research on QGIS's ability to display .las files, I came to the conclusion that it is not possible: currently any attempt at this requires converting .las to .shp as described here and here.

I am wondering if there is an ongoing effort to make direct .las viewing possible in the main window of QGIS (not a separate 3D window), and more broadly, what the main challenge is from a development perspective? (I'm no expert programmer, so I am looking for a digested explanation if anyone can provide one.)

Are there any similar efforts happening in other opensource GIS software that might be ahead of those of QGIS?


texts: How to get distance between one point and all other points in a table?

texts: I have a table in Postgres with lat and long columns. Can I use the ST_Distance_Sphere or ST_Distance_Spheroid methods to get the distance between one point and all other points in the table? Do I first have to convert my existing lat/long to the Point datatype; if so is there any method to do it?


texts: When installing ArcGIS Server 10.1, I cannot get the System/Publish service to start

texts: When I attempt to start it, it displays "System/Publishing Tools service is starting..." but never ends, time's out, or gives an error message.  Any ideas what is going on here?  I have already installed the license manager, ArcGIS Desktop and the license codes.  ArcGIS Server seems installed correctly but just won't start the publishing service.  I'm trying to develop with the ArcGIS Javascript API which is why I'm installing ArcGIS Server on my PC.


texts: How to find area of Digital Elevation Model above/below a reference plane?

texts: I am using ArcGIS10. I need to find the out the area which is below 5000 and above 5000 from the Digital Elevation Model which is available to me. How can I find out this ? My DEM ranges from 8848 to 309 m with the total area of 4100 sq.km.


texts: Erasing buffer within another buffer without erasing overlapping buffers using ArcGIS for Desktop?

texts:
I'm not an ArcGis pro, so I need help on this one, and I have tried searching for answers.

I am using ArcMap 10.2.2. I don't know how to use python.

I have 6km and 12km buffers around the sites that I'm interested in. I want to be able to clip the 6km buffer from the 12km buffer so that I get a new layer that has buffers that are only from 6km to 12km. I have tried using the clip tool, but instead of cutting out the 6km buffer from the 12km buffer, it seems to only be able to cut away the 12km area surrounding the 6km buffer (see picture). The buffers are all in the same layer by the way.

I then tried using the the erase tool. I erased the 6km buffer from the 12km buffer, which works, but then anywhere there was buffer overlap, the erase tool has cut out areas of overlap and I want the new (6-12km) buffers to remain complete doughnuts (see picture).

How can I clip or erase the 6km buffer from the 12km buffer without cutting out overlap too? Picture:


texts: How did sailors navigate before the Mercator projection (1536)?

texts: The Mercator projection (1536) must be great for nautical navigation. It allows one to find a bearing and keeping this constant one should reach the target (if not in the shortest line). But how did people navigate before this? As an example Magellan circumnavigated the globe ca. 1520.


texts:
I create a layer package from desktop 10 and the double click to decompress it on my machine (to see what the customer is going to see).
I get the initial folder (same as my lpk file) and then another one duplicated inside.  c:\documents and settings\users\myusername\My Documents\ArcGIS\Packages\mylpkfile2-17-11\mylpkfile2-17-11
There are data files in both folder structures (with shape files in both).
There is a lyr file in the root folder that points to the first dataset.
I see a lyr for v10 in the v10 folder.
However this lyr uses the same data as the other lyr.

It looks to me that there is no reason for this data structure.
I don't know if it is just my machine or ...   Does anyone else see this also?

The above occurs if I double click on the lpk. If I run the extract package tool from within arcmap it extracts only one folder structure to the location I instruct it to.


texts: Where to start learning how to make maps with QGIS?

texts:
I have no skills in GIS, "map-making," or any other map program. Been a county planner for 25 yrs, where there was always a "drafting" staff.  I'm in another agency now and need to be able to produce basic and attractive site plans and area maps.  Is QGIS the right tool? I saw on a planners' website that it's free to download, along with GIMP and SketchUP.

I need some genuine advice on this, not something impossibly technical for my situation, nor a smug brush-off.


texts: Can PostGIS return a list of filenames to rasters that intersects with a vector polygon?

texts:
I have 390 invisibility depth viewshed rasters, each with a size of 687 MB. The rasters contain the maximum value a building could attain at any given location, before it's visible from any of the 390 viewshed points.

Would it be possible to get PostGIS to query all the rasters that intersects with a vector polygon (could be converted to raster if easier) and return a list of filenames to those rasters that the vector polygon (a planned building with height in a column called "Z") intersects with?

I would like to be able to take a planned building and without to much manual work get a list/summary over which of the 390 rasters it intersects with.

If this is possible, how would you recommend me to set it up to be able to minimize manual steps?


texts: Finding the sequence/rank of lines (polylines) in QGIS

texts:
I have line shapefile with different bus routes (each route consist of multiple line), in which I want to sequentially arrange the lines for each bus route based on the direction of route.
I tried with QChainage, but couldn't achieve it.
Is it possible in QGIS?


texts: SyntaxError with cursor and Feature Selection in ArcPy?

texts:
I tried to select each feature in the layer "temp" and use the selected feature to select features in the other layer "Railway" based on location and then do some calculations with the "Railway" layer using the values in the "temp" layer.  Here is the code:

import arcpy
fc = "C:/Users/xkuai2/Downloads/temp/temp.shp"
cursor = arcpy.SearchCursor(fc)
for row in cursor:
    count = 0
    Vol = row.getValue("Volume")
    arcpy.SelectLayerByAttribute_management("temp", "NEW_SELECTION", "FID" = count)
    arcpy.SelectLayerByLocation_management("Railway", "WITHIN", "temp", , "NEW_SELECTION")
    arcpy.CalculateField_management("Railway", "Volume", "[Volume] + Vol", "VB")

But I came across with the error message:

Parsing error SyntaxError: invalid syntax (line 5)

I checked the code and have no idea what goes wrong.


texts: Django PolygonField: how to display on map?

texts:
Extremely beginner question.

I have a GeoDjango model with a MultiPolygonField. How do I retrieve the value of this field and display it as a polygon on a map?

class county(models.Model):
    cd = models.CharField(max_length=10)
    geom = models.MultiPolygonField(srid=27700)
    objects = models.GeoManager()

def lsoa(request, code):
    county = get_object_or_404(county.objects, cd=code)
    return render_to_response('county.html', { 'county': county }, context_instance = RequestContext(request))

I know GeoJson is involved somewhere along the line, but I'm struggling with to get it out of Django, and with how to display it.


texts: How to cluster Points with ST_ClusterDBSCAN taking into account their type stored in a separate column?

texts:
I was looking to use ST_ClusterDBSCAN on a set of points that also have a type. So I want the points to be clustered by type, so that points having different types won't end up in the same cluster. Since ST_ClusterDBSCAN is windowed function I thought I could simply PARTITION BY type column and that will do. Something like:

ST_ClusterDBSCAN(point, eps := 20, minpoints := 3) OVER(PARTITION BY type) AS cid,

However I still end up having same cid assigned to points of different type.

Obviously there's something I misunderstand either about clustering or windowed functions in general?

ST_ClusterDBSCAN assigns cids differently with and without of PARTITION BY type so partitioning definitely affects the outcome somehow, I'm just not sure how exactly. Obviously not the way one would expect.


texts: Digitizing complex archaeological scenarios in various levels?

texts:
I don't know if this issue has already been answered, but I can't get my head around it. I used to do this in Illustrator, and there it was simply a matter of clicking "bring forward/bring back" etc.

We are about to start digitizing a lot of orthophotos of archaeological scenes generated in Agisoft Photoscan Pro. The excavations contains A LOT of wood and other artefacts. Various orthophotos covers several vertical layers (as we have dug down, each layers has been documented).

I have been digitizing two overlapping orthophotos, but the main problem is, that many of the artefacts or pieces of wood, overlap each other.
I would really like to avoid, having to create layers for EACH artefact, and thus being able to contain fx all "Mandibles" in its corresponding Shapefile/layer. So, that the mandibles with the most positive Z-values, wouldnt be overlapped by wood with lower Z-values, simply because the layer for wood is above mandibles in the layers list.

Would it make a difference, if each polygon was assigned a Z-value (possibly extracted from a DEM)?


texts:
I need to open Google Maps at Latitude: -30.0330600 Longitude: -51.2300000.

We report that the object map or call the API?

<!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>

 <!-- Chamada a Api Google Maps -->
 <script src="http://maps.google.com/maps/api/js?sensor=false&v=3.2"></script>

 <script type='text/javascript'>

 //Variável map
 var map;

 //Função a ser chamada para criar o mapa na leitura da página
 function init() {

 // A classe Map espera dois argumentos "map_element" que o ID do elemento HTLML geralmente uma DIV e opções de mapa {key:value}
 map = new OpenLayers.Map('map_element',{
    maxExtent: new OpenLayers.Bounds(
        -128 * 156543.0339,
        -128 * 156543.0339,
        128 * 156543.0339,
        128 * 156543.0339),
    maxResolution: 156543.0339,
    units: 'm',
    projection: new OpenLayers.Projection('EPSG:900913'),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
 });

 //Função de inicialização - google.maps.MapTypeId.ROADMAP (the default)
 var google_streets = new OpenLayers.Layer.Google(
    "Google Streets",
    {numZoomLevels: 20}
    );

//Adiciona as camadas ao mapa
map.addLayers([google_streets]);

 //Camada de controle que vai mostrar as camadas no mapa
 map.addControl(new OpenLayers.Control.LayerSwitcher({}));

 center = new google.maps.LatLong(-30.0330600,-51.230000);

 if(!map.getCenter(center)){
 map.zoomToMaxExtent();

 }
 }

 </script>
</head>

<body onload='init();'> <!-- Chama a função js init() -->

<!--Elemento HTML onde o mapa é exibido -->
 <div id='map_element' style='width: 50px; height: 500px;'>
 </div>

</body>
</html>


texts: How to transform a selected multipart feature into singlepart features while editing in QGIS?

texts:
There are several questions and respective answers for how to make a multipart features Layer into a  singlepart features one using QGIS:

How to break a single polygon into multiple polys?
How to split non-contiguous feature in QGIS?

But all the solutions I saw, create entirely new layer and would transform ALL the multipart geometries into single ones.

So the questions is, how could one, during an edit session, selectively transform a multipart feature into a singlepart feature, without the need to create new layers?

It would be the opposite to the "merge selected features" tool, and similar to explode in ArcGIS.

Is there any plugin doing this? How could this be done using the python console?


texts: How to import xyz data in QGIS 3.0 Processing Modeler for creating points and raster?

texts:
I want use dsm data in QGIS 3.0 Processing Modeler. This data are in xyz-format.
I want to convert xyz data into points to create a raster. Later I want also apply this modle for the a folder with several xyz data in a folder.
I have tried tool “create points form table”.  Doesn't worked
Outside the from Processing Modeler is possible to import xyz data  with  “delimited Text”


texts: ArcGISDynamicMapServiceLayer Zoom + LocateButton

texts:
Have few questions about DynamicMap options.

1.Can i implement LocateButton?

2.Can i set the starting Zoom on the map?

This is my Code that is use only with LocateButton but wont work with him, if LocateButton is taken off the map work fine.

<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
<title>Create Map and add a dynamic layer</title>
<link rel="stylesheet" href="https://js.arcgis.com/3.16/esri/css/esri.css"/>
<link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.6/css/bootstrap.min.css">
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.12.0/jquery.min.js"></script>
<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.6/js/bootstrap.min.js"></script>
<style>
  html, body, #mapDiv{
  height: 100%;
  width: 100%;
  margin: 0;
  padding: 0;
  }

  #search {
  display: block;
  position: absolute;
  z-index: 2;
  top: 20px;
  left: 84px;
}

#LocateButton {
  position: absolute;
  top: 95px;
  left: 20px;
  z-index: 50;
}
</style>
<script src="https://js.arcgis.com/3.16/"></script>
<script>
  var map;

  require([
    "esri/map",
    "esri/layers/ArcGISDynamicMapServiceLayer",
    "esri/dijit/LocateButton",
    "esri/layers/ImageParameters"
  ], function (
    Map, LocateButton, ArcGISDynamicMapServiceLayer, ImageParameters) {

    map = new Map("mapDiv", {
      sliderOrientation : "horizontal",
    });

    geoLocate = new LocateButton({
    map: map
    }, "LocateButton");
    geoLocate.startup();

    var imageParameters = new ImageParameters();
    imageParameters.format = "jpeg"; //set the image type to PNG24, note default is PNG8.

    //Takes a URL to a non cached map service.
    var dynamicMapServiceLayer = new ArcGISDynamicMapServiceLayer("http://gis.katastar.gov.mk/arcgis/rest/services/Raster/Ortofoto_MK_color/MapServer", {
      //"opacity" : 0.5,
      "imageParameters" : imageParameters,
    });

    map.addLayer(dynamicMapServiceLayer);
  });

  require([
   "esri/dijit/Search",
   "esri/layers/FeatureLayer",
   "esri/InfoTemplate",
   "dojo/domReady!"
  ], function (Search, FeatureLayer,InfoTemplate) {

        var search = new Search({
        enableButtonMode: true, //this enables the search widget to display as a single button
        enableLabel: false,
        enableInfoWindow: true,
        showInfoWindowOnSelect: false,
        map: map
     }, "search");

     //var sources = search.get("sources");
    var sources = [];
     //Push the sources used to search, by default the ArcGIS Online World geocoder is included. In addition there is a feature layer of US congressional districts. The districts search is set up to find the "DISTRICTID". Also, a feature layer of senator information is set up to find based on the senator name.

     sources.push({
        featureLayer: new FeatureLayer("http://62.162.76.153/arcgis/rest/services/Kiosci_MKD/MapServer/0"),
        searchFields: ["HR_TITLE"],
        displayField: "HR_TITLE",
        exactMatch: false,
        outFields: ["HR_TITLE", "OBJECTID", "HR_DESCRIPTION"],
        name: "Znamenitosti",
        placeholder: "Enter Value",
        maxResults: 50,
        maxSuggestions: 10,

        //Create an InfoTemplate and include three fields
        infoTemplate: new InfoTemplate("Skopje City", "Point Name: ${HR_TITLE}</br>Description: ${HR_DESCRIPTION}"),
        enableSuggestions: true,
        minCharacters: 0
     });

     sources.push({
        featureLayer: new FeatureLayer("MY SERVICE"),
        searchFields: ["HR_TITLE"],
        displayField: "HR_TITLE",
        exactMatch: false,
        name: "Smestuvanje",
        outFields: ["*"],
        placeholder: "Enter Value",
        maxResults: 50,
        maxSuggestions: 20,

        //Create an InfoTemplate

   //Set the sources above to the search widget
     search.set("sources", sources);

     search.startup();

});

</script>


texts: Finding overlapping attributes in two different layers using QGIS?

texts:
I have two layers - one with country borders, and one with river basin borders - and I need to find a way to generate a list of basins that occur in each country.

Is there a way to find overlapping attributes (basins and countries) from two different layers using QGIS?


texts: Modelling of rent prices - which interpolation method to use?

texts:
I have national dataset of ~1,4 million households. There I have information about rent, size (number of rooms and m2) and some additional characteristics of each household.

I'd like to use this data to create surface of rent prices for the whole country and use this information as a proxy for estimation of values of remaining ~1.5 million households that are owned or do not have rent information.

Couple of questions here:

Is such approach appropriate for the this kind of problem at all?

Which method of interpolation would be most suitable to use here?

Also, would it be possible to take information about, for instance the size of the household into account?

I'm on ArcGIS 9.3 with ArcInfo license.


texts: Find Polygon in "Layer A" containing Polygon in "Layer B"

texts:
I haven't found yet someone with the same problem, so I'm posting a new question. The topic seems quite standard, so I'm sure there's a way to do it in QGis.

I have 2 layers:

Region Layer
Postcode Layer

both containing (multi-)polygon features.

I would like to cycle through all the Postcodes, find which is the Region that contains that Postcode, and eventually write the ID of the region inside the Postcode.

So far I've tried to write something in the Python Console, which I can share here:

from qgis.utils import iface
from PyQt4.QtCore import QVariant
_ID_FIELD = 'RS'
_NAME_FIELD = 'GEN'
_PLZ_FIELD = 'plz'
_NEW_ID_FIELD = 'PARENT_ID'
_NEW_NAME_FIELD = 'PARENT_NAME'
Land_layer = iface.mapCanvas().layer(1)
PLZ_layer = iface.mapCanvas().layer(0)
# Create 2 new fields in the layer that will hold the list of neighbors and sum
# of the chosen field.

#PLZ_layer.startEditing()
#PLZ_layer.dataProvider().addAttributes(
#        [QgsField(_NEW_ID_FIELD, QVariant.String),
#         QgsField(_NEW_NAME_FIELD, QVariant.Int)])
#PLZ_layer.updateFields()

## Create a dictionary of all features
PLZ_feature_dict = {f.id(): f for f in PLZ_layer.getFeatures()}
Land_feature_dict = {f.id(): f for f in Land_layer.getFeatures()}
# Build a spatial index
index = QgsSpatialIndex()
print 'PLZ Layer: %s' % PLZ_layer.name()
print 'Landkreis Layer: %s' % Land_layer.name()

for f in Land_feature_dict.values():
    index.insertFeature(f)
    print 'Inserting %s in Landkreis index' % f[_NAME_FIELD]

for f in PLZ_feature_dict.values():
    print 'Working on %s' % f[_PLZ_FIELD]
    geom = f.geometry()
    intersect_ids = index.intersects(geom.boundingBox())
    for id in intersect_ids:
        tempf = Land_feature_dict[id]
        print 'got parent Landkreis %s' % tempf[_NAME_FIELD]

Problem is, I'm getting the following error at line

intersect_ids = index.intersects(geom.boundingBox())

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "//VAIMUCO9/G114305$/Anwendungsdaten/Desktop/PLZ/ContainerFinder.py", line 34, in  <module>
    intersect_ids = index.intersects(geom.boundingBox())
AttributeError: 'NoneType' object has no attribute 'boundingBox'

Moreover, I'm not really sure that I should use these functions, and even this approach.
Can anybody please suggest me some more clever way?


texts: QGIS WFS Error - Could not commit changes to layer

texts:
I have a PostGIS database being served by Geoserver. I can load the layers I want into QGIS no worries; however whenever I try to edit the layer I get this error message:

Could not commit changes to layer X
Errors: ERROR: 1 geometries not changed.

Provider errors:
      empty response

Does QGIS-2.4 actually support wfs-t capabilities at the moment?

EDIT

QGIS

Network request http://*******.com:**/geoserver/****/wfs?SERVICE=WFS&VERSION=1.0.0&SRSNAME=EPSG:3308 timed out

Geoserver

2014-09-29 09:43:06,297 INFO [geoserver.wfs] -
Request: getServiceInfo
2014-09-29 09:43:06,330 INFO [geoserver.wfs] -
Request: getServiceInfo
2014-09-29 09:43:06,330 INFO [geoserver.wfs] -
Request: describeFeatureType
    service = WFS
    version = 1.0.0
    baseUrl = http://***.com:**/geoserver/
    typeName[0] = {*****}*****
    outputFormat = XMLSCHEMA

POST Request

<Transaction xmlns="http://www.opengis.net/wfs"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0"
xmlns:NCTPROP="NCTPROP" service="WFS" xsi:schemaLocation="NCTPROP
http://****.com:****/geoserver/****/wfs?
SERVICE=WFS&amp;VERSION=1.0.0&amp;REQUEST=DescribeFeatureType&amp;TYPENAME=*****
v&amp;SRSNAME=EPSG:3308&amp;username=****&amp;password=****"
xmlns:gml="http://www.opengis.net/gml"><Insert xmlns="http://www.opengis.net/wfs">
<nct_cov xmlns="****"><prop_name xmlns="****">test</*****><status
xmlns="****">lol</*****><geom xmlns="****"><gml:MultiPolygon
srsName="EPSG:3308"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs>
<gml:LinearRing><gml:coordinates cs="," ts=" ">COORDS HERE</gml:coordinates>
</gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember>
</gml:MultiPolygon></geom></*****></Insert></Transaction>


texts: Seeking documentation about non-iterative UTM to Lat-Long projection formulas

texts:
In

  Hager, John W., James F. Behensky, and Brad W. Drew. "The universal grids: Universal transverse mercator (utm) and universal polar stereographic (ups)." DMA technical manual 8358 (1989).

I found that the 2-6 formulas ("Conversion of grid coordinates to geographic coordinates"), in addition to not having any explicit dependency on the northing of the point to be projected, have the following note:

  Programmer's Note: The footpoint latitude $\phi'$ is normally derived by iteration based on term T1.

Now I found this implementation where no iteration is needed.

Where can I find some scientific reference for the no-iteration implementation?


texts: How to display ONLY local country names in google maps?

texts:
I am using google maps in an application where it is important that country labels display the country name in the local language.
If you pass in certain language parameters that google maps recognizes, but doesn't fully support (language=ka here), it displays localized place names with english labels underneath:

Which is great, except that I don't want the English labels at the bottom.
Is there any way to do this?


texts: Kernel Density with vector output and/or possible automation?

texts:
I'm currently working on a project that will have me calculate the kernel density of a number of shapefiles with varying parameters. My current workflow goes like this:

Import shapefiles in QGIS 2.10.1 from PostGIS 2.2.2
Run the Kernel density estimation command from the SAGA 2.1.2 toolbox
Vectorize the output with the r.to.vect command from the GRASS toolbox
Export the vector grid back to PostGIS for calculations

The problem is that the r.to.vect operation takes about 30-45 minutes per layer, and I will eventually be working with even larger files. Is there a way that I could obtain a vector grid directly instead of a raster? That would save me the useless raster-to-vector transformation. I could work with a solution in PostGIS, QGIS or even R.


texts: How to put label text in two rows using QGIS?

texts: I have a long text in one attribute field and want to place it as a label. Because the text is so long I want to divide it in two rows. How to do that using QGIS. I have a solution to divide text in two fields and make multi labeling, but I don't want that.


texts:
I own a project created with TileMill in SRS 900913. The data source is a Postgis database with the same SRS.

When I'm loading the file generated with TileMill (mbtiles) in QGIS and I'm putting under an orthophoto with the same SRS, or the original Postgis layer, the mapping displayed is moved approximately 2.5 meters moves in the Y axis relative to the PostGIS layer.

However if I'm loading QGIS directly from a PostGIS data source, works like a charm (think). Is like TileMill is displacing it because, as I said, loading the postgis layer directly from the database is working fine. In this image, you could see the Y offset between postgis (marked in blue), and mbtiles raster (marked in black)

Anyone know what could be happening?


texts:
I want to plot data from a project on a world map using choropleths. The data is at country level, and I am using ArcMap 10.2.

I had started by adding World countries (generalized) from ArcGIS Online, which essentially freezes the entire tool making unusable.

How can I import that dataset locally?

Can you point me to some tutorial/resource to draw a nice-looking world map with ArcMap?


texts: How to rotate a marker line in a polygon ring?

texts:
I've to rotate a marker line in a polygon ring by 180 degrees without rotating the same marker line at the polygon border (outline):

How can I do this in QGIS 2.18/2.99?


texts: How do I Use OpenLayers Multi-lingual capabilities?

texts:
So, OpenLayers has a i18n namespace.

But in my example at the bottom of the blog.

I am not getting the "Pan Up" which is what I would expect, but instead "PanUp", which is the string that should be translated in en.js.

I don't know enough about how this should work to debug it.  I couldn't find much in the documentation about this either.

How does OpenLayers know what language to use?  When I switch to French will it know to pull out the translations of "PanUp" from the fr.js file?


texts: Grids not originating from 2D parameter spaces

texts:
Is there any specific name for grids on sphere (or geoid, or any other manifold/surface) that do not originate from 2D parameter space?

For example generating a grid on a sphere using icosahedron tessellation do not involve any 2D parameters.

I'd like to read more about such grids. Any recommendations on literature?


texts: Using ArcMap to Split shapefile into multple shapefiles by attribute category?

texts:
I have a shapefile feature class that represents an area with signal strength per section of area.  It has signal strength from multiple radio sources but I would like to make separate file for each radio source.  For example if I start with the attribute table like this.

ALL SITES_SIGNAL LEVEL
| FID |  Shape  | Site_1 level | Site_2 level | Site_3 level |
|  0  | polygon |  -120        |   -140       |   -87        |
|  1  | polygon |  -142        |   -130       |   -77        |
|  3  | polygon |  -156        |   -111       |   -89        |
|  4  | polygon |  -60         |   -113       |   -112       |
|  5  | polygon |  -54         |   -87        |   -134       |
|  6  | polygon |  -90         |   -50        |   -117       |

I want three feature classes with attribute tables like this:

SITE_1 SIGNAL LEVEL
| FID |  Shape  | Site_1 level |
|  0  | polygon |  -120        |
|  1  | polygon |  -142        |
|  3  | polygon |  -156        |
|  4  | polygon |  -60         |
|  5  | polygon |  -54         |
|  6  | polygon |  -90         |

SITE_2 SIGNAL LEVEL
| FID |  Shape  | Site_1 level |
|  0  | polygon |  -140        |
|  1  | polygon |  -130        |
|  3  | polygon |  -111        |
|  4  | polygon |  -60         |
|  5  | polygon |  -54         |
|  6  | polygon |  -90         |

SITE_3 SIGNAL LEVEL
| FID |  Shape  | Site_3 level |
|  0  | polygon |  -87         |
|  1  | polygon |  -77         |
|  3  | polygon |  -89         |
|  4  | polygon |  -112        |
|  5  | polygon |  -134        |
|  6  | polygon |  -117        |

I've looked at split by attribute, but that doesn't do what I wasn't, in this example that gives me 7 tables.  I know I could make 3 copies of the table and delete two columns, but this doesn't seem efficient.


texts: How to add XY coordinates (latitude/longitude) to a DEM raster in ArcGIS 10.1

texts:
I would like to add xy coordinates of about 10 sites to a raster. I am using ArcGIS 10.1 and my xy coordinates are on Access 2013 database - saved as .mdb.

The steps I have taken so far:

1) In ArcCatalog - I tried making the OLE DB connection - Microsoft Jet 4.0. OlE Provider - when i clicked 'next' it says the connection is a 'success'

2) Then I clicked on the Database Connections on the arcCatalog tree, clicked on the 'OLE DB connection' and the table appears to be ok on the 'preview' screen. However when i click on the 'Description' section there is the error message:  'The item's XML contains errors' - I don't know what this means.

Anyways, I right clicked the table under OLE DB connection (on the arcCatalog tree) - clicked 'create feature class' - 'from XY'Table (i also put the right coordinate system - which is the British National Grid).

3) After that, I opened ArcMap. Added the DEM Raster (which also has the same coordinate system) and the XY feature class I just created.

But the problem is I cannot see the sites (the dots).

I have successfully done this in the past. But for some reason I just can't seem to get it to work this time. The following picture shows how I want it to appear:

I have also tried the other way: Arcmap - file - Add Data - Add XY data. And then right clicking the layer and clicking 'Data' - then 'Export data' to create a feature class.

A couple of weeks ago I did have a new version of Microsoft Access and excel installed ( the Microsoft office 2013 plus). I don't know whether this has something to do with it or whether I have simply missed a crucial step. If anyone could help me it would be massively appreciated!

Please help!!


texts: How to insert Google Maps API Lat/Long into PostgreSQL Postgis GEOMETRY(GEOMETRY, 4326)

texts:
I have very little knowledge and understanding when it comes to GIS (despite extensive experience with Google Maps API).

I've recently installed PostGIS (v2.1.5) with my PostgreSQL (v9.4) instance in order to begin storing and working with GIS data. I'm currently using Google Maps API v3's DrawingManager to draw polygons, circles, rectangles and markers which need to be stored in the database, and later retrieved to display the polygons again on Google Maps; These shapes will also be used for various aggregate queries.

I currently have two tables. The one is called "locations" which will store the polygons, circles and rectangles created with Google Maps' DrawingManager. The second, called "rawdata" will contain a set of records with a POINT (lat, long) inside of it.

First and foremost, I do not entirely understand the concept of projections and how they work, but, based on what I've read online is that calculations and accuracy may vary, significantly. I understand that Google Maps works with the Web Mercator projection which I believe is officially known as the EPSG:3857 projection (or CRS), while the preferred method for storing and working with geometry() shapes in PostGIS is in the EPSG:4326. As such and in reference to my create table statements below, I have chosen to store my geometry data in EPSG:4326 in the tables.

An easier method would to be store my geometry in the newer geographic format, but I do not wish to impose those performance and lack of functionality drawbacks as I will be working with a significant amount of data (hundreds of thousands, if not millions of records).

How should I store (and convert if applicable) the LAT/LONG coordinates which Google Maps supplies me with, assuming that my "geom" column is in the EPSG:4326 CRS? I found the simplest form was to use ST_GeomFromGeoJSON, but do I not know if this is the correct way or not.

INSERT INTO locations (label, geom) VALUES
(
    'Area Label',
    ST_GeomFromGeoJSON('{"type":"Polygon","crs":{"type":"name","properties":{"name":"EPSG:4326"}},"coordinates":[[[28.208620548248,-25.851267945548],[28.206732273102,-25.853353448277],[28.207848072052,-25.856597490521],[28.214156627655,-25.853469308462]]]}')
)

Here's where my confusion with the projections come in though. The lat/long coordinates I am passing through are retrieve from Google Maps API (which as I understand is in EPSG:3857 CRS), but I am supposedly storing this in the column as EPSG:4326... will this work?

Here's the table structures:

CREATE TABLE locations
(
    id BIGSERIAL,
    label VARCHAR(255) NOT NULL,
    geom geometry(geometry, 4326)
)  WITH (oids = false);

CREATE TABLE rawdata
(
    id BIGSERIAL,
    geom geometry(point, 4326)
);

ALTER TABLE locations_areas ALTER COLUMN geom TYPE geometry(geometry, 4326) USING ST_Transform(geom, 4326);


texts: Rule based labeling in rule based styles

texts:
I have got a layer with distance vectors between polygons. This layer has a rule based style where I can check each polygon to show only the vectors with contact to the polygon. Each vector has as an attribute the start and the endpoint (name of polygon)

I want to add a rule based labeling which should depend on the chosen stylename. Is there a function like "$get_active_stylename" which can be used as a variable?

The layerstyle Looks like this

There are some sub rules.

The polygons are buildings. So I want to chose building 86 and the label hase to be rulebased to building 86. The problem is that building 86 can be the start or the end of the vector.

QGIS 2.14


texts: How to write to geopandas dataframe directly to (Google) Cloud Storage?

texts:
I would like to write a geopandas.Dataframe directly to (Google) cloud storage, from the documentation fiona enables this, however I am not sure how to supply the correct **kwargs for the geopandas.Dataframe.to_file function. any examples?

import geopandas
gdf = geopandas.read_file("my-file.shp")
gdf.to_file("gs://my-bucket/my-file.shp", driver='ESRI Shapefile')


texts: Is it possible to build a route profile for a street with "real" heights considering tunnels and bridges?

texts: Is it possible to build a exact route profile for a street with "real" heights considering tunnels and bridges in an easy way with free data? Of course it is possible to remove thy waypoint located on ab bridge but for a longer route it's a lot of work. Is there any easy way to build a "street model"?


texts: OGR - Error trying to intersect two layers

texts:
I'm trying to intersect a buffer with a point shapefile using OGR, but when I run it, I get this error:

Traceback (most recent call last):

File "C:\Users\xxx\eclipse-workspace\xxx\Final.py", line 40, in <module>

buff_lyr.Intersection (shp2, result_lyr)

File "C:\Python27\ArcGIS10.5\lib\site-packages\osgeo\ogr.py", line 2263, in Intersection

return _ogr.Layer_Intersection(self, *args, **kwargs)

 TypeError: in method 'Layer_Intersection', argument 2 of type 'OGRLayerShadow *'

Here is my script:

import ogr

import os

import sys

path = 'C:/Users/xxx/Desktop/xxx/xxx/xxx/archivos'

os.chdir(path)

driverName = "ESRI Shapefile"

driver = ogr.GetDriverByName(driverName)

shp = 'avereproj.shp'

dataSource = driver.Open(shp, 0)

shp2 = 'Cober_point3.shp'

dataSource2 = driver.Open(shp2, 0)

numLayers = dataSource.GetLayerCount()

print 'Numero de capas: ', numLayers

layer = dataSource.GetLayer(0)

numFeatures = layer.GetFeatureCount()

print 'Numero de registros: ',numFeatures

spatialRef = layer.GetSpatialRef()

print spatialRef

memory_driver = ogr.GetDriverByName ('memory')

memory_ds = memory_driver.CreateDataSource ('temp')

buff_lyr = memory_ds.CreateLayer ('buffer')

buff_feat = ogr.Feature (buff_lyr.GetLayerDefn())

for ave_feat in layer:

    buff_geom = ave_feat.geometry().Buffer(100)

    tmp = buff_feat.SetGeometry(buff_geom)

    tmp = buff_lyr.CreateFeature (buff_feat)

result_lyr = memory_ds.CreateLayer('result')

buff_lyr.Intersection (shp2, result_lyr)

print ('Polygons: []'.format(result_lyr.GetFeatureCount()))

I have been trying to figure out how to fix this for a while but I can't find the problem.

Can you help me?


texts: Optimal projection to model solar shadows

texts:
I would like to map the shadows of mountains from a high resolution DEM. Therefore I am wondering if there is a coordinate system which would represent the projection of the rays of the sun as straight lines in the XY plane (image 1)

in order to easily compute the position of obstacles (image 2.

As mentioned in the comments, all light rays are parallel, so they should be represented as parallel lines on the map. I am not sure that it is not possible, so the less bad solution is welcome.


texts: Detect the coordinate reference system based on an image

texts:
Given a map as image, how can I find out which coordinate reference system it has?

The question is general, but here is a specific example.

Below is the map for which I want to find the coordinate reference system:

PDF source with much better quality. The map does not specify the coordinate reference system.

How do I find out the coordinate reference system used?


texts: GDAL 2.3.1 is installed but the llinux terminal is using GDAL 2.2.2

texts:
I asked this question before realising this site existed:

https://stackoverflow.com/questions/51476826/gdal-2-3-1-is-installed-but-the-llinux-terminal-is-using-gdal-2-2-2

I've installed GDAL 2.3.1 using pip on a ubuntu 16.04 in a python environment. The package is in the correct site-packages directory and yet when I run a python script, I receive this error:

Error 1: NUMPY driver was compiled against GDAL 2.3, but the current library
 version is 2.2

The error then ends in a segmentation fault (core dumped).

When using gdal-config --version, the output shows 2.2.2.

I want to know how to change the version of gdal that linux appears to be using from 2.2.2 to 2.3.1 but I have no idea how to do this.

UPDATE

As asked, here's the method I used to install gdal and the outputs. I followed the method in
Python GDAL package missing header file when installing via pip

I installed it in a python environment built to be able to run tensorflow 1.8.0 with Cuda 9.2 using a GPU - I don't know if that makes any difference

COMMAND

sudo apt-get install libgdal-dev

OUTPUT

Reading package lists...
Building dependency tree...
Reading state information...
Suggested packages:
  libgdal-doc
The following NEW packages will be installed:
  libgdal-dev
0 upgraded, 1 newly installed, 0 to remove and 3 not upgraded.
Need to get 0 B/6,496 kB of archives.
After this operation, 40.1 MB of additional disk space will be used.
Selecting previously unselected package libgdal-dev.
(Reading database ...
(Reading database ... 5%
(Reading database ... 10%
(Reading database ... 15%
(Reading database ... 20%
(Reading database ... 25%
(Reading database ... 30%
(Reading database ... 35%
(Reading database ... 40%
(Reading database ... 45%
(Reading database ... 50%
(Reading database ... 55%
(Reading database ... 60%
(Reading database ... 65%
(Reading database ... 70%
(Reading database ... 75%
(Reading database ... 80%
(Reading database ... 85%
(Reading database ... 90%
(Reading database ... 95%
(Reading database ... 100%
(Reading database ... 233309 files and directories currently installed.)
Preparing to unpack .../libgdal-dev_2.2.2+dfsg-1~xenial1_amd64.deb ...
Unpacking libgdal-dev (2.2.2+dfsg-1~xenial1) ...
Processing triggers for man-db (2.7.5-1) ...
Setting up libgdal-dev (2.2.2+dfsg-1~xenial1) ...

COMMAND

export CPLUS_INCLUDE_PATH=/usr/include/gdal
export C_INCLUDE_PATH=/usr/include/gdal

pip install gdal

OUTPUT

Collecting gdal
Installing collected packages: gdal
Successfully installed gdal-2.3.1

Could it be something to do with libgdal-dev installing 2.2.2?


texts: Geocoding Human Development Index (HDI) data

texts:
I am using QGIS for a class and I am trying to make a world map based on Human Development Index (HDI) data I got from the UN website. However, the data is in CSV and without latitude and longitude.

I already have a map of political boundaries of the world, but how would I orient the data so that it matches up with the political boundaries?

Also, the data I have is for multiple years and I would like to be able to show or not show a year's data by clicking a button if possible.

If there is HDI data available that includes latitude and longitude, I would be open to using that as long as it is recent.


texts: Clipping by mask results in blank output

texts:
After clipping by mask in QGIS 3.6 the output file has been blank (empty). Raster file and vector mask file have the same projected CRS:

USER:100029 - * Generated CRS (+proj=igh +a=6370997 +b=6370997 +units=m +no_defs)

How can I solve this problem?


texts: Possibility of converting from postgis to sqlite?

texts:
I have a large database with geodata for my country and I want to convert this into a usable format for Locus maps on Android so that it can be used for hiking.

Is it possible to convert this to sqlite? If so, how?

EDIT:

After a lengthy troubleshooting session in the comments below, ogr2gui works like a charm!


texts:
I'm trying to customize my layers from a python script. I would like to add several symbol layers to features.

# Creating symbol layers
symbolLayers = list()
lFill = QgsSimpleFillSymbolLayer.create({
    "color" : "{r},{g},{b},255".format(r = randomized_color[0],g = randomized_color[1],b = randomized_color[2])
})
symbolLayers.append(lFill)

# Create the fill symbol
symbol = QgsFillSymbol(symbolLayers)

# Adding to the current renderer
currentLayer.renderer().setSymbol(symbol)
currentLayer.triggerRepaint()

I have at this time only one symbol layer but when it will work, i will add more symbol layers.

It crash only when i add the symbol to the renderer. Before, it seems to work (setSymbol).


texts: ArcGIS layer to KML Python script gives TypeError: can only concatenate tuple (not "str") to tuple

texts:
I have some shapefiles to convert into KML file with a specific legend icon. So I tried to rearrange symbology and convert those shapefiles into .lyr files. After that, I put those files into a folder and used the following Python script.
import arcpy

arcpy.env.workspace = "E:/Layers"

# Set Local Variables
composite = 'NO_COMPOSITE'
pixels = 2048
dpi = 96
clamped = 'CLAMPED_TO_GROUND'

# Use the ListFiles method to identify all lyr and lyrx files in workspace
layers = arcpy.ListFiles("*.lyr*")

if len(layers) > 0:
    for layer in layers:
        # Strips the '.lyr(x)' part of the name and appends '.kmz'
        outKML = os.path.join(os.path.splitext(layer), ".kmz")
        for scale in range(10000):
            # Execute LayerToKML
            arcpy.LayerToKML_conversion(layer, outKML, scale, composite,
                                        '', pixels, dpi, clamped)
else:
    arcpy.AddMessage('There are no layer files in {}'.format(arcpy.env.workspace))

Which is a built-in code from the pro.arcgis.com resource site. But unfortunately, the program is showing errors and does not produce any result.
Error showed TypeError: can only concatenate tuple (not "str") to tuple
What am I doing wrong?


texts: CMV LayerControl Hide Dynamic sub-layers

texts:
I was trying to use the solution on this thread and my I am unable to get it to work. Specifically, layer._hideLayersInLegend does not exist anywhere in the scope of the LegendControl of my version of CMV, it is undefined. Do I need to add this somewhere?

CMV Configure Dynamic Layer to not display all sub layers


texts: I want migrate image tiles (osm2pgsql, mapnik, etc.) to vector tiles. I tried OpenMapTiles. But the documentation is bad. To start a local test system is no problem, but create a own server and own styling, is very laborious. What is a good alternative to OpenMapTiles?


texts:
Given files by client:

.dwg
.bak

QGIS doesn't like it; FREECAD doesn't like it; RhinoCAD doesn't like it  or Illustrator.

QGIS > Project > Import as .dwg/.dxf actually asks for .gdb, which doesn't help.

Found previous answer here: Converting .dwg to .shp using QGIS?

but it's not complete and from 2012.


texts: Writing expression syntax for Calculate Field in ArcPy?

texts:
I'm stuck on my CalculateField function with how to correctly type the syntax.

I'm not particularly tied to the variables I assigned up front.

The issue is syntax for the expression.

I am currently using ArcGIS 10.3 and Python 2.7.8.

This is the line that is getting me:

arcpy.CalculateField_management(fc, "Cohort", expression="""openYear + " to " + closeYear""")

How should I be writing the expression?

expression="'str(openYear)' + " to " + 'str(closeYear)'"

The rest is a sample of my loop statement.

#################################
for fc in fcList:
    openYear=fc[12:16]
    closeYear=fc[20:24]
    # Add "Cohort" field to *USA feature classes in scratch.mdb
    print fc
    arcpy.AddField_management(fc, "Cohort", "TEXT")
    print "Cohort added."
    # Calculate "Cohort" field to match feature class name in scratch.mdb
    arcpy.CalculateField_management(fc, "Cohort", expression="""openYear + " to " + closeYear""")


texts: Setting print scaled factor of map in QGIS?

texts:
I want to print a large area map in QGIS with 1:50000 scale but I need to fit it in 30 inches paper. How can i zoom out the map and retain its 1:50000 scale so that it will fit in my paper?

1) This is scaled in 1:50000

2) and this is zoomed out in 75% but the scale is still 1:50000


texts: How to share custom coordinate in an atlas between pages?

texts:
I know we can build custom coordinate grid (cf. url:How to create a custom coordinate grid in the QGIS 2.14.2?).

I would like to use this in an atlas but I have not succeded yet.

I get

Page 1 :
ABCDEF
1
2
3
4

Page 2 :
ABCDEF
1
2
3
4

While I would like to have a continuity in letter from page 1 to 2 :

Expected Page 2 :
GHIJKL
1
2
3
4

Any idea how to proceed?


texts: What's the strongest refresh method available to ArcEngine?

texts:
I have five layers: polygon, line, and point with a single IRasterLayer file.  I am programmatically setting the ILayer.Visible attribute for the Raster but the layer is not refreshing.  If I zoom in, or zoom out, then it refreshes properly.

I tried these two methods:

mapControl.Refresh();
mapControl.Refresh(esriViewDrawPhase.esriViewAll, null, null);

But they didn't refresh my map.  What method is stronger that this for refreshing the entire map?


texts: Is there a way to open attribute tables undocked by default in ArcGIS Pro 2.2?

texts:
I can undock attribute tables by dragging them out of the application window to another monitor. Field view, subtypes view, etc. then opens in that undocked window. However, when i open another attribute table it annoyingly opens inside the main application window. Again, I can drag that attribute table out of the application and dock it by itself or into my already open undocked attribute window.

Is there a way to change the default behaviour so that all attribute tables open undocked? Ideally, is there a way to open all attribute tables into a single undocked window as tabs - same as ArcMap.  Since I can do this manually by dragging things around I figure there should be a way.


texts: How to remove negative values from a created HAND model

texts:
I have used ArcGIS to create a HAND (height above nearest drainage) using an SRTM DEM.

However, some small portions have negative values in the model which is supposed to have zeros as the minimum value.

Can anyone help me sort this problem out?


texts:
I'm fairly new to GIS, but I've got the basics down. I'm utilizing QGIS and Python to calculate election statistics by legislative and congressional districts. I'm trying to aggregate precinct level results up to the district level but have a small problem. In some districts (about 10%) precincts are split unevenly into 2 or more districts.

The good news is that I have census blocks that are intact across the board even among precincts. Is there anyway I can take the population statistics from the block level to disaggregate the precinct election returns down to the census block level in QGIS? If not, are there tools in Python or PostGIS that will help?


texts: Deleting attribute records in one project also deletes it in another project?

texts:
Because I am learning QGIS and I had to delete records in a attribute table, I made a copy (project > save as > saving as 'project_copy' before deleting attributes. To be safe...

So I deleted the attributes in the project_copy.
But to my surprise, I see that the deleted attributes are also deleted in the original project!? There seems to be a link between both projects although I didn't ask for...

Did I really lose the deleted attributes?

Is there a way to restore the deleted attributes?


texts:
I have three points (boreholes) on the map. I created them using shape file layer. Each borehole shape file layer shows the strati-graphic section of the ground.
(i.e. between 0-50 m = Gravel , 50-60 m = Sand , 60-80 m = Clay.)
Rocks are located in different depths for different boreholes. So I need to interpolate these three points to create continuous surface with 3 boreholes.

I'm pretty new to QGIS, so it's better if someone could give me the instructions in detail since I tried to find what I look for from the questions that were asked before, but couldn't.

I want to create a continuous surface from 3 points of different borehole data to find the intermediate. I want to assume that the geological structure under the ground is same everywhere. I need to interpolate these three data sets to find an average subsurface model.

These are the boreholes:

The thing that I do was to click interpolation under "raster", and then I chose one of the boreholes & added it. For cellsize X and Y I put 300, after I entered an output file name, I clicked OK but nothing comes to the layers panel. I tried the same thing for both TIN and IDW, but it didn't work.

If I add another thing (For example Heidelberg contour) instead of adding borehole that I created, then it works, at least something comes to the layers panel after I click OK. From the image you see below, nothing comes to the layers panel when I click OK.

And this is the attribute table for Borehole_A :


texts: QGIS: Exported layer values differ from the underlying raster layer they were exported from

texts:
I´m trying to merge multiple raster layers (partial overlap) and then export them into a single layer for my required area of interest (AOI). When i use the raster merge (gdal_merge) option it merges all the rasters into one raster, in areas where there is overlap values from the last image is used(https://gdal.org/programs/gdal_merge.html), which i understand. The merged raster has only 4 bands (R, G, B and NIR) the same as the original rasters.
But when i try to export this image for my AOI using extents i get slightly different raster values in some areas. Why is this happening? Shouldn´t the raster values be exactly the same?
I export it in Raw data output mode, GeoTIFF format, with a custom extent and no compression. Let me know if anything is unclear.


texts: Do any maps exist that allow me to get lat/lng from a point with no derivative-licensing issues?

texts:
I need a map I can display that lets the user click on a point then returns the Lat and Lng.

EDIT TO CLARIFY @Igor: Yes, the user is clicking a specific feature on the map so I guess it would constitute "tracing over imagery".

EDIT: The returned lat/lng can be fairly inaccurate, say even 10m or 20m away.

The issue is that I need the Lat/Lng data to be freely usable with no legal restriction.

For instance, from http://wiki.openstreetmap.org/wiki/Legal_FAQ

  If I have data derived from OSM data, do I have to distribute it?

  The licence does not force you to distribute or make any data
  available. But if you do choose to distribute it, or anything derived
  from it, it must be under the same licence terms as the OSM data.

The area to be mapped is Edinburgh, UK. Ordnance Survey is definitely not an option.

Is this going to be possible?

(We are actually going to be doing an open data service, so we may be able to use Open Street Map data, but I would like to see if there is a freely open option first.)

Thanks,
James


texts: Make background image visible at all scales in webmap QGIS2web

texts:
I'm trying to produce a web map using the QGIS2web plug in and leaflet output, but the background image is not visible when I zoom in beyond some scale.  Is there a way to control the scales a background image is visible?

My map is set up with 5 vector layers displayed over a Google Satellite image.  It exports quickly to Leaflet, but the background layer disappears when zoomed in too much.

The background shows up fine when zoomed out

But disappears when you zoom in.

I've tried exporting when zoomed in to different scales, but that does not appear to have any effect on the output. I haven't been able to find any settings in the Export to Web window that would allow me to change the visibilty settings, as there are for the vector layers.


texts:
Is there any way to export a region on Google maps to a shapefile or any other format?

For example, I have the following region on Google Maps, and I wish to export it to shapefile.


texts:
I installed OpenLayers plugin successfully in QGIS 2.2 Valmiera, but when I add OpenStreetMap layer, the following error message appeared as a Python error popup window:

An error has occurred while executing Python code:

> Traceback (most recent call last):   File "C:/Users/I Wayan
> Mudita/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py",
> line 50, in addLayer
>     self.__plugin.addLayer(self)   File "C:/Users/I Wayan Mudita/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py",
> line 202, in addLayer
>     self.__setMapSrsGoogle()   File "C:/Users/I Wayan Mudita/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py",
> line 261, in __setMapSrsGoogle
>     extMap = coodTrans.transform(extMap, QgsCoordinateTransform.ForwardTransform) QgsCsException: forward
> transform of (-985846.534752, -198190.645686) PROJ.4: +proj=longlat
> +datum=WGS84 +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs Error: latitude or longitude exceeded limits

Python version:

  2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)]

QGIS version:

  2.2.0-Valmiera Valmiera, c3a2817

Python path:

  ['C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins\processing',
  'C:\Users\I Wayan Mudita\.qgis2\python\plugins\GeoCoding',
  'C:/PROGRA~1/QGISVA~1/apps/qgis/./python', u'C:/Users/I Wayan
  Mudita/.qgis2/python', u'C:/Users/I Wayan
  Mudita/.qgis2/python/plugins',
  'C:/PROGRA~1/QGISVA~1/apps/qgis/./python/plugins',
  'C:\PROGRA~1\QGISVA~1\bin\python27.zip',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\DLLs',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\plat-win',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\lib-tk',
  'C:\PROGRA~1\QGISVA~1\bin',
  'C:\PROGRA~1\QGISVA~1\apps\Python27',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\PIL',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\win32',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\win32\lib',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\Pythonwin',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win32.egg',
  'C:\PROGRA~1\QGISVA~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode',
  'C:\Users\I Wayan
  Mudita\.qgis2\python\plugins\DigitizingTools\tools',
  'C:\PROGRA~1\QGISVA~1\apps\qgis\python\plugins\fTools\tools']

What is wrong here?


texts:
Is it possible to use MapObjects Java with ArcIMS 10? We are starting to look at migrating our MOJava solution to run against ArcIMS 10, has anyone been successful in migrating to ArcIMS 10?

Our initial testing is indicating that there are problems with some of our feature layers. The layers don't return any geometries to the MOJava client. But the layers workes in ArcMap.

We are currrently runnig on ArcIMS 9.2 so any migrations issues for ArcIMS 9.3 might be intressting as well. According to the what's new in ArcIMS 10 there are no changes to ArcXML.

Update: Using MessageConsole in JoViewPlus gives this error message against the ArcIMS 10 layers:

> 2010-okt-31 11:12:08:    com.esri.mo2.client.aims.d
com.esri.mo2.client.aims.InvalidStreamDataException: A general exception has occurred: (Invalid data has been detected within the FeatureStream: <' "<?xml version="1.0" encoding="UTF-8"?>
<ARCXML version="1.1">
<RESPONSE>
<ERROR machine="hsovimsu5" processid="2848" threadid="2940">AXLParser: POLYGON tag incorrect.</ERROR>
</RESPONSE>
</ARCXML")
> 2010-okt-31 11:12:08:    com.esri.mo2.client.aims.d
com.esri.mo2.client.aims.InvalidStreamDataException: A general exception has occurred: (Invalid data has been detected within the FeatureStream: <' "<?xml version="1.0" encoding="UTF-8"?>
<ARCXML version="1.1">
<RESPONSE>
<ERROR machine="hsovimsu5" processid="2848" threadid="2968">AXLParser: POLYGON tag incorrect.</ERROR>
</RESPONSE>
</ARCXML")
> 2010-okt-31 11:12:08:    com.esri.mo2.client.aims.d
com.esri.mo2.client.aims.InvalidStreamDataException: A general exception has occurred: (Invalid data has been detected within the FeatureStream: <' "<?xml version="1.0" encoding="UTF-8"?>
<ARCXML version="1.1">
<RESPONSE>
<ERROR machine="hsovimsu5" processid="2848" threadid="2940">AXLParser: POLYGON tag incorrect.</ERROR>
</RESPONSE>
</ARCXML")
> 2010-okt-31 11:12:08:    com.esri.mo2.client.aims.d
com.esri.mo2.client.aims.InvalidStreamDataException: A general exception has occurred: (Invalid data has been detected within the FeatureStream: <' "<?xml version="1.0" encoding="UTF-8"?>
<ARCXML version="1.1">
<RESPONSE>
<ERROR machine="hsovimsu5" processid="2848" threadid="2968">AXLParser: POLYGON tag incorrect.</ERROR>
</RESPONSE>
</ARCXML")

More info: Talked to a colleague about this and we only have problems on some of our feature layers, all our images layers are ok. It seems to be connected to the scale of the map because sometimes it shows up when zooming in to different scale intervals. We have removed all scale dependent rendering, so it's not an renderer issue causing this. More likely that it might be some kind of generalization involved. Also the layer we are having problems with is a layer with multi part polygons. Has anyone been able to successfully migrate to ArcIMS 10?


texts: Schema error while uploading raster in PostgreSQL

texts:
Tried many ways
A.
import os
import subprocess
from subprocess import PIPE, run

db_name = 'agro_dss'
db_host = 'localhost'
db_port = '5432'
db_user = 'postgres'
db_password = 'imd123#'

os.environ['PATH'] = r';C:\Program Files\PostgreSQL\13\bin'
os.environ['PGHOST'] = db_host
os.environ['PGPORT'] = db_port
os.environ['PGUSER'] = db_user
os.environ['PGPASSWORD'] = db_password
os.environ['PGDATABASE'] = db_name

cmd = 'raster2pgsql -s 4326 -C -t auto -M D:\Advisory\Raster\*.tif -F public.polys | psql -U {} -d {} -h {} -p {}'.format(db_user,db_name,db_host,db_port)

subprocess.call(cmd, shell=True)
result = run(cmd, stdout=PIPE, stderr=PIPE, universal_newlines=True)
print(result.returncode, result.stdout, result.stderr)

The error is
1  ERROR: Unable to read raster file: public.polys

B.
from osgeo import gdal, osr
import psycopg2
import subprocess
import sys, os

fileName = r"D:\Advisory\Raster\MODIS_All_India_NDVI_2001_01_01.tif"
check_proj = gdal.Open(fileName)
proj = osr.SpatialReference(wkt=check_proj.GetProjection())
projection=str(proj.GetAttrValue('AUTHORITY',1))

input_path = r"D:\Advisory\Raster\\"

os.environ['PATH'] = r';C:\Program Files\PostgreSQL\13\bin'
os.environ['PGHOST'] = 'localhost'
os.environ['PGPORT'] = '5432'
os.environ['PGUSER'] = 'postgres'
os.environ['PGPASSWORD'] = 'imd123#'
os.environ['PGDATABASE'] = 'agro_dss'

for raster in os.listdir(input_path):
    if raster.endswith(".tif"):
       name = raster.split(".tif")[0]
       raster = os.path.join(input_path, raster)
       rastername = str(name)
       rasterlayer = rastername.lower()
       conn = psycopg2.connect(database="agro_dss", user="postgres", host="localhost", password="imd123#", port=5432)
       cursor = conn.cursor()
       print(cursor," connection suscessfull")
       cmds = 'raster2pgsql -s '+projection+' -I -C -M "'+raster+'" -F -t auto'+' public.'+name+' | psql -d agro_dss -U postgres -p 5432 -h localhost'
       subprocess.call(cmds, shell=True)
       result = run(cmds, stdout=PIPE, stderr=PIPE, universal_newlines=True)
       print(result.returncode, result.stdout, result.stderr)

The error is
1  ERROR: Unable to read raster file: public.polys

Sought many answers:

Using raster2pgsql in Python console of QGIS?
Import raster into PostGIS using Python
Upload a raster in PostGIS
Unable to read TIFF file with raster2pgsql

but did not come up with the solution
How to rectify this?


texts: Copy/Paste Layer with Attribute table to new Layer using ArcGIS Pro?

texts:
When I copy features from one layer in ArcGIS Pro to a blank layer, the polygons are transferred, but the attribute table behind those features is not transferred.

To give a specific example, I downloaded a counties file from a public data site (Ohio Department of Transportation). Small background: I work with a specific subset of counties within the state, so I want to create a basemap that allows me to filter just on those subset of counties. I also want an exploratory map that is the entire state. So, with the counties file, I want to be able to create a layer that is just my subset of counties, and another layer that merges all our counties into just the state boundary.

This first image shows all the counties along with the attribute table headers.

When I select all the features from the layer above and paste that into a blank .shp layer, I get this:

In that example, the attribute table did not carry over, instead it makes me a virtually blank table. I know that FID is associated with each feature, but then it just has Polygon for the entire shape column (which makes sense) and "0" for each row in the ID column.

I get no options on either table to add or edit fields. I feel as though I am clearly doing this operation incorrectly.

What are the correct steps to create a new layer from a set of selected features so that attributes are retained? Is this not possible in ArcGIS Pro the way it was in ArcMap?