texts: Error in Windshaft-cartodb installation step of CartoDb

texts:
I am trying to install cartodb on ubuntu 12.04. I have done all the steps mentioned on https://github.com/CartoDB/cartodb and I am stuck at the Install Windshaft-cartodb Step.
The following steps worked fine:

git clone git://github.com/CartoDB/Windshaft-cartodb.git
cd Windshaft-cartodb
git checkout master
npm install
cp config/environments/development.js.example config/environments/development.js

Whenever I hit the command node app.js development, it gives me following error:

  module.js:340
  throw err;
  ^
  Error: Cannot find module 'mapnik'
  at Function.Module._resolveFilename (module.js:338:15)
  at Function.Module._load (module.js:280:25)
  at Module.require (module.js:364:17)
  at require (module.js:380:17)
  at Object. (/home/indictrans/Documents/Windshaft-cartodb/node_modules/windshaft/lib/windshaft/server.js:8:21)
  at Module._compile (module.js:456:26)
  at Object.Module._extensions..js (module.js:474:10)
  at Module.load (module.js:356:32)
  at Function.Module._load (module.js:312:12)
  at Module.require (module.js:364:17)

I installed mapnik using sudo apt-get install libmapnik libmapnik-dev mapnik-utils python-mapnik. I also tried npm install mapnikbut the problem still persists.

I have googled lot on this issue but each time it gave the same solution.
I have also installed TileMill to map Mapnik dependecy but it is not working.


texts:
I struggle with Qgis and some drawings of water pipes for a water company. According to the program it is possible to place different signatures for different components (stopcocks, T-pieces, etc.)
From the various data located in the file "attribute table" I have to define that: the symbol is used if the value in column UNDERKLASS is equal to U22

thank you in advance
Uffe


texts: Make images/PDFs in QGIS from Vector Layers

texts:
I have a bunch of vector layers and a few polygons that I want to use to create a image out of while having my image in the background of an aerial photo.  I am unsure what commands to use to accomplish this.

Picture for reference.


texts:
I have been trying to rasterize a shapefile with polygons in it. At some point it worked but then something changed and I do not know what happened. In the code below (a bit messy, I apologize) I create a polygon of which I am sure it is within the raster that I also define. I then rasterize it but the result is a raster with just zeroes. I'm sure it is something small and stupid but I can not find it.

I have the proj4string in a function, you can find it here. Or run the following line before you run the code:

proj_string_rd_new = 'PROJCS["Amersfoort / RD New",GEOGCS["Amersfoort",DATUM["Amersfoort",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6289"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4289"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Oblique_Stereographic"],PARAMETER["latitude_of_origin",52.15616055555555],PARAMETER["central_meridian",5.38763888888889],PARAMETER["scale_factor",0.9999079],PARAMETER["false_easting",155000],PARAMETER["false_northing",463000],AUTHORITY["EPSG","28992"],AXIS["X",EAST],AXIS["Y",NORTH]]'

This is the tidy(er) code

from osgeo import gdal,ogr
import numpy as np
from PIL import Image

# open the raster layer and get its relevant properties
bbox_labs = [115733,517588,115744,517594]
xmin,ymin,xmax,ymax = bbox_labs
tileGeoTransformationParams = (xmin,0.25,0,ymax,0,-0.25)

#define a polygon which is certainly in the bounding box.
ply_wkt = '''POLYGON ((115743.63541881 517592.2735109,115742.76797302 517591.189203662,115734.852530185 517591.622926557,115734.852530185 517592.381941624,115743.63541881 517592.2735109))'''
geom = ogr.CreateGeometryFromWkt(ply_wkt)

# Create the output Driver
outDriver = ogr.GetDriverByName('ESRI Shapefile')
outDataSource = outDriver.CreateDataSource('test.shp')
outLayer = outDataSource.CreateLayer('test.shp', geom_type=ogr.wkbPolygon )
# Get the output Layer's Feature Definition
featureDefn = outLayer.GetLayerDefn()
# create a new feature
outFeature = ogr.Feature(featureDefn)
# Set new geometry
outFeature.SetGeometry(geom)
# Add new feature to output Layer
outLayer.CreateFeature(outFeature)
# dereference the feature
outFeature = None
# Save and close DataSources
outDataSource = None

# This works so I know that I can excess the shapefile and the geometry
#i = 0
#for feature in outLayer:
#    print( feature.geometry())
#    break
#print(i)

# Create the destination data source
rasterDriver = gdal.GetDriverByName('GTiff')
tempSource = rasterDriver.Create(fn_out_tiff,
                                 int((xmax-xmin)*4),
                                 int((ymax-ymin)*4),
                                 1, #missed parameter (band)
                                 gdal.GDT_Int16)

#tempSource.SetProjection(proj_string_rd_new)
tempTile = tempSource.GetRasterBand(1)
tempTile.Fill(0)
tempTile.SetNoDataValue(-9999)

# The commented lines are thing I tried but also did not work.
gdal.RasterizeLayer(tempSource,
                    [1],
                    outLayer,
#                    None,None,
                    burn_values = [1],
#                    ['ATTRIBUTE=val',
                    options = ['ALL_TOUCHED=TRUE'])
                    #options=["ATTRIBUTE=val"])

band_array = tempSource.GetRasterBand(1).ReadAsArray()
print(1 in band_array) # returns False. Should return True


texts: Snapping and topological editing in QField

texts:
Suppose I have a project with two layers, cemeteries and poly. Also, let's call cemeteries something else: hypotheses, because that makes things a bit clearer. I would like to edit these layers in QField. New features in poly cannot have overlap existing features in poly, but can overlap anything in hypotheses. New features in hypotheses can overlap other features in hypotheses, but cannot overlap existing features in poly.
Think of hypotheses as a draft layer, while poly is a definitive layer and thus has precedence, whether polygons in hypotheses are still being generated or not. Once a feature in poly has been drawn, nothing else should go there anymore, neither poly nor hypotheses.

Project Snapping Settings. Special rules for overlap. The rules are valid for all layers in edit mode. (screenshot in QGIS)

Initial situation, poly is red, hypotheses is green.  (screenshot in QGIS)

When I draw new polygons, the rules I set up above are respected.  (screenshot in QGIS)

The settings for QField are no locked geometries and both layers open for offline editing. If I try to do the same thing in QField, I get overlaps on all layers.  (screenshot in QField)

I couldn't find out how to get overlaps in hypotheses, but topologically valid geometries in poly. This way, whether or not the overlap button is active or not in QField, overlaps will always happen. If I change the Project Snapping Settings to Avoid Overlap on Active Layer, then it doesn't matter whether or not the topological-editing-button is active or not in QField, overlapping in QField will always be avoided in the active layer, instead of following the special configuration.

The active yellow feature is green and part of hypotheses. Underneath Topologische Bearbeitung eingeschaltet, you can see that the lowest corner of the new red feature has been cut off by the existing feature.  (screenshot in QField)

Conclusion: QField will respect Avoid Overlap on Active Layer, but following special rules is not yet possible. As soon as that changes, please tell me (here or wherever you can find me), as this will make my work a bit easier.

Component
Version

QGIS
3.16.5-Hannover LTR

QField Sync Plugin
3.4.2

QField
1.9.1


texts: Plotting GeoJSON data on map using ArcGIS API for JavaScript?

texts:
We are developing the portal in which we are using Arcgis javascript APIs.The requirement in our project is one of the weather provider is going to publish the webservice which will return the weather data in the geoJSON format.geoJSON will contain geometry(point,line,polygon etc)
Which is the efficient way to plot the the features on the arcGIS map?


texts: Using raster2pgsql to upload raster to PostGIS?

texts:
I am trying to upload a large number of 3 band rasters to a PostGIS database using raster2pgsql on the command line in Windows 7. Here is my command:

raster2pgsql  -I -C -s 27700 -F -t 100x100 -M F:\Filedirectory\myfile.tif mytablename | psql -h localhost -p 5942 -U postgres -d mydatabase

Sometimes (though increasingly rarely it seems) the upload works, sometimes I get the following message:

Processing 1/1: F:\Filedirectory\myfile.tif
ERROR 1: TIFFFillStrip:Integer overflow
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: F:\Filedirectory\myfile.tif , band 1:
 IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0
ERROR: rt_raster_from_gdal_dataset: Could not get data from GDAL raster
ERROR: convert_raster: Could not convert VRT dataset to PostGIS raster
ERROR: process_rasters: Could not process raster: F:\Filedirectory\myfile.tif
ERROR: Unable to process rasters

Other times I get a string of the usual INSERT 0 1 and then the same message comes up.

I have been unable to find an explanation for the meaning of this message, and I cannot find anything consistent between the files that fail and the ones that don't.

Does anyone know what this message means and how I can fix the problem?


texts: Estimate the fractal dimension of two dimensional data in R

texts:
0
down vote
favorite

I am trying to estimate the fractal dimension of a list of coordinates, such as:

          X1        X2
1  0.6073672 0.6663980
2  0.6837522 0.6080113
3  0.5109583 0.6071990
4  0.4456634 0.2165213
5  0.6773562 0.6667185
6  0.6593706 0.6774323
7  0.5778927 0.7234859
8  0.3655071 0.3496862
9  0.4694057 0.2875728
10 0.4258961 0.4424538

I was hoping to use a package like fractaldim, but it seems that I can not pass my data since it is not a square matrix.

Does anyone have an idea how I can solve this?

Background: I am trying to calculate clusters using the fractal dimension (see "Using Self-Similarity to Cluster Large Data Sets" http://dl.acm.org/citation.cfm?id=635431.635447)

EDIT1
My Data looks like this:


texts: Unable to insert data into table using geocode()

texts:
I am new to the world of PostgreSQL.
I am using PostgreSQL 9.5. I have installed PostGIS 2.2.2 plugin and along with that came Tiger, Tiger_data schemas. But the Tiger_data schema does not have any tables (not sure if it is expected). I also installed postgs_tiger_geocoder along with fuzzystrmatch. When I am trying to insert an address using geocode(), no rows are inserted. Below are the few queries I tried.

SELECT * INTO temp FROM geocode('250 WASHINGTON ST, BOSTON, MA 02108');
Select * INTO TEMP4 from geocode('343 S 500 E, Salt Lake City, UT 84102');


texts: How to find Closest Facility (euclidean distance)?

texts: I have some facilities and incidents, and I want to find the closest facility for each incident. I've done it via Network Analyst (with a network dataset) but I'm wondering if you can do it with just Euclidean distance). I'm using ArcInfo 10.0.


texts:
So when I try to export to leaflet I get this Error and I dont know how to fix it....I went in to QGIS2LEAF export and imorted all my layers and selected OSM basemap and clicked Full screen, selected the ouput folder and when I went to click okay, it stalled and gave me this error.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/RPM94/.qgis2/python/plugins\qgis2leaf\qgis2leafdialog.py", line 297, in export2leaf
    qgis2leaf_exec(self.outFileName, self.basemapName, self.basemapMeta, self.basemapAddress, self.width, self.height, self.extent, self.full_screen, self.layer_list, self.visible, self.opacity, self.encode2JSON,self.createcluster, self.webpage_name, self.webmap_head,self.webmap_subhead, self.legend,self.locate,self.address, self.precision, self.labels, self.labelshover, self.matchCRS, self.selected)
  File "C:/Users/RPM94/.qgis2/python/plugins\qgis2leaf\qgis2leaf_exec.py", line 290, in qgis2leaf_exec
    print '>> ' + crsProj4
