texts: Query containing special characters passed by pyodbc works differently in standalone ArcGIS script and toolbox

texts:
I am passing few queries through pyodbc in my script to get data from Access database and write them to layers table. Everything was fine when I have been running script from PyCharm, but when I put it in toolbox and tried to run it in ArcGIS, queries containing national diactric characters stopped working. Here is an example of such query:

expression_2 = """SELECT F_ARODES.ARODES_INT_NUM, F_STOREY_SPECIES.SPECIES_CD, F_STOREY_SPECIES.SPECIES_AGE
            INTO fill_wydzielenia_2_temp
            FROM F_GROUP_CATEGORY INNER JOIN ((F_ARODES INNER JOIN F_SUBAREA ON F_ARODES.ARODES_INT_NUM = F_SUBAREA.ARODES_INT_NUM) INNER JOIN (F_AROD_STOREY INNER JOIN F_STOREY_SPECIES ON (F_AROD_STOREY.ARODES_INT_NUM = F_STOREY_SPECIES.ARODES_INT_NUM) AND (F_AROD_STOREY.STOREY_CD = F_STOREY_SPECIES.STOREY_CD)) ON F_SUBAREA.ARODES_INT_NUM = F_AROD_STOREY.ARODES_INT_NUM) ON F_GROUP_CATEGORY.AREA_TYPE_CD = F_SUBAREA.AREA_TYPE_CD
            WHERE (((F_ARODES.TEMP_ACT_ADRESS)=?) AND ((F_AROD_STOREY.STOREY_RANK_ORDER)=?) AND ((F_STOREY_SPECIES.SPECIES_RANK_ORDER)=?) AND ((F_GROUP_CATEGORY.SUPERGR_CAT_NAME)=?))"""

db_cursor.execute(expression_2, (True, 1, 1, u"grunty leśne zalesione"))
db_cursor.commit()

Created table "fill_wydzielenia_2_temp" is empty when script is run from ArcGIS, but is correctly populated when I run the script from an IDE.
I assume the argument with special characters must be treated differently, but I don't know the cause. How can I fix this?


texts: Using aggregate function in r with function from package landscape metrics

texts:
I am trying to combine the information from two raster data sets: The first data set contains information on the numer of fotos uploaded on flickr in a given raster cell (ot_fl1k). The second dataset contains the information about landscape characteristics for the same area from the Corine Landcover Inventory (ot_cor12). The difference between the two is that they exhibit a different coordinate reference system and the Corine data display a much finer resolution. My ultimate goal is to aggregate the Corine raster data to the level of the Flickr data based on the number of patches from the landscape metrics package as aggregation method. I proceed as follows:
ot_fl1k <- raster("Osttirol_flickr_1k.tif")
crs(ot_fl1k)

CRS arguments:
+proj=aea +lat_1=32.5 +lat_2=54.5 +lat_0=51.4
+lon_0=22.65 +x_0=0 +y_0=0 +ellps=WGS72
+towgs84=0,0,4.5,0,0,0.554,0.2263 +units=m
+no_defs
res(ot_fl1k)

[1] 1000 1000
Hence, my Flickr data exhibit a 1000x1000 meters resolution. As my next step I import the Corine dataset:
ot_cor12 <- raster("Osttirol_Cor12.tif")
crs(ot_cor12)

CRS arguments:
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000
+y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs
res(ot_cor12)

[1] 100 100
So the CSR and the resolution of the two raster data sets are different. In order to align them I reproject the Flickr dataset:
ot_fl1k_rp <- projectRaster(ot_fl1k,crs = crs(ot_cor12), res=1000)
crs(ot_fl1k_rp)

CRS arguments:
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000
+y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs
res(ot_fl1k_rp)

[1] 1000 1000
So now the two exhibit the same CRS but I have two problems:
First, I seem to have almost doubled the number of cells in the reprojection process and thereby also distorted the overall number of fotos uploaded in the area:
ncell(ot_fl1k)

[1] 2156
ncell(ot_fl1k_rp)

[1] 4026
In a second step I would like to aggregate the resolution of the Corinne dataset to the one in the reprojected Flick dataset by applying the lsm_l_np command for computing the number of patches in an area. If I apply this to the area as a whole I get the following:
lsm_l_np(ot_cor12)

A tibble: 1 x 6
layer level     class    id metric value

1     1 landscape    NA    NA np       570
So basically I obtain 570 patches in the area overall. However, I do not want to compute this information for the whole area but rather for sub-areas of 1 km^2 corresponding to the cells of this size in the Flickr dataset. So basically I would like to use the aggregate command but using lsm_l_np instead of mean or mode. Something like this:
ot_fl1krp_agg <- aggregate(ot_cor12rp, fact=10, fun=lsm_l_np)

I realize that this does not work due to the function not being specified for this command, any suggestion on doing the reprojection accurately and aggregating the raster based on the indicated procedure would be greatly appreciated!


texts:
I'm very new to the world of GIS, but I'm learning a great deal from this site, as well as others.  I fell like this should be 'easy' but I need to  be pointed in the right direction.  I'm trying to map a railcar's route based on a shapefile of the railroad lines in America.

What I have:

1) I have a file containing specific railcar's routing.  Railcars have Radio Frequency ID tags (RFID) that are scanned when entering a rail yard.  This file contains all these scanned entry points as well as Lat/Log data for the yard.

2) I have a shapefile of all railroads in America.

My Problem:

I want the path of the railcar to follow the railroad shapefile exactly.  If I just tried to map the Railcar file, it would show direct lines from one yard to another, and would not follow the actual path of the railroad.  How do I do this?  Again, I'm just looking to be pointed in the right direction.


texts: Troubleshooting slow data access to shapefiles by ArcGIS Desktop?

texts:
We have just upgraded from ArcGIS Desktop 9.2 to 9.3.1

For a ~250MB .shp in 9.2 it took 1 minute to render, at 9.3.1 it now takes 50 minutes.

The format has to stay as shapefile.
The data is stored on both NetApp and Dell Servers - problem consistent on both.
Only an issue with GIS data, in particular shapefiles and only in 9.3.1.
Have already  have tested during the day and night and the times returned are consistent showing that network traffic is having no or very little effect on the poor performance.
Also tested on servers that do not have virus protection to eliminate that from the equation.
Already tried testing on several servers with differing switches (100Mbit and 1Gbit) but the results have been the same.
Have tested on a newly built pc with a vanilla install of ArcGIS so it can’t be incorrectly uninstalled components or third party products causing the issue.

I am putting together some further test-cases to carry out, to try and pinpoint the problem further before taking it to Esri:

Run the tests at 9.2 with a large ESRI sample dataset.  Run tests with it as a shp and as a file GDB.
Re-test with 9.3.1 (+SPs) and again with v10
Run tests in both ArcMap, ArcGIS Explorer and ArcReader (to see if it is specific to an application)
Compare these results with how long it takes to open when stored locally?

From what I know, there have not been any changes made to the underlying code and the way that ArcGIS connects to a network between 9.3.1 and 9.2. I also don't think that Esri  changed anything in the way that ArcGIS reads shapefiles between these versions?

Other than the above, the only other suggestion I can think of is to run WireShark to try and see how long it takes to send/receive packets from the client machine & the server.

This will tell me two things:

1) How many packets are going both directions

2) how long each packet takes to be transmitted.

My thinking is, if it is related to 9.3.1, then you would most likely see an increase in the number of packets. It could be a security configuration that at 9.3.1 we have encountered, but at 9.2 we did not.
It would be best to capture a wireshark trace for both 9.2 and 9.3.1 accessing the same shape file on the same server.

Has anyone experienced anything similar to this?

Does anyone have any other ideas on how to troubleshoot this?


texts: How to print table information using mapbasic

texts:
How do I print a tables information in Mapbasic? I have a snippet of code - I want to select all features within a buffer and then print the results as a message in mapinfo that can be copied out into word.

Select * from Land_charges_test
Create Object As Buffer From Selection Width 500 Units "m" Type Spherical
   Resolution 100 Into Table Land_charges_test Data ID=ID,Tree=Tree

Select * from tree, Land_charges_test
   where tree.Obj Within Land_charges_test.Obj into Selection
Browse * From Selection

**Print Selection**

It's the print selection function I need, but cannot find


texts:
I am trying to use the Resample tool in Spatial Analyst to  go from a 1kmx1km grid to a 500mx500m grid for increased resolution of my land use data. My cell size is 0.008, and I want a 0.004 cell size. I input this into the Resample tool, and it runs smoothly, but when I check the cell size after the resample, it is still 0.008.

Does anyone know why this is? I am using ArcMap 10.3.


texts: Distribute gridded data to corresponding subbasin boundaries?

texts: I have one ArcGIS shapefile of my river basin consisting of many subbasin. Another shapefile showing the location of gridded points over the whole basin (Its latitude and longitude are used as file name of corresponding data file). I'd like to distribute that gridded data corresponding to each subwatershed. Basically, the idea is to know which gridded points belong to which subbasin either inside the boundary or within a certain limit from the boundary. Anybody has an idea how to do this?


texts: Two different cell sizes represented in a raster

texts:
I am currently working with a DEM with a cell size of 10 meters. I only needed a certain area of this DEM and wanted it at a 1 m cell size so I extracted by mask. I set the appropriate environmental settings to obtain a cell size of 1 m and my mask layer was also set to 1 m cell size (same coordinate system as well).

My output does indeed have a cell size of 1 m, but you can see from the following image two different sized cells. Along the edges measure 1 m, but the cells in the middle measure about 8 meters (this is the cell size of the original DEM before I clipped it).

I'm wondering if my clipped DEM does in fact have a 1 m cell size or is something wrong here?

The point of setting my cell size at 1m was to obtain a smooth look. Possibly those larger cells really are 8 1m cells, all of the same values? I have tried resampling with bilinear interpolation and I get the exact same result.

I am using ArcGIS 10.5.1 Version 10.5.1.7333.

The most similar question I have been able to find, but which did not seem to have a final solution/one that solves my particular problem is Cell size (10m) in Raster properties differs from measured cell size (20m) in multiband (Sentinel) raster?


texts:
I have a problem with importing a list of coordinates to ArcMap. I have a shape file of regions of Chile and I would like to add a list of locations with coordinates I got from Google Maps.

The shape file with the Chilean regions is coming from http://sit.conaf.cl/. I open it with ArcMap and from the property tab I get the following information:
Projected Coordinate System: WGS 84 / UTM zone 19S.
Projection: Transverse_Mercator.
For instance, one point from the respective region is 223.205,971  5.545.419,924 Meters.

My locations from Google Maps look like -39,812159 / -73,178511.
I guess I just need to pick the right coordinate system when importing but I cannot not figure out which one that would be. Whatever system I choose, the points will not show up.

If I import this list using the Add XY data feature to a new project and select a random coordinate system, the points will show up. So I think that there is nothing wrong with the spreadsheet.


texts: Vector layer feature destroys after geometry transform

texts:
I want to get coordinates of a drawn polygon.

var myLatLonSquare = vectorLayer.features[0].geometry.transform(map.getProjectionObject(),
                     new OpenLayers.Projection("EPSG:4326"));

After zooming, the drawn polygon disappears because of geometry transform.

How can I solve this problem without disappearing drawn polygons?