UnboundLocalError: local variable 'crsProj4' referenced before assignment

Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.8.2-Wien Wien, 1b929ef

Python path: ['C:/PROGRA~1/QGISWI~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~1/QGISWI~1/apps/qgis/./python', u'C:/Users/RPM94/.qgis2/python', u'C:/Users/RPM94/.qgis2/python/plugins', 'C:/PROGRA~1/QGISWI~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\nose-1.3.3-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\certifi-14.05.14-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISWI~1\\bin', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISWI~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg', u'C:/Users/RPM94/.qgis2//python', 'C:\\PROGRA~1\\QGISWI~1\\apps\\qgis\\python\\plugins\\fTools\\tools', 'D:/Presentation']


texts: How to label points with two lines in tilemill

texts:
I am designing a map in Tilemill. I want to label my POIs with tow lines. If I am using

text-wrap-width:1;

then its breaking the labels with each word as a new line. Such as if the label is "One Two Three Four" then its showing like

One
Two
Three
Four

But I want to visualize like the following

One Two
Three Four

Sometimes its very hectic for viewers. Is there any way to label those with two lines always?

Thanks


texts:
I has TMS in EPSG:3857 and i want to reproject it in my local srs. My mapproxy config:

services:
  demo:
  tms:
    use_grid_names: true
    # origin for /tiles service
    origin: 'nw'
  kml:
      use_grid_names: true
  wmts:
  wms:
    srs: ['EPSG:3857','EPSG:100000']
    image_formats: ['image/png']
    md:
      title: MapProxy WMS Proxy
      abstract: This is a minimal MapProxy example.

layers:
  - name: ekb
    title: ekb
    sources: [ekb_cache_in]

caches:
  ekb_cache_in:
    grids: [webmercator]
    #disable_storage: true
    sources: [ekb_tms]
  ekb_cache:
    grids: [ekb_grid]
    meta_size: [4, 4]
    sources: [ekb_cache_in]

sources:
  ekb_tms:
    type: tile
    grid: webmercator
    transparent: true
    url: http://localhost:137/tiles/%(z)s/%(x)s/%(y)s.png
grids:
    webmercator:
        base: GLOBAL_WEBMERCATOR
    ekb_grid:
        origin: 'sw'
        tile_size: [256, 256]
        srs: 'EPSG:100000'
        bbox: [1490000, 360000, 1560000, 410000]
        #bbox_srs: 'EPSG:100000'
        res: [273.4375, 136.71875, 68.359375, 34.1796875, 17.08984375, 8.544921875, 4.2724609375,
        2.13623046875, 1.068115234375, 0.5340576171875, 0.26702880859375, 0.133514404296875, 0.0667572021484375, 0.0333786010742188, 0.0166893005371094
        , 0.0083446502685547, 0.0041723251342773, 0.0020861625671387]
        1:634 347

globals:
  srs:
    proj_data_dir: 'epsg\'

After i tried to get my reprojected layer in qgis by urls:

http://localhost:8082/tms/1.0.0/ekb/ekb_grid/{z}/{x}/{y}.png
http://localhost:8082/tiles/ekb/ekb_grid/{z}/{x}/{y}.png

i got 404 error.

I think that i have a wrong res block.
How can i generate it correctly? Or i have another mistake in config?

update
Error text:

[info] 192.168.80.85 - - [04/Apr/2018 13:45:13] "GET /tiles/ekb/ekb_grid/16/4380
3/45402.png HTTP/1.1" 404 -
[info] 192.168.80.85 - - [04/Apr/2018 13:45:13] "GET /tiles/ekb/ekb_grid/16/4380
2/45402.png HTTP/1.1" 404 -
[info] 192.168.80.85 - - [04/Apr/2018 13:45:13] "GET /tiles/ekb/ekb_grid/16/4380
4/45402.png HTTP/1.1" 404 -


texts: Georeference raster in non horizontal plane

texts:
We have raster data which is generated from (vertical) walls. We would like to keep this data in a postgis database and have the spatial reference encoded "the most accurate way".

Currently, they are saved by abusing a metric CRS and encode the z coordinate of the wall as y, and the offset from the left side of the wall as x. This gives a local reference system which works for its purpose but loses the global context.

For vector data it is straightforward to give every vertex a 3D coordinate to locate it in (global) space. This is what should be created based on the raster data (use a GIS User Interface to digitize areas of interest on top of these walls).

Furthermore, multiple walls can be situated next to each other and it should be possible to visualize these in this context (it's enough if it only works if they have the same azimuth).

There are some approaches available how this could be tackled:

Use a custom CRS in vertical space, which has its origin based on a real-world coordinate. However where exactly this "origin reference" would be stored is still unclear.

Save the information in the CRS (is that possible?) - Would require several different CRS'es for each reference plane.
Use a foreign key to a line (see red lines in sample) - Current situation, redundant information (what if the length of the line does not correspond to the width of the raster?)
Create a 3D polygon as reference plane - Redundant information, see above
Create an origin point on the line, which in combination with the azimuth of the line can be the reference plane - Would different walls share the same reference plane?

All of the approaches seem to be somehow "workarounds" and have their caveats.

The two images below show a topview of the situation and a composition of several frontal raster images. (It's ok if they are mapped to a single reference plane)

What is the most appropriate way to store the vertical raster images in the database while not losing its geographic context in horizontal space and with elevation information?


texts: Getting monthly weather means for future climates?

texts:
I am trying to generate spatially downscaled (5-10 arcminutes) future climate maps indicating monthly mean min and max temp, and precip.  Most people usually need to extrapolate monthly data into daily data (using for eg. MarkSim).  I just need the monthly.  There are various ways to do this using GCM data available on the internet, but I would like to try doing it using .cli files available from the DSSAT crop modeling program. I don't know how to get access to these files though.

Do I have to install DSSAT to get the files?


texts: Limit Cartodb map to single geographic area

texts:
This should be really simple, but I've screwed up a map several times trying to figure it out and I'm working on a deadline, so I don't have any more time to play around with it.

I just want to produce a map where the base shows only the state of Massachusetts, and none of the surrounding states. There's a FAQ link on how to do it, but I can't seem to figure it out. http://docs.cartodb.com/faqs.html#how-can-i-show-only-one-countryarearegion-on-a-map

I've applied the SQL query SELECT * FROM tablename WHERE name = 'Massachusetts'
to my map, but it doesn't seem to do anything besides screw up the filters/legend/settings I have in place.

Am I missing something critical?


texts: how to produce a RASTER representing X coordinates?

texts:
In a journal article I have read that two rasters, representing X and Y coordinates respectively (see attached figure), have been entered as predictors in a GIS-based Logistic Regression for a given reason (namely, in order to reduce possible auto-correlation issues).
I was wondering how such rasters can be produced in ArcGIS (10.1).
The reference is:

  Hu, Z., & Lo, C. P. (2007). Modeling urban growth in Atlanta using
  logistic regression. Computers, Environment and Urban Systems, 31,
  667–688. http://doi.org/10.1016/j.compenvurbsys.2006.11.001

@Luke pointed out an earlier thread. I read it, and I see that different approaches have been suggested, mainly outside arcGIS. I would rather prefer to stick with the latter. I saw something that could fit my need, but I do not understand one thing: in the suggested workaround: xmap = (FlowAccumulation(1) + 0.5)*cellsize + xmin, what that (1) stands for? The same applies to the (64) in the suggestion for the ymap (namely, ymap = (FlowAccumulation(64) + 0.5)*cellsize + xmin)


texts:
We published a project under Ubuntu Server, QGIS, LizMap web client, ..
Everything works except the "time manager".
Under the desktop client, the time manager works correctly (for example
exp: sequences between the dates 2007-01-01 to 2020-06-01) - therefore functional date format.
However, when it is placed on the server in the right place, it does not work and
the dates do not correspond (here 2021-02 -current date/month as start time- to 2021-07).
First, we worked with :

Lizmap-web-client version is 3.4.1
This is a shp.
qgis-server: Installed : 2.18.17+dfsg-1
Desktop QGIS 3.16
Fields - we tried many formats (raw data : 2007-01-01 to 2020-06-01)
Qstring as in the demo  > do not work
QDate as well yyyy-MM-dd + other date formats as yyyy-MM-dd 00:00:00  > do not work
QDate/Time yyyy-MM-ddT00:00:000. > do not work

To solve the time manager problem, we have installed on another system: ubuntu 20.04.2 LTS
qgis-server: 3.10.4 = qgis-dekstop 3.10.4
dfsg-1ubuntu2 (current version of ubuntu repository)
lizmap-web-client 3.4.1.
While running command below i got the error message below. I check, indeed the vendor directory does not exist in this version of the lizmap-web-client. Do you have an idea ?
php lizmap / install / installer.php - PHP Warning: require (/var/www/html/lizmaplive/lizmap/vendor/autoload.php): failed to open stream: No such file or directory in /var/www/html/lizmaplive/lizmap/application.init.php on line 12
PHP Fatal error: require (): Failed opening required '/var/www/html/lizmaplive/lizmap/vendor/autoload.php' (include_path = '.: / Usr / share / php') in / var / www / html /lizmaplive/lizmap/application.init.php on line 12
if not, can you confirm a qgis-server and lizmap-web-client version that works on a ubuntu server 20.04 LTS?


texts:
GEE Link here
I would like to export monthly temperature from TerraClimate, But instead of getting the 12 months period It has index of 9rows and from the year that I desire from 1995 to 2020 I got only until 2004.
This is very confusing to me as I am new to Google Earth Engine.
This also some modification from any post but I couldn't track from where.
//var geometry = ee.Point([76.50672654438921,35.35824030837368]);
var fg_points = geometry;
// Set years and month
var startYear = 1995;
var endYear = 2020;
var years = ee.List.sequence(startYear, endYear);
var months = ee.List.sequence(1,12);
// load the image collection
//var Daily = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY");
var Daily = ee.ImageCollection("IDAHO_EPSCOR/TERRACLIMATE").select('tmmx');
// make monthly summed mosaics
// loop over the years and months to get summed monthly images
var byMonth = ee.ImageCollection(ee.FeatureCollection(years.map(function(y){
  var yearCollection = Daily.filter(ee.Filter.calendarRange(y, y, 'year'));
  var byYear = ee.ImageCollection.fromImages(
    months.map(function(m) {
      var summedImage = yearCollection.filter(ee.Filter.calendarRange(m, m, 'month'))
                  .reduce(ee.Reducer.sum());
      var date = ee.Date.fromYMD(y, m, 1).format("MM_dd_YYYY");
      return summedImage.set('system:time_start', ee.Date.fromYMD(y, m, 1)).rename(date);
      //.set('month', m).set('year', y); // eventually set year and month
  }));
  return byYear;
})).flatten());
// filter the empty one out
var outputMonthly = byMonth.filter(ee.Filter.listContains('system:band_names', 'constant').not())
                    .sort('system:time_start').toBands();
print(outputMonthly);
var features = outputMonthly.reduceRegions(fg_points, ee.Reducer.first(), 30);
print(features);
Export.table.toDrive(features,
"PrecipitationHTL",
"PrecipitationHTL",
"PrecipitationHTL");


texts:
I'm trying to move from downloading geodata manually and then working with it using ogr2ogr on the local file to just reading the data in from the url directly in order to have the actual source in the code and save myself some work.

For example I'd like to grab this county boundary from the url: https://opendata.arcgis.com/datasets/4286963f248b46829ed860f4743eb85e_0.zip

this command works for me:

# ogrinfo a shapefile in a zip file on the internet:

ogrinfo -ro -al -so /vsizip//vsicurl/https://raw.githubusercontent.com/OSGeo/gdal/master/autotest/ogr/data/poly.zip

but when i change it to point at the county boundary it is unable to open the source.

ogrinfo /vsizip//vsicurl/https://opendata.arcgis.com/datasets/4286963f248b46829ed860f4743eb85e_0.zip

Do I need to supply a layer name within the zip folder?


texts: Enabling 'on the fly' CRS transformation via python console of QGIS?

texts:
I am using QGIS 2.18.

Is there a way to enable the 'on the fly' transformation project setting from the python console?


texts:
I am trying to get a full-functioning download script specifically for the Modis data (Version 6) on the USG/NASA site (AppEEARS/LPDAAC). I’m running Windows 10. I know that there is an R Modis package (RTS) available from the CRAN site and iterations on various other sites and have installed such things as wget, Rcurl, as well as the MRT reprojection tool and successfully used this manually,  however my understanding is that the USGS site has recently changed a few things including getting rid of MRT service (changing tools) and moving from FTP to HTTP which, I understand, causes issues with the R package. Python has been recommended as the way to go, with R used for any analysis required. I’ll look at ArcGIS/ R comparisons for what I want to do analytically further done the track.

I want to be able to access the NASA site, determine product availability (LST, EVI, NDVI, albedo, etc), input the area extent, mosaic required areas, select the dates/product, reproject, use scale/offset factors where required, convert from Julian to usual dates, Kelvin to Celsius and productively use the associated pixel quality data (ideally identifying cloud cover affected areas). Some of this may be better done manually, and I can use various tools in ArcGIS eg raster calculator, for scaling, etc however as I am looking at quite a bit of data (1990 to current at 5 yearly intervals) for a project, to automate this as much as possible would be the ideal way to go. There appear to be services available for some of this on the LPDAAC site.

The AppEEARS site provides API script (uCurl, R and Python) for interacting with the site and getting info and I will probably look at using Python for this. I have started using the Python API to access the site (need logon and also current token info). I have also pulled down the Modis Python toolbox (scaling and quality decoder) from the EarthData site and after some setup requirements (uses pip, arcpy, requests) have got this to run in ArcGIS, although currently not successfully in regards viable outputs re scaling and specifically getting the downloaded quality data into readable and usable format. So at the moment it is a bit bitsy with some things working and others not. Mainly accessing the Earthdata Search, Earth Explorer and AppEEARS sites for the initial searches. AppEEARS has a good “extract area sample” feature for manual, mosaic’ed, etc downloads.

My question is if anybody knows of an essentially complete python script which will allow download of required info and which I can access, similar to the R package? I’m currently looking at PyModis which appears to cover off on some requirements. I have only basic programming skills though working on it. I have looked around and asked questions on the USGS/NASA site itself but currently have had no responses.


texts: Export looped layers at once using PyQGIS?

texts:
I am trying to use one-side buffer on road centerline (layer name: small_road) by 1~40m. I've succeeded in creating one buffer (2m) file at a time using

'SEGMENTS':2,'JOIN_STYLE':0,'MITER_LIMIT':2,'OUTPUT':'C:/Users/sky/Desktop/1/buffer/b1.shp'})

Then how can I save multiple buffer files at once using loop? I've tried below code for loop, but it only creates one file named 'x.shp' with 40m buffer.

for x in range (41):
    processing.run("qgis:singlesidedbuffer", {'INPUT':'small_road','DISTANCE':x,'SIDE':0,'SEGMENTS':2,'JOIN_STYLE':0,'MITER_LIMIT':2,'OUTPUT':'C:/Users/sky/Desktop/1/buffer/x.shp'})

I also tried changing the output name into 'OUTPUT':'"C:/Path/Buffer" + str(x) + ".shp"' but it gives me error (core.QgsProcessingException: There were errors executing the algorithm).


texts: Retrieve data from table for spatial queries in PostGIS

texts:
We just want to write some spatial queries for our PostGIS tables, but we can only write queries with geometry columni like (POINT(36 43). Assume that we have two columns like geom and name, we want write our queries with the point name.

select st_astext(geom) from table
where ST_DWithin (geom, 'POINT (36 43)', 100);

Any ideas about writing these kind of queries without using geometry?


texts: ArcGIS Data Driven Pages Creates PDF in Random Order

texts: I am using ArcGIS 10.4 and Data Driven Pages to create a 130 page atlas.  The export to PDF works (exported to a single PDF), but the order of the PDFs seems to be random.  There is an ESRI shapefile that serves as the index grid, and the records in this shapefile are in numerical order.  Why does the export not export to PDF in numerical order? How do I get the PDF to build in numerical order without having to do it manually?


texts: Will ArcGIS model written in version 10.2.2 work with ArcGIS 10.1 for Desktop?

texts:
Will an ArcGIS model written in version 10.2.2 work on an instance of 10.1 desktop?

I suppose it would be somewhat dependent on the features in the model and whether they were available in the earlier version of ArcGIS.  Either way, is there some way to export the model to an earlier version to test it and see if it works?


texts: Loading shapefiles in standalone PyQGIS application without GUI

texts:
One of the related questions mentions an issue with import order, but I don't think that applies here because this is a command-line script: I'm not using Qt, only QGIS libraries.

I'm trying to run the following:

import sys
from qgis.core import QgsApplication, QgsVectorLayer

sys.path.append("/Applications/QGIS.app/Contents/Resources/python")
QgsApplication.setPrefixPath("/Applications/QGIS.app/", True)
qgs = QgsApplication([], False)
qgs.initQgis()

layer = QgsVectorLayer("/path/to/file.shp", "layer", "ogr")
layer.isValid()  # False
for field in layer.fields():
    print field  # Nothing prints

qgs.exitQgis()

Assuming the path exists and points to a valid ESRI shapefile, any ideas as to why this fails to load correctly? Does it need to be explicitly loaded into the QgsApplication context? Is there something else I'm missing?

Environment: macOS Mojave, QGIS 2.18, Python 2.7.


texts: Optimizing a nearest neighbor calculation using PostGIS

texts:
I'm using PostGIS to calculate the nearest neighbors of polygons. What I want to calculate is the minimum distance from each polygon, to the nearest polygon.

Thus far I have got great help from Mike Toews' answer (which I quote with a minor change) here:

SELECT
  a.hgt AS a_hgt,
  b.hgt AS b_hgt,
  ST_Distance(a.the_geom, b.the_geom) AS distance_between_a_and_b
FROM
  public."TestArea" AS a, public."TestArea" AS b
WHERE
  a.hgt !=  b.hgt AND ST_Distance(a.the_geom, b.the_geom) < 400

Then I calculated the minimum:

SELECT a_hgt, MIN(distance_between_a_and_b)
FROM public."lon_TestArea"
GROUP BY a_hgt

However, my challenge is to calculate this for a large number of polygons (1,000,000). As the above calculation compares each polygon to every other polygon, I wondered how I could improve the calculation so that I do not have to perform 10^12 calculations.

One thought I had was to buffer each polygon, and then to calculate the nearest neighbors of all values within the buffer for that polygon, and record the minimum. I'm not sure if that is the best approach, or whether there is a function in PostGIS that I should be using.

EDIT: Using one of Nicklas' suggestions, I'm experimenting with ST_Dwithin():

CREATE TABLE mytable_withinRange AS SELECT
  a.hgt AS a_hgt,
  b.hgt AS b_hgt,
  ST_DWithin(a.the_geom, b.the_geom, 400)
FROM
  public."lon_TestArea" AS a, public."lon_TestArea" AS b

This returns a table of the ID of each polygon, and whether it is within a certain distance or not. Is it possible to construct an IF/ELSE type statement using SQL? (I read about using the CASE condition) Or should I try joining the table I produce to the original table and then running the query again using ST_Distance?


texts: Extract Latitude/Longitude from polygon vertices on QGIS

texts:
Created a global map by adding the OpenStreetMap and WorldMap layers in QGIS. Drew polygons over that, creating what I believe is a polygon layer.

I've also created a layer of the Nodes/Vertices of the Polygon by clicking Vector->Geometry Tools->Extract Vertices (see photo)

Would like to extract the latitude and longitude of these polygon points (csv, GeoJSON, or somehow connect to MSQLServer eventually)

I've tried adding calculated fields to the Attribute Tables of the Vertices layer following but am new to GIS and have been unsuccessful.


texts: Recovering table data from table layer of ArcGIS.com

texts:
I am a researcher tracking Florida's COVID-19 tests that they report on the coronavirus dashboard here: https://experience.arcgis.com/experience/96dd742462124fa0b38ddedb9b25e429

They used to have a data navigator I could access and get the data from, but they disabled it today. I can view the table of test data in the ArcGIS viewer, but cannot figure out if I can extract/download it. I need it because it is the only place they report the age and sex breakdown of who gets tested by county.

Is there any way for me to get this data from them or is my research basically ruined?

I would normally try to do this by myself, but I don't have the luxury of time and I don't want to compromise my research by a multi-day gap.


texts: Merging separate features into single feature in shapefile programmatically?

texts:
I have a shapefile in which some features shares same values for all the attributes.

I can merge them manually by going into edit session and, selecting features to be merged, and then using merge button under Editor.

However, since number of such features are very large, a programmatic solution is required.


texts: Iterate select attribute and buffer creation

texts: I'm trying to retrace my steps from a few months ago.  I'm working on a project where I have 161 points in a shapefile.  I'm hoping there's a python script for selecting each point by the row, exporting it, and creating a 2.5 mile buffer around it.  What I meant by retracing my steps is that I've already done this in python using a script I found in here, but now I can't find it.


texts: Using OFFSETA or OFFSETB in Viewshed (ZTV) calculations in ArcGIS

texts:
I am trying to perform a Viewshed analysis in ArcGIS, and after reviewing multiple questions on this forum involving the definitions of OFFSETA v OFFSETB I am confused as to which to use.

The accepted and upvoted answer to Viewsheds and Observer Points in ArcGIS Spatial analyst appears to state that OFFSETA is the height of the target point, (in this case a turbine of 50m) and OFFSETB is the eye level height above DEM (in this case 1.65m).

Whereas the accepted and upvoted answer to Viewshed Analysis for bridge polyline observer feature appears to state that OFFSETA is the eye level height above DEM, and OFFSETB is the height of the target point.

I think the confusion is coming from the way the ESRI help documentation is worded - they appear to explain a viewshed as calculating which areas of ground you can see from an observer point, rather than which areas of ground you can see a target point from? Both are applicable, it just makes the terminology confusing if you are using the opposite approach to their example.

Their definitions are:

  There are two offset items, one defining the elevation to be added to
  the observer location and the other defining what will be added to
  each cell to be considered for visibility.

    OFFSETA
    The OFFSETA item indicates a vertical distance in surface units to be
    added to the z-value of the observation point.
    If OFFSETA exists in the feature attribute table, its value is added
    to the SPOT elevation when present; otherwise, it is added to the
    interpolated surface z-value. The OFFSETA value must be positive. If
    the OFFSETA item does not exist, the default value is 1.

    OFFSETB
    The OFFSETB item indicates a vertical distance in surface units to be
    added to the z-value of each cell as it is considered for visibility.
    If OFFSETB exists in the feature attribute table, its value is added
    to the surface z-value of each cell location when it is being analyzed
    for visibility. The value must be positive. If no OFFSETB item is
    found in the feature attribute table, it defaults to 0.*

For my own query, I have a shapefile containing target points, and a raster DEM of terrain data.

I want to see all the areas on my terrain data from which I would be able to see my points, if the points are sitting at ground height but my eye level across the terrain data is 2m.

I have used Extract Values To Points to make sure my points have ground level heights from the DEM, and I know I need to use OFFSETA or OFFSETB as an additional field in my points file to incorporate my eye level of 2m, which should add 2m to every cell of my terrain model in the visibility calc.

My understanding is I would use OFFSETB and set this to 2m, as will add 2m to each cell of my raster, rather than adding 2m to my points, which I want to keep at ground height.

Can anyone confirm or correct my understanding?


texts: How to set the cell size of a tiff to 1 meter when working in a Geographic Coordinate System?

texts:
I am working in a Geographic coordinate system in ArcGIS, and I need to create TIFFs by exporting rasters. THese .TIF file must have a resolution (cell size, or GSD) of 1 m.

Is it possible to set any parameter in the Export Raster Data window to achieve this goal?


texts: Google Earth KML Folder Point descriptions clean-up

texts:
Having trouble figuring out how to fix the Folder with my point description. The picture I added shows you 2 different point description styles. One has the main point name and popup info underneath while the other is CLEAN doesn't have the popup info underneath.

I converted both these files from ArcGis 10.2 and they both come out different. I looked at the HTML code (novice with code) and nothing seems different besides the info....

I can seem to get the code to post here.


texts: QGIS- Missing measure line box/point during editing

texts:
I've just updated my QGIS 3.4.4 to QGIS 3.10.8 and now I cannot measure distance and place point on measured distance. As soon as I click on to editing/adding point the measure tool box and line/point disappears.
I know this still works in 3.4.4 but can I somehow make it work in 3.10.8?
Attaching picture from old 3.4.4 QGIS version to better explain what I mean. As you can see on the picture there is still yellow rubberband and an end point on to which I still can snap while adding lines or points.


texts: GeoPandas: Get the maximum length inside a polygon

texts:
I'm having a a GeoDataFrame of polygons. I want to find the maximum length of a straight line that can fit inside each polygon. You could call this the characteristic length of the polygon. The length I'm looking for is shown in below image. I took the image from this SE question, but they use POSTGIS and I would like to do this with Python GeoPandas.
Is this possible with GeoPandas? If yes, could somebody explain me how to do this? Thanks in advance!
I tried polygondf.length, but this gives me the length of the perimeter, which is not what I'm looking for.


texts: Converting PDF files into ArcGIS Metadata?

texts:
I got a set of shapes with some additional informations looking like extracted/exported metadata, such as sources, extents, summary, descriptions, contacts... All in digital pdf and without XML files I could try to play around with.
Is there a way to get these information into ArcGis Metadata ? I was wondering how I could set up a metadata template with the differents types of information contained in the Pdf or if I could append the FGDC Metadata with specific field.
I will welcome any suggestions or even links to tutorial.

I would like to create an xml formular (never done that, is there some kind of wizard/tutorial for it?) with different topics such as: Name, Fullname, sources, contacts, extents, scale of production, date of creation and date of update, related legislation, "legal rights", geometry...
Then I'd like to import this formular as template for metadata and then filling the box then would be done by Hand.

I guess another way around would be to modify an existing ArcGis Metadata template. is there a tool to do that?

Here is a screenshot of the type of categories I’d like to get in a form/template. Well it's in german, but it says things like 1_General Information (sub: Editors, actualisation, contact...), 2_ sources, 3_related regulations, 4_scale, 5_ right of use, 6_technical Parameters, DCM...


texts: Fail to import gdal from conda environment?

texts:
I installed gdal from this link:
http://www.gisinternals.com/query.html?content=filelist&file=release-1911-x64-gdal-3-0-0-mapserver-7-4-0.zip
When I use the the following command from this environment and location:

**(base) C:\Users\Ibrahim Khalilullah>**from osgeo import gdal
it gives the following error:

  Traceback (most recent call last):
    File "", line 1, in
    File "C:\ProgramData\Anaconda3\lib\site-packages\osgeo__init__.py", line 21, in
      _gdal = swig_import_helper()
    File "C:\ProgramData\Anaconda3\lib\site-packages\osgeo__init__.py", line 17, in swig_import_helper
      _mod = imp.load_module('_gdal', fp, pathname, description)
    File "C:\ProgramData\Anaconda3\lib\imp.py", line 242, in load_module
      return load_dynamic(name, filename, file)
    File "C:\ProgramData\Anaconda3\lib\imp.py", line 342, in load_dynamic
      return _load(spec)
  ImportError: DLL load failed: The specified procedure could not be found.

But it successfully import if I use the following location:
**(base) C:\Program Files\GDAL>**from osgeo import gdal

How can I solve this issue? I want to use the first location for my conda environment.


texts:
I have a polygon layer and a polyline layer. I want to project the polygons onto the polyline like I show in the following images.

The first image shows the features I have.

The second image shows the projection lines (red dashed lines) as I think the concept might work.

The third image shows the resulted segments of the polyline as I expect them to be.

I don't know if this can help, the cells are from a raster grid layer and the polylines represent a road.

It is like the shadow of the polygons cast on the polyline using perpendicular lines to the polyline from the vertices of the polygon. However, the perpendicular is just a concept. I don’t know how it could better work.

I can use: QGIS version 3.10.1-A Coruña Or ArcGIS 10.4.1 for Desktop


texts: Merging overlapping polygons with original polygons in QGIS

texts:
I have a large vector shapefile dataset containing polygons, many of which overlap.

I want to assign each overlapping section to one of the original overlapping polygons. E.g. if polygon A overlaps polygon B, create a polygon C consisting of only the overlapping section of A and B and assign C to either A or B. The polygons are not continuous. Some non-overlapping polygons share boundaries and some are separated and they need to stay this way.

Is there a straightforward way of doing this in QGIS?

The data originally contained lots of invalid polygons so I started by cleaning it using 'LWGEOM Make Valid' plugin, then used the 'Saga Self-intersect' tool to create the 'C' polygons, The 'C' polygons have no attributes.

I exported the 'C' polygons into their own shapefile.
I exported the 'AB' polygons into their own shapefile.

How to merge the 'C's' into the 'AB's'?

I thought a nearest neighbour solution would work, merging or dissolving the 'C' polygons with their nearest 'AB' polygon.

I'm in the process of running the 'QGIS NN Join plugin' using 'AB' as input layer and 'C' as join layer. So far this has been running for 2 hours and has hit 2% which would consume a lot of time for 11 datasets.

I would prefer a solution in QGIS but could also access ArcGIS or try in R if necessary.


texts: Understanding CRS conversion made in LiDAR files?

texts:
I have metadata from a Lidar Point Cloud dataset and I'm having trouble understand what exactly was done to it.  It reads:

  The data (LAS files) were received in UTM Zone 18N NAD83 coordinates (meters) and vertically referenced to NAVD88 using the Geoid 12a model in meters.

This much I understand.  It is this part I'm unclear on:

  The data files (LAS files) were transformed to geographic (decimal degrees), ellipsoidal coordinates (meters).  The transform to ellipsoidal coordinates was done by reversing the application of the GEOID 12a model from a previous step.

Can someone help me out and tell me what this means?  Is the vertical spatial reference no longer NAVD88?

I'm trying to create a feature dataset and have no clue what to choose for the vertical coordinate system.


texts: How do I create Map with regions from Scratch

texts:
I wish to create choropleth map in R (or any other program).

The map I want to use cannot be found, hence I must create one from scratch (will use photo of map from 1800s). This involves

1- Drawing boundaries
2. Naming boundaries
3. Giving numerical value to boundaries

So far, I found this (https://www.qgistutorials.com/en/docs/making_a_map.html), but it doesn't answer my question.

What do I google to find this?


texts:
I'm writing a Grass Python Script (using GRASS GIS 7.6.1) for a project for university. Sadly, I never learned Python, so it's a bit tricky but so far I got everything done I wanted to.

The Script should be applicable for different regions. For the Analysis I have to reclass a DEM into 4 different classes and clip one of the classes out of my study Area. So I reclass the DEM, konvert it to vector, extract the unsuitable class and used v.overlay with the operator 'not' to keep only the study area which is not overlayed with that class.

My problem is: For the first study area I'm applying the script to, the unsuitable class is empty (as there are no such heights). But using the v.overlay with an empty binput creates an error instead of just not changing the study area and creating an output identically to the input (ERROR: No area features found in vector map . Verify 'btype' parameter.).

Is there a way to run v.overlay nevertheless? Or a way to jump that command, if one of the inputs is empty? or a way to ignore that error?
I can't just leave the v.overlay-command out of the script, as it needs to be transferable to different regions and would be needed if those heights do exist.


texts: ArcPy Insert Cursor not inserting all rows?

texts:
I have a script here that takes data from a cross-tabulated non-spatial dataset where monthly values are stored in a single table row for a given year and linearizes those data so that each monthly value is a new record in a new table. I'm using arcpy to create a search cursor to grab the cross-tabulated data and an insert cursor to write the linearized data to a new table. Here is a screenshot of the original table (obviously there are other ancillary fields as well):

The script:

import arcpy, os

months = {'JAN':1,'FEB':2,'MAR':3,'APR':4,'MAY':5,'JUN':6,'JUL':7,'AUG':8,
          'SEP':9,'OCT':10,'NOV':11,'DEC':12}
ignore = ['Arbitrary_count','TOTAL_MGY','JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEC']

ws = r'D:\Data\Users\jbellino\Project\faswam\data\water_use\SC\FromTomAb\SC_WELLS_data_jcb.mdb'
arcpy.env.workspace = ws
arcpy.env.overwriteOutput = True
tbl = 'ORIGINAL_DHEC_WELL DATA'
itbl = 'monthly_dhec_well_data'

fields = arcpy.ListFields(tbl)

rows = arcpy.SearchCursor(os.path.join(ws,tbl))
irows = arcpy.InsertCursor(os.path.join(ws,itbl))
for row in rows:
    for month in months:
        #--for each row in the original table, and for each month stored in that row,
        # create a new record in 'itbl'
        irow = irows.newRow()
        for field in fields:
            if field.name == month:
                #--if the field name refers to a month abbreviation it contains data
                # first convert the month abbreviation to month number
                irow.cn_mo = months[month]
                try:
                    # then grab the data in the field and process it into the appropriate
                    # fields of the new table in the correct units
                    irow.cn_qnty_mo_va = row.getValue(field.name)*1000000
                    irow.cn_qnty_mo_va_mega = row.getValue(field.name)
                except:
                    #--skip null values
                    pass
            elif field.name not in ignore:
                #--if the field name is not a month abbreviation, just copy the data
                # to the new table
                irow.setValue(field.name,row.getValue(field.name))
        irows.insertRow(irow)
del irows

Now, the problem is that the script will process every row in the original table, however not all rows are transferred into the new linearized table. The number of records that do make it is seemingly random and it has happened with a number of different source and destination tables. My work around has been to create/delete my insert cursor inside the for month in months: loop, but this can really bog things down with large datasets. Even then the script always "misses" the last record in the original table and I have to append the monthly values to the new table manually.

Any ideas of what may be going on here?


texts:
I currently have an ISelectionSet that has various features from different layers.

When I submit my edits to the Geodatabase, I lose my selection in ArcMap.

Is there a fast way to take the ISelectionSet that I saved and reselect that ISelectionSet?

Here is a brief overview of what my code looks like:

ISelectionSet savedSet = GetSelectionSet();

m_editor.StartEditing(GetWorkspaceBasedOnLayer(ManholeLayer));
//Make Edits To Rows Here
tempFeature.Store();
m_editor.StopEditing(true);
//Here is where the selection set is being cleared.


texts: QGIS 3.4 export DXF creating incorrect polylines

texts:
The layers are exported to dxf on QGIS. When that is imported back to QGIS or visualise on AutoCAD, incorrect Polylines can be seen.

Can you please help me to fix this?


texts: ArcGIS javascript api Rendering some parts of the Featurelayer?

texts: I'm confused how to render only some parts of a feature layer. Since all the given classbreak rendering samples are rendering the whole feature layer, is there any method to render only a few of the polygons in the feature layer (already have the Object IDs).


texts: Free dataset of world cities, including geometry

texts:
I am looking for a free dataset of world cities. I checked all the datasets suggested in the question: Seeking Country/State/City Database?, but couldn't find one where the city location is available as a polygon (or similar), I found only latitude/longitude.
Actually, I don't even need a very precise polygon, the end use is just to compute if a given point is in a city or not.


texts:
I am rewritting a custom QGIS application as a plugin for QGIS. I am testing it by writing some new code, then deactivate/activate the plugin in the plugin manager, which seems to me a bit cumbersome.
Is there a way to reload the plugin faster?


texts:
i already have an OSM server with north-america-latest.osm.pbf layer, this is rendered with mapnik.

My question is, where can I download OpenSeaMap layers and what would be the process for mounting on OSM layers that already have installed?

Would be the same process that do with the OSM layers?

Regards


texts: How erase the rest line from WGS84 projection?

texts:
I am working in world map (the Natural Earth data) which I would like to readjust for a political map of the world.

So, I have already changed the central meridian in Mercator projection, but I discover a new problem. It is a line which is the rest from orginal WGS84 projection that divide continental Russia and its island into two parts. The functions dissolve and split dont work. The line is still in the map :-(


texts:
I'm using OSM2PO (4.6.9) for "classic" routing, I mean highway network and everything works fine.

Now, I want to use it to route in a "transport" network, a bus network for instance. But I want the result to take in account bus_stops, so it will not "stop" between two bus stations.

Using PGSQL, shoud be easy to analyse source and targets coords and corresponding nodes to filter, but how can I do that in a Java custom router ?

Hope you understand my very bad english ;)


texts: QGIS (3.8 and latest master) failing to build Python plugins on Fedora 28

texts: I am building QGIS from source as I am using outdated OS (Fedora 28; upgrading is out of my hands, unfortunately).  I can get most of it built with cmake, but at around 95% it fails with segmentation fault errors when compiling the Python API.  My "pythonpath" points to Anaconda installs and is version 3.8.  I don't run into issues with CMAKE outside of it not being able to find OpenCL or GRASS to start, but after clearing the cache and manually adding those to the path it locates them but the build still fails. Have also tried checking out older versions that may more closely align with my system to no avail.  Compiler version is GNU-GCC 8.3


texts:
Is there a way to check the visibility of a label in QGIS by feature?

Could the pyqgis library help?

I need to search for points where no label is shown:


texts: Getting results from internet speed test into survey123

texts: I have a survey made in Survey123 Connect that lets people fill in some information. I would like them to run an internet speed test and then have those results (upload/download/jitter/ping/ip) populated into the feature class. How do I get that results string into the feature class?


texts: How to copy a group of layers between QGIS projects?

texts:
In ArcGIS I am used to being able to have a couple of MXD documents open and copy-paste groups of layers from one MXD's table of contents to another.  This can save a LOT of time ensuring consistent styles, but also in setting up the drawing order, pulling in the data, defining its CRS etc. for the group

I would like to do the same thing in QGIS.  I have two QGS documents open but there is no apparent option to save or copy a group of layers and dragging/dropping doesn't work (I get a little 'No Entry' symbol).

Is this possible and, if so, how?


texts:
I am brand new to mapping. I work in a hospital and I'd like to plot 55 sites of pedestrian and bike accidents onto a map of brooklyn (so we can see patterns and consider interventions). Most of my data is intersections, not addresses, plus zip.

How can I turn this into a map?

The data comes from our trauma registry, it is not public. it is in an excel sheet. i don't know what shapefiles or tabular means! it's just street addresses and interactions.

I'd like to make an online map. I got to this site from Carto and didn't realize it was separate.


texts: How to write gdal_translate output from memory directly to s3

texts:
I'm trying to convert a .tif to .mrf format using Python gdal bindings. I'm doing the conversion in memory using /vsimem then attempting to write directly to s3. The problem is only the .mrf file gets copied to s3 but none of the auxiliary side car files.
When I write to disk it writes a few side car files like .idx, lrc, .mrf.aux.xml which are necessary. It does this without me having to explicitly tell it to output those files.
How would I go about accomplishing this in s3? How do I know what files need to be copied to s3? This works when I just have a single output file.
from osgeo import gdal
import boto3

s3 = boto3.resource('s3')
object = s3.Object('test-bucket', 'testing/output.mrf')

tif_file = r"C:\new_folder\test.tif"
ds = gdal.Open(tif_file)

gdal.Translate('/vsimem/temp_file.mrf', ds, format="MRF", creationOptions=['COMPRESS=NONE', 'BLOCKSIZE=10'])

f = gdal.VSIFOpenL('/vsimem/temp_file.mrf', 'rb')
gdal.VSIFSeekL(f, 0, 2)  # seek to end
size = gdal.VSIFTellL(f)
gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
data = gdal.VSIFReadL(1, size, f)
gdal.VSIFCloseL(f)

object.put(Body=data)

ds = None


texts:
I have following two GeoDataFrames.
One with the crime count per polygon per year.

   Polygon  Year total_crime      shape_area          geometry
0        0  2009           1    3.266388e+06    POLYGON ((..))
1        0  2010           4    3.266388e+06    POLYGON ((..))
2        1  2009           2    3.266388e+06    POLYGON ((..))
3        1  2010           8    3.266388e+06    POLYGON ((..))
4        2  2009           4    3.266388e+06    POLYGON ((..))

and a GeoDataFrame where each row represents a Point of Interest.

    CATEGORY      geometry
13243   atm     POINT (..)
13244   bank    POINT (..)
13245   bar     POINT (..)
         ...
         ...
         ...

how would I join those GeoDataFrames to represent how many Points of Interests of each type are in every Polygon?

The Final GeoDataFrame should look something like this.

   Polygon  Year total_crime      shape_area    atm bank bar
0        0  2009           1    3.266388e+06      2    1   4
0        0  2010           4    3.266388e+06      2    1   4
1        1  2009           2    3.266388e+06      5    3   6
1        1  2010           8    3.266388e+06      5    3   6
2        2  2009           4    3.266388e+06      1    8   5


texts: How to get the pixels from a geotiff file in gdal python for a given polygon?

texts:
I use python and GDAL under apache to make a webservice.

I have a GeoTIFF file and i want only a subset of the data, everything inside one polygon defined by a buffer of a line (coord is the input parameter of the service):

line = ogr.Geometry(ogr.wkbLineString)
for i in range(len(coord)):
    line.AddPoint(coord[i][0], coord[i][1])

line.Transform(coordTransform)
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon = line.Buffer(10.0, 2)

I find only solutions with shapefile, not with programmatically created geometries. Any help?


texts: Is it possible to define a projection to a raster tiff map?

texts: I have a very large high resolution raster map that is a tiff file that I usually work with in Photoshop. I would like to be able to use it with ArcGIS/QGIS but am not sure how I can define the projection or if that is possible.


texts: How to build a custom raster mosaic operator?

texts:
I'm trying to find a way to build a custom raster mosaic operator. When I Mosaic different raster to a new raster ArcMap lets u chose between "first", "last", "max" etc. I need a operator that chooses not the maximum cell value but the second highest cell value. F.e. if I add 3 overlapping rasters, which contain values from 1 to 3, I want the resulting raster to be 2 instead of 3 (max) in the overlapping area. I thought this would be easily doable with a map algebra expression in the raster calculator, but can't make it work.
So I was thinking a python tool, but I can't seem to get a good start. I checked out rastertonumpyarry and this might even work. But since the raster to new mosaic tool is already almost doing what I want (except for the "second highest cell value operator") I was hoping to find an easier way.
Any ideas?


texts: Zoom to Selected Feature Keyboard Shortcut

texts:
When working with quality assessment fishnet tiles in ArcMap most of the time couple of thousands i have to double click each record to move from one fishnet to the other. Instead of that i am looking for a keyboard shortcut which can do first go to the next record line and also at the same time zoom to that feature to a fixed scale. In the ArcMap attribute table view Ctrl + Enter will take you the next row record the missing one is how to zoom to that feature in the map view with keyboard shortcut at the same time.

I do not have that much programming know how.


texts: satellite imagery basemap for webmapping application

texts:
I am currently trying to develop a small webmapping application to show some point data which covers the whole planet.
I do want to use projections different from webmercator, like Robinson projection, for example.
Additionally, I want to use a satellite imagery basemap.  The Problem is that almost all the ready-to-use (satellite imagery) map tiles are webmercator projected.
I have found (here in this forum) just one provider offering tiles in projections different from webmercator: https://earthdata.nasa.gov/about/science-system-description/eosdis-components/global-imagery-browse-services-gibs
But the projections used there are not these that i do need.

How can i add a satellite imagery basemap with projections different from webmercator?


texts: Any difference between EPSG:3003 and EPSG:102091?

texts:
I am working on QGIS 2.2.0.

The PRS I am working on is Monte Mario / Italy zone 1 (EPSG:3003), with the following parameters:

PROJCS["Monte Mario / Italy zone 1",
GEOGCS["Monte Mario",
    DATUM["Monte_Mario",
        SPHEROID["International 1924",6378388,297,
            AUTHORITY["EPSG","7022"]],
        AUTHORITY["EPSG","6265"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4265"]],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",9],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",1500000],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","3003"],
AXIS["X",EAST],
AXIS["Y",NORTH]]

A dataset I need to use on the above PRS is in Monte_Mario_Italy_1 (EPSG:102091), with the following parameters:

PROJCS["Monte_Mario_Italy_1",
GEOGCS["GCS_Monte_Mario",
    DATUM["Monte_Mario",
        SPHEROID["International_1924",6378388,297]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",1500000],
PARAMETER["False_Northing",0],
PARAMETER["Central_Meridian",9],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0],
UNIT["Meter",1],
AUTHORITY["EPSG","102091"]]

Is there any difference between the two PRS? (I see a difference in decimals in parameter UNIT["Degree", ...]).

In other words, do I need a re-projection or datum transformation to analyse the data together (e.g. buffering, calculating distances etc) ?

Thanks in advance.
Alessandro


texts: ASCII ID,Lat,Lon,Hgt to raster min. distance in QGIS

texts: I am new to QGIS, and I would like to know how I can create a raster colour map displaying (through colour) the minimum distance between each of the pixels of the raster and a set of geographic points (lat, long) stored in an ASCII file, so for instance if a pixel is right on top of one of the ASCII points would be green, if one pixels is 30 pixels (or by geographic distance) away from a point it would be red and so on.


texts: Turning off representation descriptions in ArcMap table of contents?

texts:
Using ArcMap 10.3.1, advanced, how do I turn off the display of the layer representation descriptions in the TOC?  I've highlighted what I'm trying to turn off in the screenshot below.

Nothing in the Arc help mentions this, nor have I found anything in the Layer Properties window (also shown in the screenshot).


texts:
How to make grid annotations like shown below in QGIS (v3.x)?
Sequential alphanumeric characters in the centre of grid lines. These are widely used in tourist maps.


texts: Retrieving OpenStreetMap data (ways and its child nodes) by coordinates using Overpass API?

texts:
I am new to using OSM.I am trying to download OSM data programmatically using Python. I need data of ways and its child nodes contained in a perimeter specified by coordinates. Currently, I am referring to https://github.com/mvexel/overpass-api-python-wrapper which is mighty useful, but I am facing issues with the APIs.

This does not give me any data. Not sure if it has to do with syntax.

data = api.get('way(42.819,-73.881,42.820,-73.880)')

When I try

MapQuery = overpass.MapQuery(50.746,7.154,50.748,7.157)
response = api.get(MapQuery)

I get a list of of features where I can't tell which is way or node.

I can work with XML or CSV format. I was hoping the data would look like a .osm file downloaded from JOSM.


texts: Google Earth Engine: Export an entire collection

texts:
I would like to export all the images that the collection (code) contains. How can I do that?

var colection = ee.ImageCollection ("MODIS/NTSG/MOD16A2/105")


texts:
I'm a totally beginner with QGIS and maybe I have a very easy problem, but I can't help myself.
I got an elevation shapefile with classes of elevations like that: 1-100m, 101-200m, ...
When I try to classify this dates to display them with a continuous shift of colours, QGIS always sorts the classes like that: 1-100m, 1001-1101m, 101-200m, 1101-1200m, ...
So the colours of the elevation classes are mixed and don't display the continuous shift. Is there an easy way to sort the classes in the attribut table in the right way?
Thanks a lot!


texts: Short streams and tiny basins | QGIS-SAGA

texts:
After watching a lot of tutorials and trying all approaches I still somehow manage to get streams that are too short over my DEM, resulting in very small basins. I do it the standard way:

Get a DEM
Fill sinks with the Wang & Liu algorithm in SAGA
Get the Strahler Orders

The resulting streams are very short, resulting in countless small basins. I figured that using r.watershed I can get a better model, but then the stream segments only appear by the outlets.
I want to be able to get the basins that correspond to a specific order. Can anyone tell me what I am doing wrong, or is my topography just special?


texts: Could someone let me know how I add a API key to QGIS 2.14.10. The diagonal API key required is still on my map layer. Bit of a novice on QGIS.


texts: Converting GeoTIFF to tiles using gdal2tiles

texts:
I have a GeoTIFF of weather radar data that I cannot seem to get converted to a raster tileset using gdal2tiles. When the tileset is generated, there are TONS of empty folders in the output directory, with a few holding a single transparent tile in them. I am not sure what is causing this, as I had other GeoTIFFs before (another source) and the same gdal2tiles script worked fine. I recently switched to a new process and I have to add the projection/bounds to the GeoTIFF myself, and I think I am going wrong somewhere when I am adding data to the image.

The original clean TIFF may be downloaded here for inspection :
https://www.dropbox.com/s/2dh74pjrjgkhqe0/latest.zip?dl=0

I know the projection and bounds the image must have, so as a first step before using gdal2tiles I have tried adding them.

First, I add the projection 3857 for use with Mapbox, and add the bounds. Have I done this correctly? The North extent should be 30.8, West is -104.7, East is -103.8, and South is 29.8.

gdal_translate -a_srs EPSG:3857 -a_ullr -104.7 30.8 -103.8 29.8 MRMS_RALA_LATEST.tiff projected.tiff

Next, I use gdalwarp to warp the raster :

gdalwarp -t_srs EPSG:3857 projected.tiff projected2.tiff

Last, I colorize it :

gdaldem color-relief projected2.tiff mrms_reflectivity.txt colorized.tiff

And this is my gdal2tiles command. Again, all of the tiles it does output are blank.

gdal2tiles --profile=mercator -z 0-13 colorized.tiff outputfolder

For what it is worth, I have attached a screenshot of these commands being executed.

Also, for completeness - this is the colorization file used with gdaldem:
https://www.dropbox.com/s/5v65j3z6orirlk0/mrms_reflectivity.txt?dl=0


texts: QGIS 3 keep creating a "palettes" useless folder

texts:
Fresh installations of QGIS 3.0.x systematically create a "palettes" folder in c:/ on Windows. Even if you delete it, it gets recreated after every restart/reinstall of QGIS.

Any trick to get rid of this folder once and for good? Or at least to change the path in which it is created?


texts: Where/How to Learn about Getting Started with Web GIS

texts:
Possible Duplicate:
  Steps to Start Web Mapping

Does anyone know of any classes or resources where I could learn about WebGIS?  I'd like to learn the steps and technology involved.  I'd be looking for something online or live in the San Francisco Bay Area.

Thanks


texts: Python and QGIS multiprocessing documentation?

texts:
Where can we find better documention for python QGIS multiprocessing.
Especially regarding the multiprocessing of python and QGIS.


texts: Fitting merged polygon within adjacent line features using QGIS

texts:
Is it possible to fit the polygon into the adjacent lines, when moving my object to the different location and can it be done automatically?

Once I moved feature, I am unsure about it location within the borderlines lying next to. I would like to make them fitted to the polygon edge.


texts: Expanding river width on raster based DEM in ArcMap

texts:
How do I expand the river width on a raster based DEM in ArcMap?
For example, my DEM has 20 m cell size, and I want that my river is wide 3x that everywhere (60 m in total, or 3 cell wide). I have river centerline too and I do not want to change the current elevations, I actually want to expand the center box elevation to adjacent 2 cells from each side of the river.


texts: QGIS create boundaries from background image

texts: I want to manually trace over an existing map image to replicate the boundaries using QGIS.  If this was an existing polygon and I need to split this into 3 separate polygons, I would take the existing outer boundary, and then I would have to manually trace the boundaries of the inner areas - as a polyline or region? How would I join the inner boundaries to the outer boundary, and then create the 3 separate polygons within the base map table? I hope I have explained this clearly?


texts: Field calculator to automatically update fields in ArcGIS Desktop?

texts:
Trying to get an attribute table to update based on the values in their corresponding rows. for example - If I have a selection of fields whose value is 3, I want the corresponding fields to auto populate for that value, in this case road surface types.

I am using ArcGIS Desktop 10.5.


texts: Populating comment field using Python Parser of ArcGIS Field Calculator?

texts:
I am new to using Python and seek help building a code for an ArcGIS project. I took one introductory C programming course in school, but other than that I am green.

My objective is to have the conditions below fill a new comment field in the attribute table with the following return statements. Multiple return statements can be true for one attribute.

def comments(?): *comments is the field I wish to populate with return statements.

    if (((PLOT_MASTER.HEIGHT - "BL_FinalLayer_SpatialJoin.OHT" ) >0.1) OR ((PLOT_MASTER.HEIGHT - "BL_FinalLayer_SpatialJoin.OHT" ) <-0.1)):
        return “Height does not match OHT”

    elif (((PLOT_MASTER.STOCKING_D - "BL_FinalLayer_SpatialJoin.STKG" ) >0.01) OR ((PLOT_MASTER.STOCKING_D - "BL_FinalLayer_SpatialJoin.STKG") <-0.01)):
        return "Stocking_D does not match STKG"

    elif (PLOT_MASTER.LEADING_SPECIES_CODE<> "BL_FinalLayer_SpatialJoin.OLEADSPC"):
        return "Leading_Species_Code does not match OLEADSPC"

    elif (PLOT_MASTER.AGE<> "BL_FinalLayer_SpatialJoin.OAGE"):
        return "Age does not match OAGE"

    else ((( PLOT_MASTER.ECO_CODE - "BL_FinalLayer_SpatialJoin.PRI_ECO_2" ) >0.01) OR ((PLOT_MASTER.ECO_CODE - "BL_FinalLayer_SpatialJoin.PRI_ECO_2") <-0.01))
        return "ECO_CODE does not match PRI_ECO"

comments(?)


texts:
I inherited some data that lines up correctly on the map, but in face are in different coordinate reference systems/SRID. Here is a screenshot example:

The CRS of the green layer is:

EPSG: 4326

The CRS of the purple layer is:

EPSG: 2264

They are lined up correctly on the map, but to do stuff in PostgreSQL/PostGIS with them, they need to be in the same CRS/SRID. If I change the CRS of the green layer to match the purple layer, naturally the green layer is now way off the map.

I could manually recreate the green layer, but this would be time consuming. Is there a better solution? Could I re-create the green layer as it is, but with a different CRS? Either using some tool in QGIS or in PostgreSQL/PostGIS?


texts: calculate distance of points from centroid

texts:
I have a spatialpoints object, and a spatialpolygonsdataframe, composed of rectangular cells, and which the spatialpoints fall into. What I would like is to determine which polygon in the spatialpolygonsdataframe each point from the spatialpoints falls into, the distance each point is from the centroid of the respective polygon, and the angle that the point is in relationship to the centroid (with the centroid acting as 0,0).

I know, for example, I could use geosphere::bearing(), raster::pointDistance(), and one of the many centroid functions, but I am not sure how to combine all this into working code.

Here's some sample code. SP is the spatialpoints object, and rintersect3 is the spatialpolygonsdataframe.

owin2Polygons <- function(x, id="1") {
  stopifnot(spatstat::is.owin(x))
  x <- spatstat::as.polygonal(x)
  closering <- function(df) { df[c(seq(nrow(df)), 1), ] }
  pieces <- lapply(x$bdry,
                   function(p) {
                     sp::Polygon(coords=closering(cbind(p$x,p$y)),
                             hole=spatstat.utils::is.hole.xypolygon(p))  })
  z <- sp::Polygons(pieces, id)
  return(z)
}

owin2SP <- function(x) {
  stopifnot(spatstat::is.owin(x))
  y <- owin2Polygons(x)
  z <- sp::SpatialPolygons(list(y))
  return(z)
}

test.verts <- list()
    test.verts$x <- c(518, 594, 421, 286, 104, 129, 315)
    test.verts$y <- c(273, 608, 1202, 1201,  620,  285,  135)
    test.owin <- spatstat::owin(poly = test.verts)
    ppp <- spatstat::ppp(x = c(344, 276, 535, 506, 459), y = c(194, 865, 551, 272, 438), window = test.owin)

    SP <- maptools::as.SpatialPoints.ppp(ppp)

    w3 <-  spatstat::Window(ppp)

    #below converts window, and overlays a grid with 10x10 cells
    spatpoly3 <- owin2SP(w3) #internal function

    bbow3 <- spatstat::as.owin(spatstat::boundingbox(w3))
    spatpoly3 <- sf::st_as_sf(spatpoly3)

    #---Alphanumeric Grid:-------
    # make grid
    ng1 <- sf::st_make_grid(spatpoly3, n = c(5, 5))
    ng1 <- sf::st_sf('ID' = seq(length(ng1)), ng1)

    # rasterise grid (ID becomes cell value)
    nr1 <- fasterize::raster(ng1, ncol = 5, nrow = 5)
    nr1[] <- ng1$ID

    # get X and Y from row and col numbers, convert Y to letters

    ng1$X <- raster::colFromCell(nr1, ng1$ID) # or seq(ncell(nr1))
    ng1$Y <- rev(sapply(raster::rowFromCell(nr1, ng1$ID), excel_col))

    # combine for label
    ng1$LAB <- paste0(ng1$Y, ng1$X)

    # trim back to polygon
    clip2 <- sf::st_join(lwgeom::st_make_valid(ng1), sf::st_as_sf(spatpoly3), join = st_intersects)

    geo <- as(clip2, "Spatial")
    spatpoly3 <- owin2SP(w3) #internal function
    spatpoly3 <- rgeos::gBuffer(spatpoly3, byid=TRUE, width=0)
    geo <- rgeos::gBuffer(geo, byid=TRUE, width=0)

    CRS.new <- sp::CRS("+init=epsg:27700")

    #proj4string(geo) <- CRS.new
    sp::proj4string(spatpoly3) <- CRS.new

    geo <- sp::spTransform(geo, sp::CRS("+init=epsg:27700"))

    #Finished product; leaf with alphanumerically labelled grid
    rintersect3 <- raster::intersect(geo, spatpoly3)

    sp::proj4string(SP) <- CRS.new

    points_over3 <- GISTools::poly.counts(SP, rintersect3)


texts: Does not accept in raster calculator variable s the number QGIS

texts:
I'm writing a plugin, where self.dlg.lineEdit_11.text() This number (attached image)

gives an error message

  Traceback (most recent call last): calc = QgsRasterCalculator("rastr_raznica1@1" + s, var, 'GTiff', layer2.extent(), layer2.width(), layer2.height(), entries ) TypeError: cannot concatenate 'str' and 'float' objects

through processing, I appeal to the raster calculator, but he does not want to be considered QgsRasterCalculator("rastr_raznica1@1" + s) It does not take s the number

def edit_raster(self):
    entries =[]
    rastrway1=self.dlg.lineEdit_11.text()
    var= os.path.join(os.path.abspath(os.path.dirname(rastrway1)),'rastr_sko.tif')
    tif2=  os.path.join(os.path.abspath(os.path.dirname(rastrway1)),'rastr_raznica1.tif')
    r = self.dlg.lineEdit_10.text()
    s=float(r)
    print s
    layer2 = QgsRasterLayer(tif2, "rastr_sko")
    her_2 = QgsRasterCalculatorEntry()
    her_2.ref = 'rastr_raznica1@1'
    her_2.raster = layer2
    her_2.bandNumber = 1
    entries.append(her_2)
    calc = QgsRasterCalculator("rastr_raznica1@1" + s, var, 'GTiff', layer2.extent(), layer2.width(), layer2.height(), entries )
    calc.processCalculation()


texts: the layer imported from csv I'm trying to add in qgis is in the middle of the ocean

texts:
I'm trying to add sampling points to a layer which was initially created.
This layer is in Bangalore and the points I'm trying to add is in the middle of the ocean, but if I search for the points on the maps generally, the coordinates seem to be right.
I've also tried changing the lat and long for x and y axis respectively.
I've checked for the CRS systems also which are set to the same system.
Please help me identify the issue cause I fail to see where the error is.


texts: Simple way to convert GeoJSON from MultiPolygon into Polygon

texts:
I need a GeoJSON file that simply has the coordinates of the polygons for the boundaries of the provinces of Canada.
I have a GeoJSON file but it has multipolygons and my software does not support multipolygons. How can I convert the multipolygons to simple polygons?
Of course, if the geometry of the province consists of disconnected regions with pieces of other provinces in the interior of a given province, then it wouldn't make sense to try to convert it to a single polygon. However, I am not aware of Canadian provinces having pieces of other provinces in their interior.
I have seen posts with what looks like SQL server code; I don't know what application they are using. I do not have SQL server on my machine.


texts: How to handle interpolation of multiple (quality) tiers of irregular data?

texts:
I'm going to try to provide a high level overview of the issue - realizing that more details/info will be necessary - but I don't want to write a novel that seems too overwhelming for people to respond to :)

Background:

Let's say I have 3 sets (A, B, C) of irregular point data of varying quality; A is "better" than B and B is "better" than C. One or more data sets might have data for a given point.

Goal:

Create a grid/map that contains the "best" possible data.

Example/Possible Solution:

Use some sort of inverse distance weighted algorithm in combination w/somehow weighting the better datasets more than the lower quality data sets.

Please feel free to post follow-up questions as I know my original post is rather vague.  Thanks in advance.


texts: Routes through centroids and array functions

texts:
I have a table stopcentroids4 like in figure with the following fields: gid, routes, geom

and i want to draw a line for each route that runs through centroids. In my code there is a mistake. Actually the centroid on the lowest and right point of the picutre (where you can read the associated routes 47,57,76A, 92C) becomes the starting point of the routes 2,5,6,7.

The postresql code follows below:

CREATE OR REPLACE FUNCTION walk_centers(bigint, text)   RETURNS geometry AS $BODY$ WITH RECURSIVE next_stop(geom, idlist) AS (
                 (SELECT geom AS geom,
                    ARRAY[gid] AS idlist
                   FROM stopcentroids4
                  WHERE gid = $1)
        UNION ALL
                ( SELECT s.geom AS geom,
                    array_append(n.idlist, s.gid) AS idlist
                   FROM stopcentroids4 s,
                    next_stop n
                  WHERE string_to_array(s.routes,',') @> ARRAY[$2]
                  AND NOT n.idlist @> ARRAY[s.gid]
                  ORDER BY st_distance(n.geom, s.geom) ASC
                 LIMIT 1)
        )  SELECT ST_MakeLine(geom) AS geom FROM next_stop;
        $BODY$
        LANGUAGE sql VOLATILE   COST 100;
        ALTER FUNCTION walk_centers(bigint, text)   OWNER TO postgres;

This is my code:

CREATE TABLE lines AS
-- Distinct route identifiers! WITH routes AS (   SELECT DISTINCT unnest(string_to_array(routes,',')) AS route   FROM stopcentroids4 ORDER BY route ),
-- Joined back to stops! Every route has all its stops! stops AS (   SELECT s.gid, s.geom, r.route   FROM routes r   JOIN stopcentroids4 s  ON (strpos(s.routes, r.route) <> 0) ),
-- Collects stops by routes and calculate centroid! centers AS (   SELECT ST_Centroid(ST_Collect(geom)) AS geom, route   FROM stops   GROUP BY route ),
-- Calculate stop/center distance for each stop in each route. stops_distance AS (   SELECT s.*, ST_Distance(s.geom, c.geom) AS distance   FROM stops s JOIN centers c   ON (s.route = c.route)   ORDER BY route, distance DESC ),
-- Filter out just the furthest stop/center pairs. CREATE TABLE brighton_lines AS
-- Distinct route identifiers! WITH routes AS (   SELECT DISTINCT unnest(string_to_array(routes,',')) AS route   FROM stopcentroids4 ORDER BY route ),
-- Joined back to stops! Every route has all its stops! stops AS (   SELECT s.gid, s.geom, r.route   FROM routes r   JOIN stopcentroids4 s  ON (strpos(s.routes, r.route) <> 0) ),
-- Collects stops by routes and calculate centroid! centers AS (   SELECT ST_Centroid(ST_Collect(geom)) AS geom, route   FROM stops   GROUP BY route ),
-- Calculate stop/center distance for each stop in each route. stops_distance AS (   SELECT s.*, ST_Distance(s.geom, c.geom) AS distance   FROM stops s JOIN centers c   ON (s.route = c.route)   ORDER BY route, distance DESC ),
-- Filter out just the furthest stop/center pairs. first_stops AS (   SELECT DISTINCT ON (route) stops_distance.*   FROM stops_distance )
-- Pass the route/stop information into the linear route generation function! SELECT   row_number() over () AS gid, -- QGIS likes numeric primary keys   route,   walk_centers(gid, route) AS geom FROM first_stops;

-- Do some housekeeping too ALTER TABLE lines ADD PRIMARY KEY (gid);

Maybe the problem is in the function included in the join: "strpos(s.routes, r.route) <> 0"
and i should replace the regular function "string_to_array" but I'm getting syntax errors.


texts: Can I export QGIS basemaps as raster in ArcGIS 10.0?

texts:
I'm mapping a National Park in Malawi using Quantum GIS as with the plugin I can get the images from G.Earth (I need to see the vegetation).

Since my supervisor is using ArcGIS 10, do you know if it's possible to save the file with the maps I need (Malawi etc) so that she can open it in ArcGIS?
I'm wondering if there's a way to save the maps as raster or something...
Thank you.


texts: Interpreting results in legend, SAGA GIS or QGIS?

texts:
Can somebody explain to me what units of measurement do I get in the results of analysis, ground morphometry in SAGA GIS software or through a plugin in QGIS?
for example:

Terrain ruggedness index (TRI)? (m,%, feet, etc)
Topographic Position Index (TPI)? ...
TPI Based Landform Classification
Terrain surface Classification (Iwahashi and Pike)?

etc...

I have read the works related to the algorithms shown, but how do I interpret the results in the legend?

Topographic Position Index (TPI)

This is numeric units of Topographic Position Index (TPI)


texts: Create Cased Line Symbology for Feature Service in ArcGIS Online

texts: I'm trying to render a line feature service with a cased symbol.  For those of you unaware, a cased symbol is essentially like drawing an outline around the line.  It creates a very nice visualization.  My current workaround is to create a copy of the feature service and draw the copied layer slightly larger.  This creates the effect but now I have two layers in my contents pane and this is undesirable.  I'm surprised by ESRI's shortcomings in regards to the rendering options they provide in ArcGIS Online.  From what I've seen it looks like they provide very few options.


texts: Generating series of cross sections (profiles) along line using QGIS?

texts:
I have a set of points (xyz) along an axis.

Is there any way in QGIS to generate cross profiles along this axis each 10 meters with, for example, +-15 meters right and left from the axis?


texts:
I am mapping animals and have to display their locations in QGIS (2.12.3, Win7 64) with gender symbols (♂/♀). The outline of the symbols should be rendered in different colors corresponding to the dates of survey. The filling should be either colored or transparent respective to classes of behavior. So I generated the svg symbols with INKSCAPE and then replaced the whole „style“-line in the svg file text

style="fill:#000000;fill-opacity:0;stroke:#000000;stroke-width:10.00027084;… … ...“

with

fill="param(fill) #000" fill-opacity="param(fill-opacity) 0.0" stroke="param(outline) #000" stroke-width="param(outline-width) 3"

In QGIS I used attribute-driven rendering for the outline (attribute field “color“) and an expression to choose between either transparency or fill (case when “behav“ <> '+' then '0,0,0,0' else “color“ end)

The outline, both color and width, was displayed as I intended. In contrast, I had no control over the fill. Sometimes all symbols were filled, no matter what the attribute says, sometimes all fillings remained transparent and in the third case the symbols that should be transparent had a black (or white) fill. When I closed and reopened QGIS, the rendering of all fills changed sometimes. The outline, however, was displayed correct all the time.

Instead of '0,0,0,0' in the case-when expression above I experimented with other statements, e.g. '#000000 0.0' or 'default' to achieve correct rendering of color or transparency, but without success.
I also varied the default values in the “param(…)...“ statements and omitted the whole “fill-opacity...“ statement. Again without gaining control over the fill.

In addition, it was strange to see that sometimes the fill of all symbols changed from transparent to color or vice versa when I varied the outline width. I then introduced fill-rule=”evenodd”, but again without visible effect.

I have no idea what's wrong.


texts: Imported .nc rasters not editable in ArcGIS Pro

texts:
I imported .nc files from the OB.DAAC website and used the make NetCDF a new raster tool in ArcGIS Pro. These show up fine but they are not editable.
I right clicked on the geodatabase and added a new feature layer and it showed up under the gdb but I have no idea what to do with that or how to fix this problem.
So basically I have 2 rasters of the whole US of chlorophyll data. I only want to focus on the east coast of FL. One raster is from 2004 during a red tide bloom and the other is from 2010 where is was normal. All I want to do is create two rasters of just that coastal spot and compare the difference in chlorophyll amount. When I used extract by mask just now I may have put the whole rasters in the wrong spot but it created a layer that's greyscale and doesn't show up as anything when I turn all other layers off.


texts: Is it possible to set points without coordinates onto a map?

texts:
I have the following Situation:
csv-file with 1000 lines of measurements. These measurements were done indoor in a building without GPS. The measurement-system allows to set "markers" at any time.
Afterwards i would like to "georeference" the measurement with these markers in QGIS 2.10.1.
The csv looks like this (fake data, just to get an idea what i mean):

id time event value
0 12:00:00 marker 0
1 12:00:01 ? 2
2 12:00:02 ? 5
3 12:00:03 ? 2
4 12:00:04 marker 1
5 12:00:05 ? 9
6 12:00:06 ? 2
7 12:00:07 marker 7
...

Is it now possible with QGIS 2.10.1 to define the Position of the markers (red dots) through a mouseclick whereever I want on a room plan? Also would it be nice when the program distributes the measurements equally in between two markers as a line.


texts: Geocoding in ArcGIS desktop and use the shapefile in ArcGIS Online

texts:
I want add 60 addresses in point format to my map in ArcGIS online. So what I did was use ArcGIS Desktop (10.2) to geocode them, exported them as a shapefile, and imported them to the online map. First of all, is this the way I have to do it? I mean can I do it directly from ArcGIS online?

The second problem is I wanted to have 60 points. But I end up with fewer points, like 15-20. I found the problem is that somewhere during this procedure ArcGIS Desktop decided to round the X and Y values to 6 decimal places. I checked the geocoding results and changed the decimal places to 10. Then I checked the shapefile and I found in the shapefile's attribute table, X and Y just have 6 decimal places. I was not able to change the decimal place in attribute table. I tried to add new field and set the decimal places to 10 and then use filed calculator. It did not work!

I tried to geocode them using batchgeo.com, then I imported them as X and Y coordinates. I made sure Xs and Ys are set to have 10 decimal places. Agian the same problem when I exported it as a shapefile.

Any suggestions? The original file is addresses, including address, city, state and ZIP.


texts: Garmin gpx - wrong timestamp when importing it into Qgis

texts:
When I copy a track (day) from a Garmin Oregon 550 to the PC an import it into Qgis a the timestamp is 2 hours earlier than it should be.

When I import it via Garmin Basecamp, export it and open it in Qgis the time stamp is still wrong. When I open the track in Basecamp the timestampes are correct.

When I import the same file in Geosetter (http://www.geosetter.de/) and synchronize images with the track it works fine (right positions).

The time on the device is the right time zone.

Question:
Is there a simple way to import the track directly to Qgis with the right timestamps?


texts: Is there any support for Rasterlite in the Android QGIS?

texts:
Is there any support for Rasterlite databases in the Android port of QGIS?

I have tried adding a raster layer as I would on desktop QGIS but there are fewer GDAL file formats listed.

Is this because rasterlite is not available on Android or is there another way to bring it in?

Thanks

Tom


texts:
I am looking for a way to transform WKT-strings from EPSG:31467 to EPSG:25832 using a ntv2-grid (and without using a gui).

Example:

POLYGON((3413732.42 5318051.77,3413747.65 5318049.04,3413750.05 5318062.33,3413734.96 5318065.13,3413732.42 5318051.77))

For shapefiles I use ogr2ogr but have no clue how to transform WKT strings.


texts: Selecting by attribute using Python and a list

texts:
I have been using this script to select a feature using Python:

layer = iface.activeLayer()
layer.selectByExpression('\"Declividad\"= value', QgsVectorLayer.SetSelection)
selection = layer.selectedFeatures()

But now, I need to select the values for Declividad from a list.

For exemple: list = [10,11,12], so I want to select the values 10,11 and 12 for Declividad.

How could I do that?


texts: Opening MOD17A3 in ENVI so it is properly georeferenced?

texts:
I am using ENVI 5.1 on Windows 8. I am trying to open a MOD17A3 (MODIS Net Primary Productivity standard product - L4) .hdf file which I downloaded from the USGS website. I am still a beginner at ENVI, IDL, Remote Sensing, and GIS.

My problem is when I try to use the ENVI GUI and open the file.

When I open the .hdf file as a EOS>MODIS file, it appears properly georeferenced when overlaid with other vector layer that I am using, BUT it does not show the data for each pixel. What I mean is, the shapes of the landmasses are present, but they are all in black, not showing any NPP data which should have made the pixels appear in varying shades of grey.

When I open the .hdf file as Generic Format>HDF4, ENVI asks me to select which datasets in the .hdf file I should open, and I select NPP. The values of each pixel are displayed, such that they appear in varying shades of grey, BUT they don't seem to be properly georeferenced since they don't align with other georeferenced layers. The map has shrunk, but is in the same general location as the other layers.

It might be possible to use the IDL controls to solve this problem but I wouldn't have any idea where to start.


texts: I have been trying to use the pygeostat, it seems to me that they need userid/password when trying to download package from Geostatistics Knowledge Base. But I could not find the contact information for this.


texts: Showing global temperature data in Leaflet layers by GRIB2 files?

texts:
Wind data already shown in map by danwild/leaflet-velocity

Now I want add temperature to map by tiles.

How can I convert GRIB2 files to tiles?


texts: How to add a time component in a heatmap?

texts:
I have a point file extracted from an AIS network. 1 year is about 900 million points. Each point has a timestamp.

I want to make a heatmap for a longer period of time, but I want the heatmap to take the timestamp into account. I only want the heatmap to show me hotspots if the points are within a certain distance within a limited timeframe.

Edit 20190815:

I'll try to explain my needs further with a hypothetical example:
Imagine a school with a centre hallway and 2 classrooms on either side of the hallway (4 classrooms in total).

Situation 1:

At 1500 the bell rings, all the children jump out of their seat, run to the door into the central hallway, and run to the main door that leads to the playground.

Situation 2:

At 1500 the bell rings, and one by one, the children get up from their seat, walk to the door into central hallway, 1st classroom 1, then 2, then 3, then 4. Forming a nice neat line to the main door.

Situation 3:

At 1500 the bell rings, the two classrooms on the left side of the centre hallway behave as in Situation 1, and on the right side of the central hallway behave as in Situation 2.

 In a normal heatmap, that doesn't take a timefactor into account, you
 will see similar results on all three situations.

What I am trying to accomplish is show congestion.
I think that:
If I am able to get the timestamp of every datapoint into the mix the heatmap will reveal in:

Situation 1:
Hotspots in front of all the doors that lead to the central hallway and a hotspot in front of the main door to the playground.
Situation 2:
No hotspots.
Situation 3:
Hotspots in front of the doors on the left side of the central hallway and a smaller hotpot in front of the main door (because the chaotic behavior of the left classrooms, influences the behavior of the right classrooms)

Is there a workflow, methodology or algorithm available that can build a heatmap with a timefactor?


texts: Conditional label expression gives error in ArcGIS

texts:
I have a polygon dataset and I would like to label the area of all polygons which have "kamer" in the column "type".

I tried some expressions but I get different Errors.
Two examples of what I tried is the following:

Function FindLabel ( [Area], [type] )
  if ([type] = 'kamer') then
   FindLabel = UCase([Area])
  end if
End Function

This gives me the following error:

  The expression contains an error.
    Modify the expression and try again.
    Error 1 on line 16.
    Syntaxisfout.

Another expression I tried is:

def FindLabel ([Soort_kamer], [Oppervlakte]):
  if int([Soort_kamer]) = 'kamer':
    return [Oppervlakte]


texts: Merging Raster Datasets of Different Resolutions using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop.

I have two DEMs in .tif format and I am looking to merge the two. The main problem with merging these are that they are different resolutions (one is 10m and one is ~20m). Also, another problem with merging the two is that they overlap and I am not sure of the most effective way of dealing with that overlap. Obviously I want to keep the values of the better resolution.

So what are the best methods here?

Is there a python script available that works best with merging rasters? I would like to automate this with a script if possible.
Which process of resampling is best for DEM data? Bipolar, Cubic... ?
What is the best way to merge overlapping DEMs?


texts: Doing OpenLayers.Request in OpenLayers 3

texts:
I have been searching for solutions to call GeoServer REST api from OpenLayers 3 for a long time.

I found that for OpenLayers 2, there is a OpenLayers class OpenLayers.Request for doing this:
http://docs.openlayers.org/library/request.html

But for OpenLayers 3, I can't find a similar class.
I tried a lot of online code for PHP proxying (andjsonp) trying to do the same without OpenLayers.Request, to no avail.

What is the correct/most reliable way to call GeoServer Rest from OpenLayers3 (or even general JavaScript with jQuery/AJAX)?