texts: How to avoid duplicate GIDs when copying features from one PostGIS layer to another?

texts:
When I enter new elements (features) in another Postgres layer, I
can do it in two ways:

Drawing new element(s) (with 'Add feature') which  I do  rarely or
Copying (or cutting) some elements from another Postgre layer
(Source layer), and paste it into target layer which I do frequently

In first example, saving of edits works normal because this layer
gets gid from postgre database sequence
*nextval('layer_name_gid_seq'::regclass)*

In second example, I got an error during saving edits, because
while copying element from source layer to target layer qgis copied
gid of element from source layer. When attempting to save edits
this error is returned:

  Could not commit changes to layer „Cjevovodi“
  Errors: ERROR: 1 feature(s) not added.
  Provider errors:
  PostGIS error while adding features: ERROR: duplicate key value violates unique constraint "cjevovodi_okill_pkey"
  DETAIL: Key (gid)=(5) already exists.

I tried to copy *nextval('layer_name_gid_seq'::regclass)* in
field gid, but this sequence can`t be pasted in field gid as field is
defined as numeric.

Does anyone know simple way to copy elements from source layer
(with existing gid) assign new gid?

Thanks!


texts: QGIS georeference a Mesh layer with wind vectors?

texts: I have an .nc file with u and v components to create a map of wind speed and direction. When I add the layer as a mesh in QGIS, I am able to get the arrow vectors correctly (image below), however, the mesh does not map to the underlying basemap. Likely because the lat and lon are incorrect. Is there a way to georeference this layer manually?


texts: Vector Union function performance in QGIS and ArcGIS DEsktop

texts:
I am trying to find the union of two vector layers (both spatial indexed) - one land use polygons and another grid. In QGIS it is running since 4 days ago, but in ArcGIS Desktop it took less than 10 min.

Does anyone know if ArcGIS Desktop pre-processes vector(s) before actual union function?

The time difference cannot be this drastic!


texts: QGIS - Export PDF map with georeferenced data

texts:
I am used to working with ArcMap but have been playing around with QGIS.
When I export a PDF map with ArcMap, the PDF is spatially referenced. This is a handy tool, I can put the PDF map on my iPhone and it tracks my location on the map (even when I'm out in the woods away from cell coverage). It's very handy for forestry work.

Are the PDF maps exported with QGIS also spatially referenced? Is there a plugin that allows for this feature?


texts: Can't find feature classes when creating SDE connection file with credentials

texts:
EDIT: When I try to access the workspace properties of the connection file, it says that the method workspaceType doesn't exist for describing data. When I use the regular 'Database Connections\myDB.sde' connection I can access the workspace properties just fine. When I run desc.dataType on it it comes back as "File"...how do I make it a workspace? The connection file is obviously not being treated as a workspace, and that's the first issue..

Using ArcGIS 10.5. SQL Server 2014, version 12.0.5207.0 (64 bit)

I am using the CreateDatabaseConnection geoprocessing tool to create an SDE connection file and list feature classes within it. I am using this tool instead of the more typical Database Connections\myDB.sde because I need to supply a username and password for connection by external users when this project is published up on a server and people run this script via a web interface. They won't have my windows authentication, thus the need for database authentication.

When using this script, I receive an exit code of 0 but it will not list the feature classes that I know exist in the database. I can see that it creates the SDE connection file. I have given the database user read and write membership using SQL Server Management Studio, so I'm not sure if it's a permissions issue.

import arcpy

arcpy.CreateDatabaseConnection_management(out_folder_path='C:\Users\xxx\Desktop\GIS_Testing', out_name='HBM.sde', database_platform='SQL_SERVER',
                                          instance='sde:sqlserver:myServerName', account_authentication='DATABASE_AUTH', username='HBMWriter', password='xxx',
                                          save_user_pass='SAVE_USERNAME', database='HBMTest')
arcpy.env.workspace = r'C:\Users\xxx\Desktop\GIS_Testing\HBMTest.sde'
featClasses = arcpy.ListFeatureClasses()
for fc in featClasses:
    print fc

It's worth noting that when I try to get some workspace properties on this SDE connection file, I get back: AttributeError: DescribeData: Method connectionProperties does not exist.

If I try to access connectionProperties of the regular Database Connections\myDB.sde it'll give me the properties so it's as if it's not treating the connection file as a workspace.


texts: proxy.cgi for GeoServer is not accessible in browser

texts:
I have a problem with WFS request:

  var wfs = new OpenLayers.Layer.Vector("WFS",{
  strategies: [new OpenLayers.Strategy.Fixed()],
  projection: new OpenLayers.Projection("EPSG:4326"),
  protocol: new OpenLayers.Protocol.WFS({
  version: "1.1.0",
  url : "http://localhost:8080/geoserver/wfs",//
  featureNS : "http://www.opengeospatial.net/cite", // workspace namespace URI
  geometryName: "geom",  // field in Feature Type details with type "Geometry"
  featurePrefix: 'cite',    //geoserver workspace name
  featureType: "InfoGeometric"        //geoserver Layer Name
}),

In particular in chrome console i have the error POST http:// localhost:8081/cgi-bin/proxy.cgi?url=http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs 404 (Not Found). I have create in C:\Program Files\apache-tomcat-6.0.33\lib the folder cgi-bin and in it the file proxy.cgi (i added in allowedHosts variable values'localhost','localhost:8080','localhost:8081',               'localhost:8080/geoserver/wfs'). In OpenLayers.js i set OpenLayers.ProxyHost="/cgi-bin/proxy.cgi?url="; Now, if i go on http:// localhost:8080/geoserver/wfs/DescribeFeatureType?version=1.1.0&typename=cite:InfoGeometric i have a description in xml of my InfoGeometric but if i go on  http:// localhost:8081/cgi-bin/proxy.cgi i don't connect to it. I don't know i must do. Thanks very much.


texts: QGIS: Style Point displacement brings too much together

texts:
I want points which are close together, about 10 mm from each other on the printed map or 500m in real world to be arranged in a circle around their center point. In other words I use Point Displacement in Layer Properties | Style. I have played around with the point distance tolerance number (mm, pixel and map units), but if I at all allow points to be displace (>0) points from a much too large area are displaced. If I for instance use 0,5mm as the Point distance tolerance, points from at least 10km away are put together. Is there something else wrong / som other settings that I need to adjust, or what is it that I do not understand.

How it should be, with these points centered around their real place (the red star):

How the points are shown, if I enter any number bigger than 0 in the Point displacement tolerance box (points from a point a few km away are included, and the red star is shown somewhere between these two points) :


texts: QGIS -- Extract Values from KDE raster to orginal point data layer

texts:
After performing a heatmap on QGIS of my cheetah location dataset, I would like to extract values from my kernel density estimate (KDE) raster data layer to my orignal point data layer used to create it in order to extract 50% of my values from my KDE.

I went on the POINT SAMPLING TOOL in the puglin toolbars. I only have one choice for the Layer containing sampling points --> Cheetah_Coalition (discrete points layer), and if I click on the "Clipped (mask): Band 1 (raster)" within the Layers with fields, the answer is that "An error has occurred while exucuting Python code; Null geometry cannot be converted to a point."

Can anyone help me to extracting my raster value to my original point data layer?


texts: Duplicate geometry column in geometry_columns

texts:
I have a PostgreSQL/PostGIS table with one geometry column. When I connect it to QGIS I see two of the same column:

If I make a request in PostgreSQL
SELECT * FROM geometry_columns WHERE f_table_name = 'mmune'

I can see, my view geometry_columns have two columns

I try to run function Populate_Geometry_Columns but the result is the same
SELECT Populate_Geometry_Columns('public.t_adm_commune'::regclass);

I try to redefine the column but the result is the same
ALTER TABLE mmune ALTER COLUMN geom type Geometry(Multipolygon,2) USING geom::Geometry(Multipolygon, 2);

If I question the tables pg_class and pg_attribute I see once column, that's ok.
Have you got any idea ?


texts: Database access pattern using PostgreSQL and ArcGIS

texts:
I am planning to connect to a PostgreSQL DBs on a server from ArcGIS on my laptop. My server admin is interested to know "The access pattern: How, when, how often, (would ArcGIS access the DBs) during the execution of a job. In order words, the workflow of your data analysis process." Does anyone know the details behind the execution?

Edited: Some ArcGIS tools I will be using include buffer, clip, spatial join and detect feature changes


texts: Why do I need to know what "endianness" is?

texts:
I'm just learning about GIS and shapefiles.  I am a total newb right now.  I'm probably going to get into some SQL Spatial stuff and I'm just reading about what shapefiles are.

I keep on seeing a lot about endianness.  I realize it has something to do with byte order.  My question is why do I need to know this information?

I'm looking for some context.


texts: postgresql trigger with FOR EACH ROW affect nonaffected rows

texts:
Im making a INSERT/UPDATE trigger to populate a geometry column from a text column containing a WKT string. But now its overwriting all rows!

CREATE TABLE public.teams
(
  name text,
  coverage_wkt text,
  coverage_geom geometry(Polygon,4326)
)

CREATE TRIGGER trigger_insert_or_update_on_teams
AFTER INSERT OR UPDATE ON teams
FOR EACH ROW
WHEN (pg_trigger_depth() = 0)
EXECUTE PROCEDURE set_geom_columns();

CREATE FUNCTION set_geom_columns() RETURNS trigger AS $set_geom_columns$
BEGIN
    IF NEW.coverage_wkt IS NOT NULL THEN
        NEW.name := 'im changed forever';
        NEW.coverage_geom := ST_GeomFromText(NEW.coverage_wkt,4326);
    END IF;
    RETURN NEW;
END; $set_geom_columns$ LANGUAGE plpgsql;

INSERT INTO teams (coverage_wkt) VALUES ('POLYGON((17 59, 19 59, 19 60, 18 60, 17 59))');
INSERT INTO teams (coverage_wkt) VALUES ('POLYGON((18 59, 19 59, 19 60, 18 60, 18 59))');

Running the commands above leads to:

 name |                 coverage_wkt                 |   coverage_geom
      | POLYGON((17 59, 19 59, 19 60, 18 60, 17 59)) | 010300002...
      | POLYGON((18 59, 19 59, 19 60, 18 60, 18 59)) | 010300002...

The populated geometry column holds identical geometries: the last one inserted!    That's leading me to believe its affecting all rows each time it is triggered. But from the docs I read:

  A trigger that is marked FOR EACH ROW is called once for every row that the operation modifies

Furthermore, while the coverage_geom is actually persistent after "RETURN NEW" I cannot set NEW.name the same way. I don't know if it is related.


texts: Cannot set input into parameter Coordinate_System using arcpy.AddGeometryAttributes

texts:
I have a feature class of points: AL_CB_Parcel_Intersect. This feature class is already in GCS_WGS_1984. I am trying to use arcpy.AddGeometryAttributes to create X, Y attributes, but I keep getting this error:

  ExecuteError: ERROR 000622: Failed to execute (Add Geometry
  Attributes). Parameters are not valid. ERROR 000628: Cannot set input
  into parameter Coordinate_System.

This works in ArcMap.

What am I doing wrong in my script?

import arcpy
import os
parcels_dir = r"...\Processed"
stateabb = {'Alabama': 'AL',....}

coord_sys = "GCS_WGS_1984"

for k, v in stateabb.items():
if k == 'Alabama':
    print(k)
    parcel_gdb = "{}.gdb".format(k)
    parcel_cent = "{}_CB_Parcel_Intersect".format(v)
    parcel_path = os.path.join(parcels_dir, parcel_gdb)
    parcel_feat = os.path.join(parcel_path, parcel_cent)
    arcpy.env.workspace = parcel_path
    arcpy.AddGeometryAttributes_management(parcel_feat, "POINT_X_Y_Z_M", "", "", coord_sys)
    print("{} is done!".format(parcel_cent))


texts: Volume calculation analysis with Quick Terrain Modeler 8.0.2. (QTM)

texts:
I have been unsuccessful in finding a tutorial, guide or example of the 'Volume Calculation' analysis tool in Quick Terrain Modeler.
http://appliedimagery.com/

I was wondering if anyone else has any success stories they would wish to share and/or how (in)/accurate they have found it to be when comparing (for example) a building while empty versus a building while occupied. How accurate were the measurement(s) between an empty structure and one that is filled, and how long does the calculation generally take (obviously it depends on the size though)? Unfortunately, I'm currently having issues with my LiDAR datasets, so I have yet to test it out myself, nor do I currently have a 'bare' dataset in which to reference to one that is occupied/filled.

Is anyone aware of better (or free) software that also does volume calculation analysis?

I'm using Quick Terrain Modeler v8.0.2. (x64) Build 80777.


texts: CartoDB map does not react to tap on Chrome Android

texts:
I Created a map with CartoDB

When you click the polygon, an info window should appear.
It works fine on desktop browsers, it also works fine on my phone with Firefox.
It doesn't seem to work on my phone's Chrome.

Any ideas?

My phone is Nexus 4, Android version 4.4.4
Chrome version 37.0.2062.117


texts: I have a point shapefile with several fields. Among these fields is a comment field. I would like to be able to zoom into each point and add a comment into this field in real time without having to export the attribute table, etc. Is there any way to do this in ArcMap?


texts: select street intersection with 3 oor more legs and at least one suitable for cycling

texts:
I have a bicycle and road network and I need to count the numbers of intersections with 3 or more legs and where at least one is suitable for cycling according a parameter that I assigned to the all network.
Any suggestions?


texts:
having a binary raster map of pixels of interest,
a road network,
and a target point,
my goal is to get a raster map with network pixels having as value the number of pixels of interest they allow to reach.

Would you know how i could use GRASS GIS tools to perform this?

The steps i foresee are:

rasterize the network using v.to.rast
compute pixels of interest distance to the network using r.grow.distance,
compute a direction map with r.cost,
add the target point to the network with v.net operation=connect,
alter network lines directions toward the target,
rasterize the network again using v.to.rast use=dir,
overwrite the direction map with network direction pixels with r.mapcalc,
build a directed graph from direction map,
set a value of 1 to nodes that are on pixels of interest,
perform a graph traversal toward the target, cumulating nodes values.

By now, i still don't know how to revert lines directions or build the graph.


texts: How to join lines in QGIS for later using "lines to poly"?

texts:
After successfully merging lines (using snapping) with merge_selected_features i got nice result - one merged line. Next thing was to make polygon from created merged line. Result is not nice - it creates multiple polygons connecting ends from separate lines which were merged previously together. Somehow merging geometry is not working for me or i am missing something.

..
Solved. I imported lines in Grass, then i used v.build.polylines and got nice line from which i was able to create correct polygon.


texts:
We are a land survey company and are considering using QField to perform asset surveys for local government agencies.

Looking through the documentation it appears that QField can do most of what we require. However, it is critically important that we are able to record a z-value for the surveyed position of an asset/data point. So some questions:

Does QField have the ability to record the z-value as well as the x and y values from the GPS position?
If so, is it able to utilise a geoid file in order to determine the orthometric (rather than ellipsoidal) value for the recorded z-coordinate or is this something that would require post-processing?
Given that the instrument we would be using would be attached to a survey pole, is it possible to apply a 'pole-height' in order to determine the actual elevation of the point at ground level?

This may be going beyond the scope of QField but the product we are currently using does all of these. However, it requires an always-on internet connection to function and some locations we work in do not allow for this.


texts: Access the new QGIS 3D library with pyqgis

texts:
I just found the new 3D map feature in QGIS 3, and information about the API in C++ here.

I would like to access to the API using python, for example to change camera view point and save 3D view to image file. Is there any way to do it?


texts:
What are the correct world bounds in the Spherical Mercator projection (EPSG:3857 and EPSG:900913)

In some source I got those bounds see this link in the maxExtent parameter:

minx: -20037508.3427892
miny: -20037508.3427892
maxx: 20037508.3427892
maxy: 20037508.3427892

However in FME, I got those bounds:

minx: -40075016.6855785
miny: -20037508.3427892
maxx: 40075016.6855785
maxy: 20037508.3427892


texts: How to visualize texture along borders of multipart polygon

texts:
I came across this map and would like to replicate the textured styling found along the borders.  I'm starting with a multipart polygon but my process has perhaps too many steps. I'm using QGIS 2.8.2.

Processing algorithm: Convert original polygon to lines
Processing algorithm: Fixed distance buffer
Clip buffer using the original multipart polygon
Style the clipped buffer with Line Pattern Fill, adjusting transparency, line width, etc.

Here is my output:

Any ideas on how to achieve this effect in fewer steps?  I sense I'm overlooking an option like color burst, but where you can use a line pattern as a fill.  Or, perhaps there's a way to do an inverted or negative buffer, which would save me the clipping step.

Update:
Buffers with negative values create smaller versions of the original multipart polygons, not a sort of boundary area.  (Note: I'm using Vector/Geoprocessing/Create Buffer.)

I think the buffer layer itself does need to be clipped too, presuming you have outermost boundaries to consider.


texts: Selecting feature types when using ogr2ogr to convert to shapefile?

texts:
Shapefiles can only have one geometry type, but many other formats (kml, geojson) can have multiple types.

When converting to shapefiles, is it possible to to tell ogr2ogr to create multiple shapefiles, one for each type?

It seems like the "-where" option might be useful but the man page does not explain the syntax.

For example, here's a command that fails:

$ wget http://a841-tfpweb.nyc.gov/jackson-heights/wp-content/themes/tfp/kml/transit.kml
$ ogr2ogr -f "ESRI Shapefile" transit.shp transit.kml
Warning 6: Normalized/laundered field name: 'Description' to 'Descriptio'
ERROR 1: Attempt to write non-linestring (POINT) geometry to ARC type shapefile.
ERROR 1: Terminating translation prematurely after failed
translation of layer Transit (use -skipfailures to skip errors)


texts: How to extract Z data from a shapefile from points on another separate layer QGIS

texts:
I have two layers; one with elevation data (m), and the other with roads (derived from points via points to polyines). I want to run a spatial query that basically extracts the elevation data from the start and end nodes of each link. The core function sounds like a common scenario many users may face. Is there a way of doing this?

I have provided an image to help. The blue lines are roads, and the red lines is the elevation data.


texts:
I am using ArcGIS 10.3 for Desktop.

I received a feature class with attributes, as illustrated in the first table below.

I am seeking a geoprocessing tool or process that would help me reformat the feature class into the arrangement in the second table.


texts: Default expression won't write to database

texts:
Made this expression and put it in the defaults value box.:

CASE
---cardinal
WHEN ("CATCAM" = '1')THEN'2,6'
WHEN ("CATCAM" = '2')THEN'2,6,2'
WHEN ("CATCAM" = '3')THEN'6,2'
WHEN ("CATCAM" = '4')THEN'6,2,6'
---lateral
WHEN ("MARSYS" = '1' AND "CATLAM" = '1')THEN '3'
WHEN ("MARSYS" = '11' AND "CATLAM" = '1')THEN '3'
WHEN ("MARSYS" = '2' AND "CATLAM" = '2')THEN '3'
WHEN ("MARSYS" = '1' AND "CATLAM" = '2')THEN '4'
WHEN ("MARSYS" = '11' AND "CATLAM" = '2')THEN '4'
WHEN ("MARSYS" = '2' AND "CATLAM" = '1')THEN '4'
WHEN ("MARSYS"='11' AND "CATLAM"='15')THEN '3,1,3,1'
WHEN ("MARSYS"='11' AND "CATLAM"='16')THEN '4,1,4,1'
WHEN ("MARSYS"='11' AND "CATLAM"='3')THEN '3,4'
WHEN ("MARSYS"='11' AND "CATLAM"='4')THEN '4,3'
WHEN ("MARSYS"='11' AND "CATLAM"='8')THEN '3,4,3,'
WHEN ("MARSYS"='1' AND "CATLAM"='3')THEN '3,4,3'
WHEN ("MARSYS"='2' AND "CATLAM"='4')THEN '3,4,3'
WHEN ("MARSYS"='1' AND "CATLAM"='4')THEN '4,3,4'
WHEN ("MARSYS"='2' AND "CATLAM"='3')THEN '4,3,4'
---special
WHEN ("ATONTYPE"='3')THEN '6'
---safewater
WHEN ("ATONTYPE"='4')THEN '3,1,3,1'
---isolated_danger
WHEN ("ATONTYPE"='5')THEN '2,3'
---new_wreck
WHEN ("ATONTYPE"='6')THEN '5,6'
END

The preview under the expression gives the right value, but when I make a new feature the value doesn't write to the designated column, i selected the apply on update. (field stays empty)
What am I doing wrong?

I use QGIS 3.6, the layer is stored in a Geopackage, the field is a text field with enough room to write the value.


texts: How to create Spatial Reference System in QGIS?

texts: I have some Attribute and the coordinate and the coordinate reference system. I want to create a DBmanagamet form my spatial data. To do that I want to create one table with name Spa_ref_sys, but I don't know how to populate my table. The point of this is that I want to connect my spatial data with QGIS.  How to do that?


texts: "System not initialized. Programmer forgot to call G_gisinit()." message when loading table into GRASS via QGis

texts:
I am trying to load a CSV table into Grass with v.in.asciifrom the interface of Qgis 2.0.1 on Windows XP. As the progress bar indicated 74% of the data already loaded, I got the message "System not initialized. Programmer forgot to call G_gisinit().". From that moment on, the loading stops. There is no further error message, the progress bar just gets stuck at 74%.

The exact command was:
v.in.ascii input=N:/my_path_here/AREA_NOAS04_17_131004.csv output=t7j format=point "fs=;" skip=1 x=1 y=2 z=9 cat=0

I also tried loading the CSV table into QGis, then importing the resulting vector layer (points) into GRASS with v.in.ogr.qgis, with the exact same problem (same message, loading stops at 74%).

If that helps, I'm using the CSV file contained in the Zip archive at the bottom of this page.

I'm able to load other vector layers from Qgis into GRASS without problems. Might the large size of the dataset (4'128'498 points) be a cause of the problem? I tried to google the error message, but without any helpful results.

Does anyone know what this message means and how I can solve this problem? Thanks a lot!


texts: convert geotiff to simple xyz elevation file?

texts:
I have looked everywhere and can not find a tool to convert a geotiff file to a simple ascii elevation file... did I miss something?

thanks for any pointers!


texts: Google Earth historical imagery overlayed

texts:
Over the weekend something appears to have changed with the way Google Earth's historical imagery is displayed.

Previously if I set the slider to e.g. 1970 then most of the map would be blurry brown and as I moved the slider forward in time very obvious blocks would appear over areas where there was a satellite image.

Now some sort of overlay seems to appear over the top of the actual satellite images, completely obscuring them.  This overlay will disappear and reveal the satellite images if you zoom in far enough, seems to be around 50% of the zoom bar.

There is one overlay per year and the date on all of them is December 31st of whatever year they relate to.  The copyright on the overlays all seem to be "Landsat/Copernicus", which is different to the actual satellite images.  The slider bar for changing the date of the image/s you are seeing also changes as you zoom in/out to reflect either the satellite images or the overlays.

Below is an image showing what I used to see, blocky satellite images with brown fuzzy areas where there was no image to show (top right and mid left):

And here is an image where I have zoomed out 1 step.  The blocky images are gone to be replaced by some kind of overlay:

This change basically makes it impossible for me to continue working as I have been.  Does anyone know what has caused it (Google Earth updating itself?) and how I can disable it?  I have looked through all the options in Google Earth and tried googling the problem too but with no luck.


texts: Constructing LS factor for RUSLE without outliers using ArcGIS for Desktop?

texts:
I'm using the following formula in RUSLE via raster calculator to obtain a LS factor

Power(Flowacc)*1/22.1,.4)*Power(sin(slope)*.01745/.09,1.4)*1.4

My expected outcome is a raster from 0 to approximately 53.  The range of my output raster is 0 to 6,0XX.  All of the raster is 0-53 except for concentrated flow areas which have the anomaly values.

Is there a way to exclude the anomaly concentrated flow values in formula or other?


texts:
Im using QGIS 2.6 on Ubuntu 14.04. Suddenly a strange caching behaviour appeard in QGIS. When I change the plugin code and try to reload it with for example the Plugin Reloader plugin or when I reload it manually nothing happens. Even after a QGIS restart there is still the old code in the cache. At the moment the only possibility to get the new code into QGIS is to delete my plugin out of the QGIS folder, restart QGIS, put my plugin folder into python/plugins and reload the plugin via Plugin Reloader.

I tried to reload it with qgis.utils.reloadPlugin('MyPlugin') and also with

qgis.utils.unloadPlugin('MyPlugin')
True
qgis.utils.loadPlugin('MyPlugin')
True
qgis.utils.startPlugin('MyPlugin')

The functions returning a True value but nothing happens.

Any ideas what is causing this cache behaviour?


texts: Some cities refuse to georeference in CartoDB?

texts:
My dataset includes two columns, "City" and "Country." The cities are all major world cities, nothing obscure. All are georeferencing fine in CartoDB, except Rome, Italy, and Moscow, Russia. I've even tried georeferencing in a different account, starting from scratch with a fresh, manually created dataset. But those two cities just return null values and error messages. (I've also tried variations like Roma, Moscova, and Italia.)

Any idea what's going on here? Yes, I know I could type in a lat/long manually, but I'm really curious if there's a glitch here, or if I'm doing something wrong.


texts: Which Map does ArcGIS Online use for geocoding?

texts: I would like to ask which Map provider does ArcGIS use for geocoding? My addresses are unstructured (with Floor #, Block #,...) but the precision is quite high when I manually check the locations in Google Maps.


texts: Unreasonable Raster Statistics from NetCDF file in ArcGIS

texts:
I am trying to do an analysis in ArcMap 10.3.1 using chlorophyll data from Ocean Color Web. I downloaded the data as a NetCDF file and used the Make NetCDF Raster Layer tool in order to add it to my map. When I do this, I get min/max values that are  unreasonable (like 0.13 - 0.14), and furthermore, when I click around on that layer with the "identify" tool, it frequently reports cell values well outside of the maximum. These at least seem to make sense, but I was hoping to record the mean values that show up in the source tab as part of my analysis, and now I can't do that since I cannot trust those statistics. I feel like there has to be something stupid I'm missing here. Does anybody have a clue what's happening?

Here is an image showing what parameters I entered into the Make NetCDF Raster Layer Tool. Ignore the error about the file path- I just changed that because it had my name in it.

Here is an image showing the raster statistics from the source tab and the identify window, which shows a value that is clearly larger than the supposed maximum. You can see that they are the same layer as well.


texts: What is math behind entropy calculation in ArcGIS Voronoi diagram?

texts:
ArcGIS for Desktop 10.3.1, Geospatial Analyst Toolbar, Voronoi Diagrams

How is the entropy value calculated? What is the math?

I have gone through ESRI's E-learning exercise. The exercise has a data set that contains temperature for various location.

temp 1, at (x1,y1)
temp 2, at (x2,y2)
temp 3, at (x3,y3)
etc.

One temperature reading is an outlier. I can see the 'variation' in the Voronoi diagram with the outlier (70 deg), but if I change the outlier data to be similar to the cells nearby (20 deg C), I cannot understand the diagram. I have similarly changed the values of other temperature readings, and cannot interpret the Voronoi results.

With the semivariogram, once I found out the equation, it was very simple to understand the result. Hoping for the same with the entropy Voronoi method.

I found this link, but can't quite figure out how I do the calculations.
How to find Entropy in land use (Voronoi polygons)?

I have attached the two Voronoi diagrams and the data behind the one with the outlier.


texts:
I am working with qgis 2.8.2 c++ api and qt4.8.5
in my project i had to implement the pinchzoom because i use a 27" touch screen
i am looking for information on how to use the qgis lib with the qt qtouchevent of Qt

does anyone have information ?


texts:
I'm using qgis2web plugin to create a simple web map. I have a huge attribute table and I would like to see some of fields from attribute table in web popup filed.

When I set popup fields and export map, all popups are same and show a value 1. How can I fixed that problem?

Also, in QGIS I set those points as triangle (not circles) but in web all of them are circles, plugin does not recognize point style?


texts: How to compare a SHP file and a GeoJSON file?

texts: I am working on a project to convert a mass of GIS data in Shapefile format to GeoJSON. I am able to complete the conversion using ogr2ogr but I would like to be able to test the resultant GeoJSON against the originating Shapefile. Is there an equivalent of diff that would tell me if the geometries encoded in the two files are different?


texts: Making spatial view in PostGIS and adding it as layer in QGIS?

texts:
In PostGIS 2.0, I am seeking an efficient way to create a spatially enabled view that is viewable in QGIS.

I followed the context of Part 2 of @Mike Toews answer to How to create a table from existing tables (applying new schemas) in postgis DB?, and the view shows up as a point layer (in my case) in the "Add PostGIS Layers" dialog.  But after creating the view and spatially enabling it, I received the following error when trying to add the view to QGIS:

  You must select a table in order to add a layer

Can someone explain why the spatial view will not load into QGIS, and what I need to do in order to make it loadable as a layer in QGIS?

I'm currently using QGIS v1.8 from the OSGEO installer

I view is built from the geonames database, loaded into PostGIS.  The geoname table/spatial data can be selected/viewed in QGIS just fine.

When I try to select the "mt_view", QGIS sees the view in the dialog, but won't let me select it.  It gives a dotted box around the area, but does not highlight and select the view.

I just noticed that when I run my first query:

CREATE OR REPLACE VIEW data.test AS
SELECT * FROM data.geoname
WHERE admin1='MT'

the view shows up in the QGIS Add Layer Dialog, but is unselectable.  Then when I run:

SELECT Populate_Geometry_Columns('data.test'::regclass);

I get the result of "0", where according to the documentation, I think it should be returning "1".


texts: How to change the order of layers in a Leaflet map?

texts:
On this map, there are four layers, the baselayer, segments, intersections and fatalities. I'd like them to appear in that order, bottom to top. By default for some reason the intersections layer (red circles on the screenshot) always appear above the fatalities, (black skulls).

URL for the map: https://leafletmap.neocities.org/

Why is one icon layer placed above another, in any case? Is there some default mechanism at play? The intersections were added later than the fatalities -

My latest attempt has been to create panes as in https://jsfiddle.net/abrurqwn/5/
however when I do this, all layers but the skulls/ fatalities vanish. I've attempted every re ordering of the panes, z values, etc. I can think of-

js code enclosed:

var map = L.map('map',{
    center:[45.5230 , -122.6676],
    zoom: 13,
    minZoom: 13
    });

//BASEMAP

L.tileLayer ('https://{s}.tile.thunderforest.com/cycle/{z}/{x}/{y}.png?apikey=b60c02dc2f1c476983e78e9b863fc6d1',{
    attribution: 'Map data © Thunderforest'}).addTo(map);

//FATALITIES

var skull = L.icon({
    iconUrl: 'images/blackSkull.svg',
    iconSize: [38, 95],
    popupAnchor: [-10, -30],

    });

function fatalitiesLayer (feature,
layer) {
    layer.bindPopup("
Cyclist Fatality

 "+feature.properties.name+
    "
"+feature.properties.description+
    "
Date: "+feature.properties.date+
    "
Location: "+feature.properties.location+
    "
Cause: "+feature.properties.cause+
    "
");
    layer.setIcon(skull);
    }

L.geoJson(fatalities,{
    onEachFeature: fatalitiesLayer
    }).addTo(map);

//INTERSECTIONS

var intersectionsLOGO = L.icon({
    iconUrl: 'images/intersectionLOGO.svg',
    iconSize: [38, 95],
    popupAnchor: [-10, -30],
    opacity: .5,

    });

function IntersectionsLayer(feature,
layer) {
    layer.bindPopup("dangerous intersection:
"+feature.properties.name+"
");
    layer.setIcon(intersectionsLOGO);
    layer.bringToBack;
}

L.geoJson(Intersections,{
    onEachFeature: IntersectionsLayer
}).addTo(map);

//SEGMENTS
function SegmentsLayer(feature,
layer) {
    layer.bindPopup("dangerous road segment, proceed with caution");
}

function styleLines(feature) {
    return {
                color: "red",
                weight: 9,
                opacity: .5,
                dashArray: '10,7',
                lineJoin: 'round',  //miter | round | bevel
            };
}

L.geoJSON(segments,{
    onEachFeature: SegmentsLayer,
    style: styleLines
}).addTo(map);


texts: Accessing item in list of lists using Google Earth Engine

texts:
How do I access an item in a list of lists in GEE?
I know that using the method get() returns the element at the specified position in list. But the result of list.get(0) is [1, 2, 3], and I want the first element of each list in the list.
var list = ee.List([
  [1, 2, 3],
  [2, 4, 5],
  [3, 4, 5]
]);


texts: Arcpy Spatial Join- points within a distance of lines

texts:
I am trying to perform a spatial join to join attributes from a point feature layer to the attributes of a line feature layer where the points are within a specified distance of the lines (so that these attributes can later be used in an update cursor). Currently, I am receiving a '000561: Relationship invalid for selected layers' error while trying to run this piece in Visual Studio. A spatial join based on distance between point and line features appears to work without a hitch in ArcMap and ModelBuilder. I would like to avoid field mapping if possible, but maybe that is something I may need to add in.

What am I missing?

import arcpy
from arcpy import env

arcpy.env.OverwriteOutput = True

defaultGdbPath = 'C:\Topo_Check_Tess_V5.gdb'

#lines
transLayer='C:\Topo_Check_Tess_V5.gdb\Transmission'

#points
ppLayer='C:\Topo_Check_Tess_V5.gdb\PP'

#Make feature layer
arcpy.MakeFeatureLayer_management(ppLayer, "ppLayerGU")

#Make feature layer
arcpy.MakeFeatureLayer_management(transLayer, "transLayerGU")

#output join fc
pp_trans_GUJoin ='C:\Topo_Check_Tess_V5.gdb\pp_trans_GU_SpatialJoin'

#spatial join of feature layers
arcpy.SpatialJoin_analysis ("transLayerGU", "ppLayerGU",
pp_trans_GUJoin,'JOIN_ONE_TO_ONE','KEEP_ALL','WITHIN_A_DISTANCE',.000002)


texts: How to digitize building footprints with orthogonal edges?

texts:
Using QGIS, I am trying to digitize an old raster map, and because I ultimately want to load these 2d layers into sketch-up to turn into 3D models, I need everything to be square.

Here's a sample of the map I'm working with:

The green lines you can see are a vector layer of railway tracks I've already digitized using line geometries, the rest is a raster based layer that I've geo-referenced correctly and is the background behind every layer.

What I'm aiming to do is to capture the building footprints, taking my example, and using the normal capture polygon tool, something like this:

While this works sufficiently for producing a 2D vector map (EG: to display as web-map tiles or use in QGis) it's no good for exporting as KML, then importing into Sketch-up, as everything is not 100% square, and what's more it's extremely difficult to get things 100% square.

As you can see from my tools in QGis, I have the digitize rectangles, and cad-tools plug-ins installed, but these are no good for what I need to do, for example if we take the orthogonal polygon capture tool:

I can start off drawing my first line at the angle I need, but then as soon as I press CTRL to get things at a right angle and/or square to the first line, it all goes wrong because the plug-in doesn't take the rotational angle of the building into account, meaning that any additional points are true by dedicated north/east horiz/vertical, and 90 degrees to those axis.

If the building I'm trying to trace out is on those axis, then no problem, everything's perfect, but a lot of them are not, there at various angles around the complex (It's a large factory site I'm digitizing)

The perfect tool, would be similar to the origin tool in sketch up, what you do with that is place the origin start on one corner, then stretch out an X and a Y which are at 90 deg to each other, one across the front face, one down the left face, anything you then snap to those origins, or any guides derived from it are automatically at 90 deg to the other lines, but with the rotation factored in as needed as you can see in the next image:

Capturing and digitizing the polygons in sketch-up is certainly possible without using QGis first, but if I do them that way, I can't georef them, they just come out with standard world co-ordinates centered around 0,0

Digitizing them in QGis, is done so that the location in space, is correct and conforms to the OSGB36 projection.

Sketch-up also will not obey the geo referencing in the background map, when I try to import that, it anchors the bottom left at 0,0 then asks me to scale it appropriately.

The only way I've found of getting a geo-referenced map into sketch-up is if I use the capture area from Google earth tool.  However, as of SK version 8 or higher, when doing this it only imports the base Google earth Ariel photography, and not any custom overlays you might have displayed, meaning if I load my georef in, it gets ignored when I try to capture the geo-referenced location.

Sketch-up V7 and lower did not ignore custom imagery, but V7 is no longer available for download, and even if you find a copy it wont run because the auth server won't authenticate keys for it anymore... :-(

The Ultimate Question

So essentially what I'm asking is does anyone know of a plug-in for QGis that gives me similar functionality to that in sketch-up, and if not can anyone suggest a way (or workflow) that would allow me to capture the building footprints while keeping the edges orthogonal to each other, but without being locked to the horizontal and vertical as the current tools are.

Even someway of providing a snappable grid where I could at least align the points to building corners would be an improvement.

If not, then it looks like I'm going to have to learn the QGis plug-in api and write my own plugin.


texts:
I have a pair of lines, which I like to find the intersect of each of one them with each other. I know for two lines I can do the following:

from shapely.geometry import LineString

line1 = LineString([(0,0), (1,0), (1,1)])

line2 = LineString([(0,1), (1,1)])

print(line1.intersection(line2))

Let's say I have the following multilinestring, and I am looking for finding the intersection between each pair?

MULTILINESTRING ((0 0, 0 20), (20 0, 20 20), (7.49080237694725 11.97316968394073, 5.67208473076472 5.229834926470677), (19.01428612819832 3.12037280884873, 25.96681259449994 4.024633854021577), (14.6398788362281 3.119890406724053, 21.51766222234933 1.803845175041164))


texts: Create geometries from data-only tables using virtual layers

texts:
In QGIS 3.18, I have two data-only layers (csv) - the table coordinates containing lat/lon values, the table data containing attribute data. Both tables have a common unique id: 1 (data) : n (coordinates). I now want to visualize the data on the map.
What I tried: I made a table join to add the coordinates to the data table - that works. Than I created a virtual layer with this query:
select d.field_1, d.field_2, d.field_3, d.field_4, d.field_5,d.field_6, d.coordinates_lat, d.coordinates_lon, make_point(d.coordinates_lon, d.coordinates_lat) as geometry
from data as d

This works fine.
Now I would like to do table join in the virtual layer as well. How is the correct syntax to achieve that?
Screenshot: creating point geometries from table join:


texts: Changing the width of a categorized line layer with PyQGIS?

texts:
I'm trying to change the width of a categorized line layer with PyQGIS.

Here is my code :

NOX_Actuel=iface.activeLayer()
target_field = 'NOX_Actuel'
myRangeList = []
symbol = QgsSymbol.defaultSymbol(NOX_Actuel.geometryType())
symbol.setColor(QColor("#67000d"))
myRange = QgsRendererRange(15.0, 500.0, symbol, '> 15')
myRangeList.append(myRange)
myRenderer = QgsGraduatedSymbolRenderer(target_field, myRangeList)
symbol.setWidth(5)
myRenderer.setMode(QgsGraduatedSymbolRenderer.Custom)
NOX_Actuel.setRenderer(myRenderer)
NOX_Actuel.triggerRepaint()

Whatever I try, I've the following error : "AttributeError: 'QgsGraduatedSymbolRenderer' object has no attribute 'setWidth'"


texts: Parsing error : invalid syntax (line 1) from ArcGIS Raster Calculator?

texts:
I have a (Aspect)raster (32bit, floating point,overall 2.9GB) of a landscape and need to devide it with π/180, in order to get the mean aspect later.

My raster calculator (Spatial Analyst --> Map Algebra --> Raster Calculator ) says:

  Python-Syntaxfehler: Parsing error : invalid syntax (line 1)

when I click onto it to calculate it.

What does that mean and how can I work with this raser?


texts: Seeking Open Source alternatives to Google's StreetView engine?

texts:
I am looking to create an open source javascript app using either Leaflet or OpenLayers.

Is there any way to implement panoramic views (similar to Google's StreetView) into either of these, or any other open API?

My panoramic views would be custom made from personal photos so I will not need to grab panoramic images from an existing service.

I understand that it is possible to incorporate Google StreetView into non-Google apps, but is there anything out there that would give me similar functionality without using any Google API?


texts:
I have all my data and services (WMS, WFS, WMTS) in GeoServer 2.8 (workspaces, shapefiles, etc.) and I have to migrate to GeoServer 2.12.1.

What things do I have to consider to do that?

I mean files (data_dir), PostgreSQL databases, PostGIS.


texts:
I wanted to create new field in QGIS to store date and time from a datetime string field named dateposte. Here is the code I used to make it:
import datetime
layers = QgsProject.instance().mapLayers()
for layer_id, layer in layers.items():
    #creating new field
    layer.dataProvider().addAttributes( [ QgsField("date", QVariant.Date) ] )
    layer.dataProvider().addAttributes( [ QgsField("time", QVariant.Time) ] )
    layer.updateFields()
    with edit(layer):
        for feature in layer.getFeatures():
            thedatetime = datetime.datetime.strptime(feature['dateposte'], '%Y-%m-%d %H:%M:%S')
            feature.setAttribute(feature.fieldNameIndex('date'), thedatetime.date())
            feature.setAttribute(feature.fieldNameIndex('time'), thedatetime.time())
            layer.updateFeature(feature)

But I got error on this line:
feature.setAttribute(feature.fieldNameIndex('time'), thedatetime.time())

and getting error message
Erreur OGR lors la création du champ time : Can't create fields of type Time on shapefile layers.

How can I do this?


texts: Adding shapefiles to ArcGIS Pro that only contain xxxx characters

texts:
I have a vast number of shapefiles downloaded for Forest Management Unit data and I only want to add the ones that contain the characters "epso" to my map (without having to scroll through hundreds of files to select).

For example, the files download with several file types for each unit:

t1409epbo.shp
t1409epla.shp
t1409epma.shp
t1409epso.shp
t1410epbo.shp
t1410epla.shp
t1410epma.shp
t1410epso.shp
etc.

How do I add only the epso ones?


texts: Arcpy: Describing spatial reference system

texts:
I am writing the following code to make a report describing a feature class. SO far, I got their names and shape type, but now I need to get SRS/coordinate system for each of the feature class.
I think I can use

spatial_ref = arcpy.Describe(dataset).spatialReference

but what would the parameter be if I use this in my code?

def data_report(workspace):

elements = defaultdict(list)

walk = arcpy.da.Walk(workspace, topdown=True)

for dirpath, dirnames, filenames in walk:
    for dirname in dirnames:  # 'FeatureDataset'
        desc = arcpy.da.Describe(os.path.join(dirpath, dirname))
        elements[desc['dataType']].append(desc)

    for filename in filenames:
        desc = arcpy.da.Describe(os.path.join(dirpath, filename))
        elements[desc['dataType']].append(desc)

report = []

for element_type, element_list in elements.items():
    report += [f"{element_type}:"]
    for element_description in element_list:

        if element_type == "FeatureClass":
            report +=  [f"\t\t{element_description['baseName']}"] # Here I just output the path, but there's lots more info in the Describe object
            report +=  [f"\t\t\tShape type is: {element_description['shapeType']}"]


texts: Strange output while calculating raster in QGIS

texts:
I've encountered the following problem frequently when I am using the raster calculator, even though sometimes the calculator works. However, I cannot figure out the difference between the times it works and the ones it does not.

The raster layer I am using is a gridded population. I want to create a mask layer with 1 meaning that there exists a population larger than 100 in the grid and 0 otherwise. The operation I use is ("lspop2010@1" > 100) * 1.

The output range is depicted in the picture below. Basically, the calculation produces a layer with no-data values only. Does anyone know the source of this problem?


texts: CSV file projection not aligned with Basemap

texts: I uploaded a CSV file with DMS coordinate systems in it into QGIS. The points appeared on the white background. When I add a GoogleImage background in the projected coordinate system WGSPsudoMercator, the points show on a black background. When I zoom to layer on the Google basemap, I see the points off the east coast of South Africa, they are supposed to be in Arizona, USA. Last time something similar to this happened I believe I changed the basemap from the WGS84Psudeomercator to WGS84 to align it with the points and the points showed up where I wanted them. This time, when I change the basemap projection to WGS84, the basemap disappears. I am using Degrees, Minutes, Seconds from a Juno Trimble 3B GPS. In terms of my data, I have tried switching the x field to the y-field and vise versa, adding optional numbers in front of the coordinate system such as N (for north) and W (for West), adding a negative infront of the west DMS column (the one that begins with 112). I have also tried using the other coordinates (I think LatLong or UTM), but the same problems occur. Here is the first few lines of my data in a CSV file.


texts: Calculating distance to 200m depth in bathymetric raster layer using ArcGIS Desktop?

texts:
I have a bathymetry layer and I would like to be able to calculate the distance to the 200m depth in the bathymetry raster in arcmap 10.2.

How would I go about this?

I've tried the Euclidian distance tool but it did not come out as expected and there were many NA.


texts: Import ASCII (X, Y, Value) to Feature Attribute in ArcGIS

texts:
The output from my landscaping modeling program is simply way too big to be rendered by  Paraview. So I am wondering is there anyway that I can import it into ArcMap. The data is divided into several files, with one .nodeinfo file which contains x-y coordinates and other files contain data but in the same order. I can re-arrange them into a |x|y|value| format from Matlab, but I don't know how can I import this into ArcMap.

This is just like the reverse version of "Export Feature Attribute to ASCII" function in ArcMap.


texts:
I closed the window and now I can't find it again. I've tried the Editor toolbar dropdown and Create Features is not even on there (at all, not even greyed out).

Suggestions?


texts: QGIS 3.0 can't load plugin processing, no module named 'osgeo', on Manjaro Linux

texts:
I have installed QGIS3 from the AUR package, as well as all optional dependencies:
https://aur.archlinux.org/packages/qgis/

However, when launching, I receive the following error:

Couldn't load plugin 'processing'

ModuleNotFoundError: No module named 'osgeo'
Traceback (most recent call last):
  File "/usr/share/qgis/python/qgis/utils.py", line 311, in loadPlugin
    __import__(packageName)
  File "/usr/share/qgis/python/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/__init__.py", line 28, in
    from processing.tools.dataobjects import *          # NOQA
  File "/usr/share/qgis/python/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/tools/dataobjects.py", line 51, in
    from processing.algs.gdal.GdalUtils import GdalUtils
  File "/usr/share/qgis/python/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/share/qgis/python/plugins/processing/algs/gdal/GdalUtils.py", line 35, in
    from osgeo import gdal
  File "/usr/share/qgis/python/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'osgeo'

Python version: 3.6.4 (default, Jan  5 2018, 02:35:40) [GCC 7.2.1 20171224]
QGIS version: 3.0.0-Girona Girona, exported

Python Path:
/usr/share/qgis/python
/home/thiago/.local/share/QGIS/QGIS3/profiles/default/python
/home/thiago/.local/share/QGIS/QGIS3/profiles/default/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python36.zip
/usr/lib/python3.6
/usr/lib/python3.6/lib-dynload
/usr/lib/python3.6/site-packages
/home/thiago/.local/share/QGIS/QGIS3/profiles/default/python

It seems to be a path issue, but I am not sure how to solve it. Pacman shows gdal and python-gdal both as installed.


texts: Clipping multiple shapefiles using ArcPy

texts:
I am fairly new to Python.
I created a folder with 4999 bufferfiles. Specifically, there 5km radius buffer circles, and I am trying to use the clip function with the already made bufferfiles to clip them to an environmental layerfile (called NPP2).
For the life of me I cannot figure out a way to dynamically import each of the 4999 bufferfiles into the clip function, specifically the rectangle variable.
#partially how I made bufferfiles/where they're stored
randompoint = "randompoint_" + k
bufferpoint = "buffer_" + k
randompointbufferfile = bufferpoint + ".shp"
randompointshapefile = randompoint + ".shp"
shapedirectoryString = "C:\\Users\\72sno\\Documents\\RandomPython\\FixRandomPoints\\Slope\\ShapeFiles"
bufferdirectoryString = "C:\\Users\\72sno\\Documents\\RandomPython\\FixRandomPoints\\Slope\\BufferFiles"
full_path_shape = os.path.join(shapedirectoryString, randompointshapefile)
full_path_buffer = os.path.join(bufferdirectoryString, randompointbufferfile)
arcpy.analysis.Buffer(full_path_shape, full_path_buffer, "5 Kilometers", "FULL", "ROUND", "NONE", None, "GEODESIC")

#Where I need help

#Access shapefile folder?
aprx=arcpy.mp.ArcGISProject(r'C:\\Users\\72sno\\Documents\\ArcProjects\\Prairie.aprx')
map = aprx.listMaps()[0]
map.addDataFromPath("C:\\Users\\72sno\\Documents\\RandomPython\\FixRandomPoints\\Slope\\BufferFiles")

#dynamically create 4999 clips
clippoint = "clip_" + k
clipdirectoryString = "C:\\Users\\72sno\\Documents\\RandomPython\\NPPRandomPoints2\\ClipFiles"
full_path_clip = os.path.join(clipdirectoryString, clipoint)

#How to add each 4999 buffer to the clip function (rectangle)???
arcpy.management.Clip(r"C:\Users\72sno\Documents\RandomPython\NPP2.lyrx", rectangle, full_path_clip, full_path_buffer, "0", "ClippingGeometry", "MAINTAIN_EXTENT")


texts: Python, create a polygon from range of three lat/lon points

texts:
I have three lat/lon points, I need to use those in order to create a range of all possible points within the polygon that the first three create.

Test every lat/lon if it matches a condition I have (if that point creates an angle of between 50 - 150 deg with the first two points, or not).

I have my code for the condition ready.

lons = [run.gdt1[1], run.gdt2[1], run.target_area[1]]
lats = [run.gdt1[0], run.gdt2[0], run.target_area[0]]
min_lat = min(lats)
max_lat = max(lats)
self.list_a = list(arange(min_lat, max_lat, 0.01))
min_lon = min(lons)
max_lon = max(lons)
self.list_b = list(arange(min_lon, max_lon, 0.01))

if len(self.list_a) < len(self.list_b):
    self.list_b = self.list_b[: len(self.list_a)]
elif len(self.list_a) > len(self.list_b):
    self.list_a = self.list_a[: len(self.list_b)]

self.lat_lon_list = tuple(zip(self.list_a, self.list_b))
return self.lat_lon_list

What I tried to do is to create a range between the three points and check my condition (happens on another method) on every point in the new list.

My output was:

I realized my problem is that using the min/max to create the list, creates a list only between those two points (which I don't really understand because I figured that if a point is in the range of the other two it will create a polygon of the points).

So, how can I create a list that consists of all the points that lay between those three points? (I'm using jumps of 0.1, but if too many points will return it can be narrowed down)

My desired output is a list of all points in this triangle (gdt1, gdt2, target_area are lat/lon.


texts: Convert Ordnance Survery (GB) Easting/Northing data to WGS84 lat/long

texts:
This is my first venture into GIS datasets.  I've decided to start a little project which will eventually hook into the Google Maps api to overlay/plot/pin from data in my database.

I've started loading different datasets from various sources in a MySQL database.  In particlar, I'm working with Ordnance Survey data which appears to be accurate enough for what I'd like to achieve.  However, the X/Y co-ordinates in my dataset are Ordnance Survey (GB) Easting and Northing.

Ordnance Survey provide a very useful conversion spreadsheet (https://www.ordnancesurvey.co.uk/docs/support/projection-and-transformation-calculations.xls) to convert between different datums (is that the correct term? or it Projections?)  and I've attempted to reverse engineer the formulae in the spreadsheet so I can manually do these calculations.  I totally failed at this bit.

In an ideal world, I'd like a simple calculation which can be dumped straight into a (My)SQL statement which will run through the dataset held in my table, populating additionals column for both latitude and longtitude with the converted values (in WGS84) of both easting and northing co-ordinates.

Can this be achieved using pure SQL? (i.e. What is the underlying maths etc that can be transposed into a SQL statement).

Note: This question seems to have been asked a couple of times but the solution seems to be to employ the use of an external coding language such as ASP, PHP, JS etc.  I don't really have much front end coding experience - I'm a database administrator so I'm more familiar with SQL.  I'm also pre-empting that someone will ask what front end language I will eventually be using for my project as surely I must know at least a little bit in one or most of the web dev languages.  To answer that: I'll be developing the application in "Oracle APEX" (using ODBC link from Oracle to MySQL) which is essentially pure PL/SQL, no web coding experience required.


texts: Clipping lines vector layer to polygon layer in QGIS?

texts:
I have a lines layer (a contour map of mean annual precipitation in Africa) that I want to clip to a polygon layer (shape of the African continent). However every time I use the clip tool (Vector-->Geoprocessing tools-->Clip), only an empty shapefile without attributes is created.

After looking at other questions in the site, I've already addressed two potential issues: (1) I've put the layers and the project in the same Coordinate Reference System (EPSG: 4326 - WGS 84) and (2) I converted the lines layer from multipart to singleparts after seeing this solved the problem for another user.

Something else to note, and I'm not sure if it even makes a difference: the lines vector layer was converted from a raster file using Raster-->Extraction-->Contour.

Below are a couple screenshots to add context:

Project overview:

Clip tool:

I found a solution by clipping the original raster file using the Africa polygon as a mask. I then converted that clipped raster file to a contour map and it solved my problem. As @ChrisW pointed out below, there was a lot of data off the coasts and I may have not waited long enough for original file to process (i.e., 40+ hours).


texts: Annotating north arrow degree rotation in ArcMap

texts:
I am using ArcMap 10.1.

I rotated my data frame for a map, and the north arrow dutifully is rotated as well to show north. I've been told it's possible to annotate the degree rotation on the north arrow symbol itself, but I cannot find the option in map settings or symbol properties.

How can I annotate the degree rotation on a north arrow itself?


texts: Creating file geodatabases using list as filenames in ArcPy?

texts:
I'm trying to create multiple file geodatabases from a list in arcgis. When I run the script I get a box that pops up asking for a name to be input (not what I want).

How do I reference each item in the list as a name for each fDGB?

I'm fairly certain the problem is in out_name, but, since I'm a novice, there could be other issues.

My code is as follows:

import os
import sys
import arcpy

arcpy.env.workspace = "myfilepath"
arcpy.env.overwriteOutput = True

fd = []
for i in range(1980, 2101):
    fd.append(i)

out_folder_path = "myfilepath"
out_name = (".gdb" %fd)

arcpy.CreateFileGDB_management(out_folder_path, out_name, "Current")


texts:
This may be an easy question so please forgive me if it is (or if it makes no sense) as this is all very new to me.

I have 3 files that I am trying to cross reference. One file has columns category, shape area, and shape length. Another file has category and value. The third has a location id, latitude, and longitude.

Is it possible to match the values to the locations id's? What I can't figure out is if Latitude and Longitude can be cross referenced to Shape Area and Length.

There are gdb files downloaded from FEMA that contain flood zone data.  I can open them in QGIS and see the flood zone panels as well as the attribute data that shows the flood zone name.  I have a second file in Excel with geocoded addresses.  I was hoping to pull the flood zone name for each address into Excel systematically. Is this possible?


texts:
For a project I have to use four sets of JSON datas (and one of the files is bigger than 10M).

I am loading these files asynchronously using queue.js. Once it's fully loaded and displayed I am meeting two problems.

Firstly, there is no animation when I zoom in or out the window. It's jerky and laggy. Secondly when I tried to zoom with the mouse the the zoom is systematically centered in the top-left corner.

I self-hosted my code here:
http://adrienrahier.com/test/SFmap/

you can found the sources here: https://github.com/FracArt/ChallengeSF


texts:
My project requires me to control number of stores by population density in each neighborhood.
For instance, suppose I know neighborhood A has 10 stores and a population density of 100 people/sq mile meanwhile neighborhood B has 20 stores with a population density of 500 people/sq mile.
Can I simply divide number of stores by population density to compare these values? Suppose I multiply ('number of stores'/'population density') * 100,000, could I say this is number of stores per 100,000 people/sq?


texts:
Does anybody know how to export the layer symbology from QGIS into a format that can be used by with ESRI products?

I have been successful in exporting an XML but it quotes the file type as: !DOCTYPE qgis_style. Can anybody answer whether this will import into ArcGIS developers account?


texts: How should I manage PostGIS Raster data with different projections?

texts:
I have a requirement to store and manage archaeological geophysics data that is collected as a rectangular array of samples -- a raster image.

Each raster will usually 20x20 or 30x30 floating-point samples, typically sampled at 1m intervals.
A survey will consist of one or more of these images in a given location.
It is possible that two different surveys may take place in different countries, or areas that use different projections, but each survey will use one and only one projection.
They're never likely to be viewed together, each survey will usually sit by itself.
The data will only be accessed by a custom front-end, so there will be no users getting direct control of it through psql or similar.
Each sample needs to be stored as it was collected, so I can't reproject it into a common CRS such as Web Mercator because one sample could end up covering more or less area than in the original projection, and analysis will need to be performed on the data.

How should I best store the data in a PostGIS Raster database? The options I have come up with are:

Ignore SRID constraints and store all the data in one table, writing my front-end code to deal with manipulating the data in a consistent manner.
Store all the data in one table, and rewrite the SRID constraint as a compound of SRID and survey ID.
Through table inheritance, create a new table for each new SRID.
Through table inheritance, create a new table for each survey.

1 and 2 break some of the nice automated parts of PostGIS, but will be otherwise hidden in front-end code. But queries will probably take slightly longer.

3 and 4 could end up with an explosion of tables that would make it harder to manage FK constraints and so on.

Practically, the number of rasters per survey is anywhere from 1 to 100 or more, and the number of surveys is likely to run into the hundreds. But the number of distinct projections is likely to remain very low, which favours 3.


texts:
I have some routable data for testing my pgRouting functions, and everything is fine. Now I have to download some more countries from OSM. I downloaded Germany shape files, but OSM maps are not routable. OK, next step, I found osm2pgrouting Import Tool but it works only on Linux, I need Windows version.
Is there any possibility to make OSM map routable on Windows?


texts: Binary insert into Postgis geography column results in value being inserted as geometry?

texts:
When writing WKB as a value (and as a binary parameter) to a geography column, the persisted value is not geography.

Example Java is also available here:  https://github.com/ayuudee/issue-pad/blob/master/src/com/jesusthecat/im/pggeog/BinaryGeogTest.java

Here's what it does:

Creates a table with an ID and Geography(Point, 4326)
inserts a row using WKT.
inserts a row using WKB (as bytes).
Prints out id, point, and ST_SUMMARY(point).

SQL result of step #4 is:

 id |                         pt                         | st_summary
----+----------------------------------------------------+------------
  1 | 0101000020E610000009C6C1A5E3E662406BB75D68AEED40C0 | Point[GS]
  2 | 0101000020E610000009C6C1A5E3E662406BB75D68AEED40C0 | Point[S]

Example Log of the binary insert (where pt is Geography(Point, 4326))

  LOG:  execute : insert into px(pt) values($1)

  DETAIL:  parameters: $1 = '\x0101000020e610000009c6c1a5e3e662406bb75d68aeed40c0'

You'll notice that:

The WKB for both #1 and #2 are the same; and
The flags in the ST_Summary result for #1 are [GS], where for #2 (the binary) they're [S].

The documentation for ST_Summary would indicate that #2 has spatial information, but is not geodetic (i.e. not geography).

I'm writing a Java library that seeks to persist geography as binary, but this would seem to indicate that it isn't possible. Also, is it normal that this should be allowed to happen anyway (i.e. writing a value to a geography column that's not geography)?

Versions:

  POSTGIS="2.1.3 r12547" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23
  September 2009" GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.7.8"
  LIBJSON="UNKNOWN" TOPOLOGY RASTER

  POSTGRES PostgreSQL 9.3.5 on x86_64-unknown-linux-gnu, compiled by gcc
  (Ubuntu/Linaro? 4.6.3-1ubuntu5) 4.6.3, 64-bit


texts: Leaflet & Esri-Leaflet: JSON objects not drawing correctly when using 'fields' filter

texts:
I am working with Leaflet 0.7.3 and esri-leaflet beta.6. Using L.esri.Layers.FeatureLayer works fine until I add the fields filter. At that point, the geojson objects draw sporadically (mostly not at all). I might see one or two shapes on the screen.

Example code at http://codepen.io/BBurnworth/pen/fhBrG/

Lines 32-36:

 //The line below works.
 var parcels = new L.esri.Layers.FeatureLayer(url + '24');

 // **The below line works (Checked JSON attributes in firebug) but the parcels flicker and don't display properly**
 // var parcels = new L.esri.Layers.FeatureLayer(url + '24', { fields: "CITY" });

Is this a bug or am I leaving something out? The JSON (in Firebug) reads correctly, with the geometry data coming in. It just doesn't draw correctly.

I would love for someone to point out something simple that I am missing. I would rather  fetch just the bare minimum of data and then use the tasks.IdentifyFeatures to get the attribute data when the user wants it.


texts: ArcObjects access to ArcMap's Geoprocessing Results

texts:
Is it possible to get ArcObjects access to the geoprocessing results that show up in the geoprocessing Results window in ArcMap?
I'm looking to remove results from geoprocessing operations that are run as part of a custom .Net addin.


texts: Raster statistics per latitude cross-section

texts:
I want a profile graph with biomass volume across a landscape's latitude. Im using ArcGis Pro 2.3.3. It's hard to explain so I'll provide a rough sketch of the desired result:

The sketch shows a polygon (placeholder), but I will use a raster with biomass per hectare instead. Every point on the x-axis should represent the mean of biomass per hectare for that corresponding cross-section on the map. I've been playing around with summary statistics but it won't do the job.

Do I have to divide the raster into cross-section and run summary statistics on them separately or is there a better way?


texts: Zonal stats for two polygon layers: calculate proportion of feature cover within polygons

texts:
In order to calculate land cover diversity within defined patches, I need the coverage proportion of each landcover type for each patch. So I am searching for something similar to the zonal stats tool available for raster-within-polygon calculations.

In contrast to
Calculate proportional area of polygon within another layer's polygon, I do not aim to calculate the proportional coverage of single polygons but the total proporional coverage of all polygons representing the same land cover type. Note: land cover types can be represented by more than one polygon within a patch.

The screenshot exemplarily shows two patch-polygons (dotted lines) with land cover polygons (different forest and open land types, each represented by a specific signature). For each patch polygon, I need a summary of proportional coverage of each land cover type.

The screenshot shows just a few land cover types, however there are many more - so I seek an "automatic" solution.

the result should be a "proportional coverage of landcover types"-table like this (where the "sum" culumn is not necessary):


texts: ogr2ogr extract shp-file by extent with bash script

texts:
I'm batch extracting features from shp-files in a directory with this bash script:
Is there a way to avoid the loop? I somewhere read that you can process folders but I was not able to get the right syntax..

# Author: Kay Cichini
# Date: 30-04-2013
# Purpose: Extract features by input extent for all shp-files in current dir
#          and save to new file
# Extent Tirol in EPSG:3035 - ETRS89 / ETRS-LAEA
# from D:\GIS_DataBase\GIS_Tirol\Tirol_Extent_LEAE-ETRS.shp
# xMin,yMin 4328054.73,2617730.23 : xMax,yMax 4547691.32,2739524.35

ext='4328054.73 2617730.23 4547691.32 2739524.35'
mydir=D:/GIS_DataBase/CorineLC/Vector_Data

cd $mydir

for x in *.shp
  do
  echo trying to write new_$x to $mydir
  echo extent of input file is
  ogrinfo -al $x | grep Extent
  echo new file is extracted by extent $ext
  echo ..........
  ogr2ogr new_$x -spat $ext $x
done

Here's the working example:

thanks to the links shared by @DMCI

# Author: Kay Cichini
# Date: 30-04-2013
# Purpose: Extract shp-file features in {infolder} by input extent
#          and save new shp-files to {outfolder} that is created on the fly
# Extent Tirol in EPSG:3035 - ETRS89 / ETRS-LAEA
# from D:\GIS_DataBase\GIS_Tirol\Tirol_Extent_LEAE-ETRS.shp
# xMin,yMin 4328054.73,2617730.23 : xMax,yMax 4547691.32,2739524.35

infolder=D:/GIS_DataBase/CorineLC/Single_Zips
outfolder=C:/Users/Kay/Desktop/Test_GIS/new

ext='4328054.73 2617730.23 4547691.32 2739524.35'

ogr2ogr -overwrite $outfolder -spat $ext $infolder


texts:
I am new to OpenLayers and I want to add OpenLayers.Control.WMSGetFeatureInfo to get the feature info to website, but I can't figure out right method.

I also see the source code of GeoServer OpenLayers page and still I can't find the syntax of OpenLayers.Control.WMSGetFeatureInfo code?

I use GeoServer, Apache-Tomcat-6.0.35 and OpenLayers 2.10. Should I use proxy to call OpenLayers.Control.WMSGetFeatureInfo?

This is my GeoServer WMS URL:

http://localhost:8080/geoserver/sac/wms?
service=WMS&
version=1.1.0&
request=GetMap&
layers=sac:contour&
styles=&
bbox=187754.5,204980.703125,192031.40625,210022.890625&
width=434&height=512&s
rs=EPSG:5235&
format=application/openlayers

Please give me sample code which calls OpenLayers.Control.WMSGetFeatureInfo?


texts: Digitizing overlapping hatched range maps

texts:
So I have a book of range maps for species within a taxon. As they are presented, many of them are hatched (or stippled, etc.) and overlap ("cross-hatching" shows two species coexisting in the same spot). How would I go about digitizing and extracting the range maps accurately and separate for every species?

Also could this be done with non-ESRI software, preferably R?


texts: Expression in field calculator to populate a field with the same attributes as other field

texts: I am just after a super simple expression to copy the attributes of one field directly into a new field with a different name - to export to AutoCAD with Z height data.


texts: Create MultiLineStrings segments between points with Python Shapely

texts:
I have a layer of LineString and MultiLineString geometries. I'm trying to create a single geometry of the LineStrings split by certain snapped points.
This is an illustration of the current data that I have.

I would like to create a single MultiLineString element that merges all the LineStrings between the green points. An example of a single MultiLineString (yellow line) would look like,

I'm using the following Python code to try to do this.
import geopandas as gpd

import shapely
from shapely import geometry, ops
from shapely.ops import split, snap
from shapely.geometry import Point

def split_line_by_nearest_points(line, gdf_points, tol=5):
    coords = gdf_points.geometry.unary_union
    split_line = split(line, coords)
    segments = [feature for feature in split_line]
    gdf_segments = gpd.GeoDataFrame(
        list(range(len(segments))), geometry=segments)
    gdf_segments.columns = ['index', 'geometry']

    return gdf_segments

merged_line = ops.linemerge(df_prim['geometry'].unary_union)
df_prim_split = split_line_by_nearest_points(merged_line, df_prot, tol=15)

The split line function is based on this post
Where df_prim is a geopandas dataframe of LineStrings and df_prot is a geopandas dataframe of points.
However, when I use the code described above, the geometry segments are created between each point and each vertex. The segments from this code look something like. Where the red line indicates an example segment which is split at verticies.


texts: PyQGIS getting the upper left coordinates with numpy

texts: my task is to create a QGIS plugin which loads a .tiff file and then snips the whole map in specified Raster like 1:25000 equals every 5 km. Im trying to do this with the gdal_warp function but for that I need shapefiles and i can't create shapefiles with random coordinates so I have to get upper left corner of the layer to start with.. One idea was with numpy, but i've got no idea how to load the layer in the 2D array If you guys got any other ideas im open for it..


texts: How to best fit in certain scale 500,1000,2500,5000,10000 in Data Driven ArcGis

texts:
I have more than 2000 parcels. I added this layer in Data Driven.

Need:
One of the option that data driven has about scale of the Map is "Best fit". What it does its fits parcels inside data frame in random scales.
What i need is that this "Best fits" option to be forced to certain scale 500,1000,2500,5000,10000 automatically.

For example: if parcel 1 best fit is 1:556 i need it to go to 1:1000, and if parcel 2 best fits in 1:1300 i need it to go to 1:2500

P.s.
I can use a field too to force the scaling, but this field needs to be populated manually with the scales mentioned above for every parcels, this is not good for me.


texts: How to utilise single 'randomForest' model for several rasters in R?

texts:
I am using randomForest package in R to derive species distribution models. When I predict the model to the RasterStack containing the predictors in which the values were extracted to build my model, everything goes fine. The problem is that I would like to predict this model I generated in the present scenario to the same set of predictors in the past and future. However, R gives me an error message.

Past and present predictors RasterStack are the same extension, resolution, etc. This is not the problem.

Reading the help of the 'predict' function {raster package} in R, I found the following:
"The names in the Raster object should exactly match those expected by the model. This will be the case if the same Raster object was used (via extract) to obtain the values to fit the model".

So, according to this, I can't predict my model to another RasterStack that was not used to obtain the values for the model. So, I would like to know another way that I could predict values from a present scenario to past or future scenarios.

Thanks a lot, Mariana


texts: Can anyone point me to a dataset (shapefile or other format) for the North American Continental Divide?

texts: I'm looking for a more detailed dataset than the 1:2,000,000 shapefile available from the National Atlas, preferably one suitable for use at 1:24,000.


texts: ST_Intersect is giving wrong polygon back

texts:
I imported the french IRIS dataset  by shp2pgsql.

Added LAMB93 to the PostGIS spatial_ref_sys.

shp2pgsql -W latin1 -I -s 97527:4326 CONTOURS-IRIS.shp iris20 | psql -d iris20 -v ON_ERROR_STOP=1

So far so good, when I am trying to get the polygon for the point
45.827515,4.972910  (4326)
by

SELECT insee_com, code_iris FROM iris20 WHERE ST_Intersects(geom,ST_SetSRID(st_makepoint(4.972910, 45.827515), 4326));

My only result is: 013760000

But you can see on this homepage, the coordinate is on the border of the shape:
https://www.geoportail.gouv.fr/donnees/iris

Any help?


texts:
Please somebody help because I am on the edge of madness:) I have no idea what is wrong with my code...
Only what I want is adding layers to my leaflet map: 1 point layer, 1 line, and one polygon. They all come from geoserver. here is the map: http://gis.xyz/ngrid.html
So I have added 3 layers, I can see only two, and in the layer control this two layer are swapped. It would be great if somebody look at the example, or show me a working example...


texts: DotSpatial reprojection introduces offset compared with QGIS

texts:
Im writing a converter from shape file to Geojson using C# and DotSpatial to do coordinate conversion. Currently this mostly works, with newly created shape files in QGIS the coordinates are perfect but when it comes to other files not made by me there is always an offset.

They are using this .prj

PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

That I am importing and converting coordinates like so:

private ProjectionInfo fromProjection;
private ProjectionInfo toProjection = KnownCoordinateSystems.Geographic.World.WGS1984;
private Encoding encoder = Encoding.ASCII;

public ShpPrj(byte[] prjFile)
{
    fromProjection = ProjectionInfo.FromEsriString(encoder.GetString(prjFile));
}

public ShpPoint Convert(ShpPoint point)
{
    var data = new double[] { point.X, point.Y };
    var z = new double[] { 1 };

    Reproject.ReprojectPoints(data, z, fromProjection, toProjection, 0, 1);

    return new ShpPoint { X = data[0], Y = data[1] };
}

It results in this offset:

Anyone know what could be the issue?


texts: How to generate a map showing more than one feature from the coverage layer of Atlas Generator in QGIS?

texts:
I have a shapefile (point) with occurrence records of species. There are at least two features for each species, one for currently record and a second for past record.

I need to produce one map for each species and in this map I need both types of records with different styles. I tried to follow instructions of Automating Map Creation with Print Composer Atlas at http://www.qgistutorials.com/pt_BR/docs/automating_map_creation.html
But what I got was two maps per species: one with currently records, another with pas records. I think that the problem is at the expression that I used in style rule-based: $id  =  @atlas_featureid which gives me on map per feature.

I think that I need an expression that selects features based on the field that contains the species name and show only these features on the map.


texts: Pass execution control from Context Menu BaseCommand to custom app

texts:
I got a context menu to work by following this tutorial: http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//0001000000v9000000

The interesting thing about the tutorial is that the context menu items are COM objects and they are added to the registry via the EsriRegAsm.exe tool at build time.

I add them programmatically to my context menu at runtime via their GUIDs.  I right click my map, then I display a context menu popup, then I select a menu item which is derived from a BaseCommand (here I can't debug into because it is in the COM world), now, according to the tutorial, I am supposed to do some map operations here since I have a IMapControl4 object. But what I really want to do is to pass the control up to my WPF project.

To summarize I would like to call a custom function of mine from the BaseCommand.

I am using ArcEngine 10 in VS2010 with C#.

I found a dumb hack to do it, in the OnClick method, I set the
IMapControl4instance.DocumentFilename = "SomeCommand";

Then in my main app, I have an mouse clicked event handler.  From the COM object I simulate a mouse move event from within the OnClick() via the user32.dll.  Then in my main app the mouseclick event handler checks if mapControl.DocumentFileName is one of my commands, and then sets DocumentFileName to the original value, and executes some operation.

I'm not really proud of this hack, but I'm pleased it works, I'd rather do it the normal way, whatever that is.


texts: I am trying to downscale sentinel-2 SWIR band (11) from 20 m to 10 m in ArcGIS. The create Pan-sharpened Raster Dataset tool requires colored RGB image. I only want to downscale SWIR band because red green and blue bands are already at 10 m resolution. What is the possible workaround?


texts: Which of these is the proper convention?

texts:
Possible Duplicate:
  Difference between lon and lng

I run over these three abbreviations:

"lat/lon"
"lat/long"
"lat/lng"

Are they equal?


texts: Can't view fill of small leaflet layer when zoomed out

texts:
I have a leaflet map with lots of little circles to mark particular points. When zoomed in, they look great. I can see both the orange fill and the black outline.

When zoomed out, however, I can only see the outline.

Does Leaflet provide a good way to always have both the outline and fill visible, no matter the zoom? I don't want to dynamically resize the layers on a zoom event (a la
"change width of outline for leaflet map depending on zoomlevel") because that adds a ton of complexity to my code (there are tons of these circles, zoom changes for a lot of different reasons, etc.). But I can, if that's the only way.

Here's the code where I initialize the layers

  const options = {
    color: 'black',
    fillColor: '#FF802B',
    fill: true,
    fillOpacity: 0.8,
    radius: 100,
    weight: 10,
  }
  item.layer = L.circle([lat, lng], options)


texts: Predicting nature of Result object returned by using ArcGIS geoprocessing tool in ArcPy?

texts:
In the Code Sample from the ArcGIS Pro help for the ArcPy class named Geometry there is something that catches my attention:

# Run the Copy Features tool, setting the output to the geometry object.
# geometries is returned as a list of geometry objects.
geometries = arcpy.CopyFeatures_management("c:/data/streets.shp",
                                           arcpy.Geometry())

Here the Result object (geometries) returned by running arcpy.CopyFeatures_management() is a list of Geometry objects, which is very useful for subsequent processing in that sample code.

However, other than by coming across this sample code, is there any way that I could have learned what that Result object would be?

I have reviewed what I think are the two most obvious places to look in the ArcGIS Pro help:

Copy Features
Result

I think it may be hinted at in the Result documentation where it says that the Return Value is:

  The output, either as a recordset or a string.

  If the output of the tool, such as MakeFeatureLayer is a layer,
  getOutput will return a Layer object.

  Result outputs can also be accessed by index, so result.getOutput(0)
  and result[0] are equivalent.

but I was/am hoping to gain a clearer understanding.


texts:
I have a MapInfo TAB file which looks like this:

!table
!version 300
!charset WindowsLatin1

Definition Table
  File "hp.tif"
  Type "RASTER"
  (400000,1300000) (0,0) Label "TL",
  (500000,1300000) (10000,0) Label "TR",
  (500000,1200000) (10000,10000) Label "BR",
  (400000,1200000) (0,10000) Label "BL"
  CoordSys Earth Projection 8, 79, "m", -2, 49, 0.9996012717, 400000, -100000
  Units "m"

I also have an image file for the referenced area. I'm trying to convert this to an mbtiles file which is almost working but it's slightly out of alignment.

I'm obtaining the bounding box from this file which is 400000, 1200000, 500000, 1300000 which I'm then converting to latitude / longtiude which is -2.0018997431234, 60.670742105011, -0.17172476489587, 61.580961113449.

I'm then using the following code to convert this to an MBTiles grid using zoom 12:

$x = (int)floor((($longitude + 180) / 360) * pow(2, $zoom))
$y = pow(2, $zoom) - (int)floor((1 - log(tan(deg2rad($latitude)) + 1 / cos(deg2rad($latitude))) / pi()) / 2 * pow(2, $zoom)) - 1;

This produces a grid of 2025 - 2046 and 2921 - 2943.

So far so good. Next, I scale the associated image to 256px * tiles X x 256px * tiles Y and export 256x256 tiles mapped to the grid tiles/columns.

This produces an mbtiles file which works but is slightly out of alignment with the the real map.

I have attempted to export the same files using MapTiler Desktop and this produces an mbtiles file with the same grid of 2025 - 2046 and 2921 - 2943 but it aligns perfectly. I can't for the life of me figure out why.

What is the -100000 in the CoordSys Earth Projection line for? MapTiler shows the geographical location as 400000 10 0 1300000 0 -10. I'm wondering if the -10 is an offset or something that I need to take into account.