symbology - QGIS just installed on MAC OSX - vector layer styles do not work properly

I am a newbie with QGIS - I have just installed QGIS 2.10.1-Pisa. I am following a youtube tutorial by Steven Bernard and I just found that the styles do not work properly (I run Mac OS X Lion 7.5.1. on a 1.7 GHz Intel Core i5 MacBook Air).

When I want to change the sybology of my map, this does not work properly. If you have a look at step 3.2.3. Changing Symbol Structure in the manual on 3.2.1. basic Follow Along: Changing Colors , QGIS does not provide me with any choices. It tells me that 'This layer doesn't have any editable properties'. I am attaching a screenshot just to be clear.

Can anybody tell me what is wrong?

Answer

I had a similar error with pseudocolor properties on macboook air mid 2011 OS x 10.10.4

This solved my problem:

openstreetmap - Creating columns in PostGIS using osm2pgsql hstore tags

While I feel this may come down to a database question, it feels GIS-related enough that I feel comfortable posting it here. I'll be happy to take it elsewhere if it belongs there, though!

I've started the arduous journey of mangling OSM data to work in an RDBMS format that is useful to me. One of the major problems I am running into is the hstore data type. While it is quite useful for being able to hold an extremely large dictionary of values, they are hard to "get at", in my experience. Specifically, the osm2pgsql tool creates addr:flats, addr:housenumber, and addr:interpolation fields, but leaves out some that are quite obvious to me, such as city, postcode, and street. I would like to parse this data to populate new columns with it. What methods would you recommend for accessing such data?

Thanks

Answer

You can transfer the desired hstore key/value data to new columns:
1-Create the desired columns (ex. addres, city, key1, key2, keyn)
2-Run:

UPDATE table 
SET 
address=hstorecolumn->'address',
city=hstorecolumn->'city',
key1=hstorecolumn->'key1',
key2=hstorecolumn->'key2',

keyn=hstorecolumn->'keyn';

Or you can simply get used to hstore, it works pretty well and I have so many good results with it that I can say IMO "hstore changed my life". lol

EDIT:
An example on how to use hstore:
SELECT
streetdatatable.hstoredatacolumn -> 'address',
streetdatatable.hstoredatacolumn -> 'city'
FROM someschema. streetdatatable

It will show you all the address and cities values.

hide attribute with null value using leaflet

http://jsbin.com/hinekaj/edit?html,output

this is a link to a map am working on. It will be great if someone could help me to hide the attributes that have "null" as their value. I am aware that I need to use the IF command but have been struggling with it.

        onEachFeature: function (feature, layer) {
          layer.bindPopup('Address: '+ feature.properties.address 
                          + '
Sample 1: '+ feature.properties.sample_1 + ' Pb'
                          + "
Sample 1: "+ feature.properties.sample_1 + ' Pb'
                        + " , Sample 2: "+ feature.properties.sample_2 + ' Pb'
                         + "
Sample 3: "+ feature.properties.sample_3 + ' Pb'

                         + " , Sample 4: "+ feature.properties.sample_4 + ' Pb'
                         + '
Sample 5: '+ feature.properties.sample_5 + ' Pb'
                         + ' , Sample 6: '+ feature.properties.sample_6 + ' Pb'
                         + '
Sample 7: '+ feature.properties.sample_7 + ' Pb'
                         + ' , Sample 8: '+ feature.properties.sample_8 + ' Pb'
                         + '
Sample 9: '+ feature.properties.sample_9 + ' Pb'
                         + ' , Sample 10: '+ feature.properties.sample_10 + ' Pb'
                         + '
Sample 11: '+ feature.properties.sample_11 + ' Pb'
                         + ' , Sample 12: '+ feature.properties.sample_12 + ' Pb'
                         + '
Sample 13: '+ feature.properties.sample_13 + ' Pb'

                         + ' , Sample 14: '+ feature.properties.sample_14 + ' Pb'

         ); 
        },

Snapping points to lines in ArcGIS Desktop and automate using VBA?

I use ArcGIS Desktop and I want create a network with a point feature and a line feature. I need my point features to be on the line features so I should do Snap point to line.

However, I do not know how do it in ArcGIS, and then automate it with code for this in VBA.

Converting PostGIS table to Shapefile in Python?

I want to convert a PostGIS table into a shapefile (without using pgsql2shp).

In order to create a geometry in the shapefile I have to give the Xmin,Ymin and Xmax,Ymax, and the geometry which I have in my PostGIS table is an irregularly shaped one (I can get the exterior using bounding box but that will include some extra area more than my area of ineterest). Is there any method by which I can get the task done?

I want to do the thing programmatically and using Python.

ValueError: array larger than output file, or offset off edge Python GDAL

I have been working on this iterator that would count instances where data conformed to a set of conditions over a number of rasters. I had people help me with other aspects of this script (Counting events from a GRIB raster), but now I have hit an error where the return is "ValueError: array larger than output file, or offset off edge", Which I can't seem to be able to work arround. The code is provided bellow. It requires the user to input the path to folder and provide a name for the output GeoTIFF raster. Files read are GRIB.

import gdal
import numpy
import osr
import os

del_mapa = str(input('Vpiši pot do datotek v obliki "C:\\Mapa1\\Mapa2\\...\\Zadnja mapa": '))
os.chdir(del_mapa)


pix_vel = 1000
#Defines pixel size
ime_rez = str(input("Poimenuj datoteko z rezultati: ") + ".tif")
#User can provide a name for the output picture
format = "GTiff"
driver = gdal.GetDriverByName(format)

vir_mere = gdal.Open('inca_20140101-0100.grb')
trans = vir_mere.GetGeoTransform()

geotrans0 = trans[0]
geotrans1 = trans[1]
geotrans3 = trans[3]
geotrans5 = trans[5]
#Extracts transformation data from one of the GRIBs
dst_ds = driver.Create(ime_rez, 401, 301, 1, gdal.GDT_CInt32)
dst_ds.SetGeoTransform([geotrans0, geotrans1, 0, geotrans3, 0, geotrans5])
#and pastes it directly into the newly created GeoTIFF
projInfo = vir_mere.GetProjection()
srs = osr.SpatialReference()

srs.ImportFromWkt(vir_mere.GetProjectionRef())

dst_ds.SetProjection(srs.ExportToWkt())

raster = numpy.zeros((401, 301), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)

def krog(baza,novi):    
    x = 0
    while x < 25:

        baza = gdal.Open(ime_rez)
        if x < 10:
            novi = gdal.Open('inca_20140101-%s00.grb') % ('0' + str(x))
            #This should open the next file in line
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0

        elif x > 10:

            # novi = gdal.Open('inca_20140101-%s00.grb') % (str(x))
            if novi > 15 and novi < 25:
                baza = baza + 1
            else:
                baza = baza + 0
        x += 1 
dst_ds.GetRasterBand(1).WriteArray(baza)        
dst_ds = None

arcmap - How to create parameter windows for custom geoprocessing tools?

I would like to create a geoprocessing tool with a custom GUI for its input data. To put it in other words: imagine that your client works mainly with arcmap documents and he already has an arcpy script which deals with the geoprocessing. Unfortunately he doesn't like building blocks which take parameters for the script (drop-down list, onelineinput etc.) and would like to have a custom input method like let's say a scrollbar (f.e. for selecting time interval instead of typing it manually or selecting from a list).

Is it achievable in ArcMap to have custom input GUI elements? I'm pretty sure it's not, but my colleague swears on everything that he's seen something like that.

Before I completely cross out that idea I wanted to have a confirmation from people who have spent more time working with ArcGIS. To my understanding the only solution is to create a standalone application using Runtime SDK or ArcObjects/ArcEngine.

Answer

I think that if you are creating a Python script tool, your options are limited to what user interfaces ESRI provides for each type of parameter.

I think it is possible to create custom geoprocessing tool controls, but not without ArcObjects and COM.

If you don't need to use the geoprocessing framework, you could use whatever UI library you want (Windows Forms, WPF, Java Swing, Qt, Wx, Tk, etc.).

If you are comfortable using Python and want to try creating a Python add-in with a custom GUI, this example should get you going:

Note: I dunno what happened there with the button being cut off :)

Calculating sequential numbers into sorted table using ArcGIS Desktop?

Is there a way to calculate a sorted field with sequential numbers? I have seen Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator? that outlines how to calculate sequential numbers, but this is always calculated on FID order, not on sorted order.

#Pre-logic Script Code:
rec=0
def autoIncrement(): 
    global rec 

    pStart = 1  
    pInterval = 1 
    if (rec == 0):  
        rec = pStart  
    else:  
        rec += pInterval  
    return rec

#Expression:
autoIncrement()

An example of what I'm trying to do. I've used an advanced sort to sort by year, month, day, and now want to have sequential numbers in the Seq field. You'll see that my OBJECTID field is not in order, so the above code won't work.

Can this be done either in the Field Calculator or using an Update Cursor in arcpy?

Answer

"Solution" with 2 sorted fields (ascending):

mxd = arcpy.mapping.MapDocument("CURRENT")
lr=arcpy.mapping.ListLayers(mxd)[0]
tbl=arcpy.da.TableToNumPyArray(lr,("oid","A","B"))
bs=sorted(tbl, key=lambda x: (x[1], x[2]))

def sortSeq(fid,a,b):
 for i,ent in enumerate(bs):
   if ent[0]==fid: return i

sortSeq( !OID!, !A!, !B! )

UPDATED VERSION:

mxd = arcpy.mapping.MapDocument("CURRENT")

lr=arcpy.mapping.ListLayers(mxd)[0]
tbl=arcpy.da.TableToNumPyArray(lr,("oid","A","B"))
bs=sorted(tbl, key=lambda x: (x[1], x[2]))
aDict={}
for i,row in enumerate(bs):
 aDict[row[0]]=i
def sortSeq(fid):
 return aDict[fid]

sortSeq( !OID!)

takes 1.5 seconds to complete task on 10000 records. Original takes slightly more than 2 minutes

GetFeatureInfo Request with GeoServer and OpenLayers.loadURL Not Working

I'm making a simple GetFeatureInfo request against a GeoServer WMS, but it returns empty responses. The request utility I'm using is OpenLayers.loadURL.

The funny thing is, I can paste the generated request URL from FireBug into a new browser tab and get a valid response. So, it seems to be something with the OpenLayers.loadURL. See the examples below.

I can get a valid response with this map event where the GetFeatureInfo request is in a new window:

  map.events.register('click', map, function (e) {
    var url = "http://localhost:8080/geoserver/MYLAYER/wms" 
    var params = "" 

      + "REQUEST=GetFeatureInfo"
      + "&SERVICE=WMS"
      + "&EXCEPTIONS=application/vnd.ogc.se_xml"
      + "&VERSION=1.1.1"
      + "&BBOX=" + map.getExtent().toBBOX()
      + "&X=" + Math.round(e.xy.x)
      + "&Y=" + Math.round(e.xy.y)
      + "&INFO_FORMAT=text/html"
      + "&QUERY_LAYERS=NERD:NERD Risk Rates"
      + "&LAYERS=NERD:NERD Risk Rates"

      + "&FEATURE_COUNT=50"
      + "&format=image/png"
      + "&SRS=EPSG:900913"
      + "&STYLES="
      + "&WIDTH=" + map.size.w
      + "&HEIGHT=" + map.size.h;

    console.log(url+"?"+params);
    window.open(url+"?"+params,
      "getfeatureinfo",

      "location=0,status=0,scrollbars=1,width=600,height=150"
    );

   });

However, this map event using OpenLayers.loadURL always returns an empty response and no error. But if I cut and paste the GET request url from FireBug into a new browser window, I get the proper response:

  map.events.register('click', map, function (e) {
    var url = "http://localhost:8080/geoserver/MYLAYER/wms" 

    var paramHash =  {

      REQUEST : "GetFeatureInfo" ,
      EXCEPTIONS : "application/vnd.ogc.se_xml" ,
      SERVICE : "WMS" , 
      VERSION : "1.1.1" ,
      BBOX : map.getExtent().toBBOX() ,
      X : Math.round(e.xy.x) ,
      Y : Math.round(e.xy.y) ,
      INFO_FORMAT : "text/html" , 
      QUERY_LAYERS : "NERD:NERD Risk Rates" , 
      LAYERS : "NERD:NERD Risk Rates" ,

      FEATURE_COUNT : 50 , 
      format : "image/png" ,
      srs : "EPSG:900913" ,
      styles : "" ,
      WIDTH : map.size.w , 
      HEIGHT : map.size.h
    };

    OpenLayers.loadURL(url, paramHash, this, setHTML, setHTML);
        OpenLayers.Event.stop(e);


   });

Someone have any ideas?

Answer

Sounds like a same origin problem to me. Are you serving your page from a different host/port combination than the map server? See http://trac.osgeo.org/openlayers/wiki/FrequentlyAskedQuestions#ProxyHost for details of how to fix this using a proxy or move your html page to data/web in your geoserver installation.

python - QGIS docutils install (path,pythonpath,pythonhome)

New to QGIS and Python. Have QGIS installed and working great on Windows, want to install and use Docutils -- Python Documentation Utilities (presumably all simple and small .py programs).

I see that QGIS has a Python Console built in. Do I just use that? What is the command line to run install.py?

QGIS obviously installed Python. However it looks like the customary path variables were not set during the QIS install, such that Python can be run from a windows command (cmd.exe) window.

That seems to be the task here, so I proceed to try to enable Python to run independent of QGIS.

First I search for python.exe. I find it at: C:\Program Files\QGIS Lyon\bin

So I set the windows path variable to include that.

After searching (and a lot of confusion as to where the Python distribution actually is, is it standard, or is it custom, questions as to how libraries, etc are laid out, and noting its all installed "underneath" C:\Program Files\QGIS Lyon) so it appears custom(?) Will I break that by using it for "generalized" Python programs?

I try the following environmental values:

PYTHONPATH C:\Program Files\QGIS Lyon\apps\Python27\Lib

PYTHONHOME C:\Program Files\QGIS Lyon\bin

Have no idea if those are even close to correct. I then run the install for docutils

C:\Users\rick\Desktop\docutils-0.12>python install.py

Program starts, but then errors per below.

Seems like a "documentation utility" should be a snap to install. Plus there probably a straightforward way to run Python programs independent of QGIS? I also note an additional install of Python with Inkscape (a vector drawing program used in the QGIS training manual.) So already two installs of Python on the system already. (Correction, three installs. FWTools has yet another Python install.)

If a user wants to run Python programs (such as docutils and other small learning programs) should they download a fresh Python "stand alone" install? (looks like on this machine that would be the 4th install?)

---

errors

C:\Users\rick\Desktop\docutils-0.12>python install.py

This is a quick & dirty installation shortcut. It is equivalent to the command::

python setup.py install

However, the shortcut lacks error checking and command-line option processing. If you need any kind of customization or help, please use one of::

python setup.py install --help
python setup.py --help

Traceback (most recent call last): File "install.py", line 27, in dist.run_commands() File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 95 3, in run_commands self.run_command(cmd) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 97 2, in run_command cmd_obj.run() File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install.p y", line 575, in run self.run_command(cmd_name) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\cmd.py", line 326 , in run_command self.distribution.run_command(command) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dist.py", line 97 2, in run_command cmd_obj.run() File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install_l ib.py", line 97, in run outfiles = self.install() File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\command\install_l ib.py", line 115, in install outfiles = self.copy_tree(self.build_dir, self.install_dir) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\cmd.py", line 377 , in copy_tree dry_run=self.dry_run) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dir_util.py", lin e 139, in copy_tree mkpath(dst, verbose=verbose) File "C:\Program Files\QGIS Lyon\apps\Python27\Lib\distutils\dir_util.py", lin e 76, in mkpath "could not create '%s': %s" % (head, exc.args[-1])) distutils.errors.DistutilsFileError: could not create 'C:\Program Files\QGIS Lyo n\bin\Lib': Access is denied

Answer

Don't make life more complicated, it's much easier:

• On Unix systems (Linux and Mac OS X) QGIS use the standard Python, installed by default. This is not the case in Windows where Python is not installed by default and there is no PYTHONPATH variable.


• Therefore QGIS installs its own version (with libraries in C:\Program Files\QGIS Lyon\apps\Python27 in the Standalone version and a different location in the OSGeo4W installer)

So:

1) You can use the Python console as a classic Python shell using exclusively the QGIS Python with PyQGIS and the other modules installed by default, as Numpy (in C:\Program Files\QGIS Lyon\apps\Python27\Lib\site-packages) = the PYTHONPATH of the QGIS Python. You can control it by typing in the console

import sys
print sys.path

2) if you want to install an other Python module as Docutils, the problem is that Setuptools (setup.py, easy_install) or pip are not installed by default, even with the OSGeo4W installer, -> look at How to install 3rd party python libraries for QGIS on Windows?, QGIS Standalone and the Python Modules, Installing Python setuptools into OSGeo4W Python or OSGeo4W External Python Packages

At this stage, you have a full Windows Python version: you can use setup.py, easy_install or pip to install new modules, you don't need to fix the PYTHONPATH and you can use PyQGIS and the other modules installed in the QGIS site_packages folder

1. via commands in the Python console


2. via Python scripts in the Processing Toolbox or with the ScritRunner plugin of Gary Sherman


3. via the development of custom Plugins

If you want to use other simple modules installed in other Python versions, simply use:

sys.path.append(path_of_the_module)

3) But if you want to use PyQGIS from outside (in other version of Python), it is another problem: search for on the GIS SE or the Web (as in Configure PyScripter to use with QGIS (and still use arcpy) on Windows, Standalone applications using QGIS and environment variables, and many others)

Why is this GRASS script (from GRASS WiKi) causing an error?

I'm trying to run a script from the GRASS WiKi to create a simple map using ps.map. I'm confused because it's generating an error that the output was not set, when it appears to me that it is set at the beginning of the script. The script and error follow. The first code box contains the contents of the file spearfish.map, which is referenced as the input value in the second code box.

# Simple ps.map example using the Spearfish dataset
g.region rast=elevation.10m

ps.map output=spearfish.ps << EOF
paper us-letter
  end
scale 1:120000
raster elevation.10m
colortable
  where 2 6.75

  end
text 50% -7% Spearfish County, ND (Mount Rushmore)
  fontsize 16
  end
mapinfo
  where 4.5 7.25
  end
vlines roads
  where label ~ 'highway' OR label = 'interstate'
  color grey

  end
vpoints archsites
  symbol basic/triangle
  end
end
EOF

The command for running the script and the error it generates are below:

> ps.map input=spearfish.map                                                     
ERROR: Required parameter 
 not set:

    (PostScript output file)

UPDATE: Based on ShaunLangley's comment, I tried ps.map output=out.map. This seems to have started ps.map in interactive mode and I was able to enter the rest of the script one line at a time. That did create a map as a postscript file, so I'm headed in the right direction. I'll close this thread as solved and open a new one more targeted on the question of how to run the commands as a script.

Answer

the ps.map function requires you specify output. The example you have:

> ps.map input=spearfish.map

doesn't specify the output variable. I'm not sure what the distinction between the code boxes is. The former should work because the input raster is defined with "raster" and output defined at the top. It looks like you're trying to replicate an example; is this correct?

Perhaps you could try:

ps.map input=elevation.10m output=spearfish.map

You can specify the additional options as well just as the code example did.

Extracting centerline of a Complex-Polygon in PostGIS/Python

Any idea how to extract the centerline of the following polygon (Fig1) using PostGIS functions or Python modules? I was trying Voronoi Function but thats not helpful for this kind of shape. I want something like Fig2 blackline.

shapefile - merging 2 polygon adjacent to each other using R

I have several shp files which are the components (or tiles) of a map, I have joint all these tiles together to have a large shp file. Then I realized some of the polygon is now split into 2 halves, as the upper part belongs to a tile, while the lower part belongs to another tile.

This output disturbs me a lot as I need to fill some polygon, but not part of the "split polygon", so I need to find a way to merge these split polygon back to one (preferably mechanical way, as I have 10000+ polygons, not knowing which are split, which are not).

I need to know if the two or more polygons originate from the larger polygon. And I think I can somehow use union and spRbind to create an algorithm for this task. But I doubt if I am really the first one who encounters this problem, so I wonder if there is any existing tool in R that can solve this.

Thanks.

Answer

I can only provide a partial answer to your question.

There is a function in the maptools R package called unionSpatialPolygons, see here for a manual.

The function expects an ID field, and all polygons with the same ID value will be merged. If you have an ID that identifies your polygons to be dissolved, use it. Otherwise you could try something like

unifiedPolygons <- unionSpatialPolygons(myPolygons, rep(1, length(myPolygons))

but that will probably also dissolve all other polygons that touch within you area.

Spatially joining polygon fields onto point layer in QGIS?

I have gridded depth points in CSV format that I'm trying to add other map data to, in additional columns. I've tried joining fields, hoping that I could link the Shapefile data to the csv, but when I try to display the csv by graduated symbols based on the new fields, it goes blank. I used Nathan's tips from here, but imported the csv as csv since I want to add map info to that, rather than adding csv info to a map.

See this picture below for the sea bottom temperature & salinity voronois by month, underneath the csv points (a sample subset), with the joins already made. I added columns on the csv for each of the fields i'm intending to append.

I also tried spatial joins as per this however it says the layers aren't the same CRS (I checked, they are), then gets to 15% and hangs.

I also also tried spatial join through mmqgis but it hit a python error. I can try this again and look to debug the problem if this is the right way to go.

In short: appending polygon data as fields on a points-grid csv: what's the best way in QGIS 2.0?

Answer

The best tool for this job in my experience is Add polygon attributes to points in the Processing toolbox. If it does not work directly with the CSV, just save the points to a Shapefile before you run the spatial join.

landsat 7 - remove one image from imageCollection in google earth engine

I work with google earth engine on an imageCollection from landsat7. Some of those images are 'corrupted'. So I would like to remove one image from an image collection.

Actually I have :

var collectionLS7NDVI0 = ee.ImageCollection('LANDSAT/LE7_L1T_32DAY_NDVI')
    .filter(ee.Filter.calendarRange(6,6,'month'))
    .map(addTime);


var listOfImages = collectionLS7NDVI0.toList(collectionLS7NDVI0.size());
print(listOfImages)
0: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/19990610 (2 bands)
1: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20000609 (2 bands)
2: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20010610 (2 bands)
3: 

Image LANDSAT/LE7_L1T_32DAY_NDVI/20020610 (2 bands)
4: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20030610 (2 bands)
5: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20040609 (2 bands)
6: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20050610 (2 bands)
7: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20060610 (2 bands)
8: 

Image LANDSAT/LE7_L1T_32DAY_NDVI/20070610 (2 bands)
9: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20080609 (2 bands)
10: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20090610 (2 bands)
11: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20100610 (2 bands)
12: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20110610 (2 bands)
13: 

Image LANDSAT/LE7_L1T_32DAY_NDVI/20120609 (2 bands)
14: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20130610 (2 bands)
15: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20140610 (2 bands)
16: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20150610 (2 bands)
17: 
Image LANDSAT/LE7_L1T_32DAY_NDVI/20160609 (2 bands)

And I want to remove images 6 (2005).

Any ideas are welcome

Etienne

Answer

You can remove images of a collection by constructing a filter that selects the images you want to remove, then negate the filter using ee.Filter.not().

var MAX_LIST_SIZE = 100;

var collectionLS7NDVI0 = ee.ImageCollection('LANDSAT/LE7_L1T_32DAY_NDVI')
    .filter(ee.Filter.calendarRange(6,6,'month'));
print('listOfImages', collectionLS7NDVI0.toList(MAX_LIST_SIZE));


// Filter out image(s) from 2005.
var filtered = collectionLS7NDVI0.filter(
  ee.Filter.date('2005', '2006').not()
);
print('filtered', filtered.toList(MAX_LIST_SIZE));

QGIS Difference not working?

I'm trying to simply get the Difference between two sets of lines, but the algorithm returns a copy of either one, depending of the combination of INPUT/DIFFERENCE layer I try. Never the difference.

Both of my shapefiles are in the same projection, and have been cleaned (v.clean).

I tried to dissolve them, but no luck.

Any ideas?

(Using QGIS 2.18.1 over Win7)

---

Here is a picture of the situation. I have 2 sets of roads (red and green). The green network overlaps everything from the red one, but has some other roads in it. I'm trying to get a shape of these extra segments (hence, the Difference). But the result gives me either an exact copy of the red or green network.

wfs - How to edit feature attributes with openlayers?

I want to edit object's attributes on WFS layer. I see this example http://gis.ibbeck.de/ginfo/apps/OLExamples/OL26/examples/styles_unique_with_group_wfs.html but popup is horrible. I want to use ExtJs form to editing.
So how to edit attributes in this case?

Answer

At first need to create store

                          st = new GeoExt.data.FeatureStore({
                                    layer: myVecLayer,
                                    fields: [
                                        {name: 'id', type: 'int'}
                                    ],
                                     proxy: new GeoExt.data.ProtocolProxy({
                                        protocol: new OpenLayers.Protocol.WFS({
                                            version: "1.0.0",
                                            srsName:"EPSG:900913",

                                            url: "http://localhost:8080/geoserver/wfs",
                                            featureType: "filedata_temp",
                                            featureNS: "http://www.opengeospatial.net/cite",
                                            filter:Ffilter
                                        })
                                    }),
                                    autoLoad: true
                                });

Then EditorGrid

    var cm = new Ext.grid.ColumnModel({
                columns: [
                    {header: "id", 
                     dataIndex:"id",
                     editor: new Ext.form.TextField({allowBlank: false})
                    }]
    })


var ed_tab = new Ext.grid.EditorGridPanel({

                store: st,
                region: 'center',
                    cm: cm,
                height: 200,
                title:'Аттрибуты',
                    autoScroll: true,
                bbar:[{
                     xtype: 'tbbutton',
                     cls: 'x-btn-icon',
                     tooltip: "Сделать что то",

                     id: "something",
                     icon: url_servlet+'externals/gxp/src/theme/img/silk/delete.png',
                     handler: function(){
                                        ed_tab.store.proxy.protocol.commit(
                                               myVecLayer.features
                                          );
                                        ed_tab.store.reload();
                                    }
                    }]
                })



var x = new Ext.Window({
        title:'Аттрибуты',
        items:[ed_tab]  
        })
 x.show();

That all about attributes. To modifity geometry htere: Modify WFS with OpenLayers

python - Is there a way to locally add links to an OSM data set

Here's what I'm trying to do:

• I'm modeling some transportation logistics data


• I want to test the effects of a "new highway" on existing travel demand


• The new road will be created from an arbitrary point (usually not on the existing roadway) to any point on the road network (road line layer)


• The new road should acquire the same fields / attributes as the road it connected too

Is this possible with OSM data?

How to export and import gps points from Garmin 62S and QGIS

I'm pretty new to using QGIS and I'm curious about how to import and export GPS points created in QGIS. I have QGIS 2.0.1 Dufour on Windows XP and a Garmin 62s. Any suggested walk throughs would be great.

javascript - Access GeoJSON from GitHub in Leaflet

I am new to JSON, GitHub, and Leaflet.

I am trying to create a simple application which will render my points which are in GeoJSON up on GitHub.

From everything I have read this should work. Here is what I have so far but I do not see my points displayed.

What am I missing?

attribute table - Color in QGIS based on record of events in CSV file?

I have a shapefile with the following attributes:

ID    Region
1     Region1
2     Region2
3     Region3

Also I have the following record of events in a CSV file

EventID     Year     Month     RegionID     Quantity
001         2015     12        1            6
002         2015     12        2            7
003         2015     12        3            3
004         2015     11        1            4
005         2015     11        3            3
006         2015     10        2            6
007         2014     12        3            3
008         2014     11        2            2
009         2014     10        1            7

...         ...      ...       ...          ...

I would like to color my map acording to a filtered record. For instance, say I want to display all the events of November 2015. Then I would have:

ID    Region    Quantity
1     Region1   4
2     Region2   0
3     Region3   3

Or say I would like to display all the events of 2015. Then I would have the sum of the events:

ID    Region    Quantity

1     Region1   10
2     Region2   13
3     Region3   6

Is there a way to solve this problem dynamically, without rewriting the attribute table of the shp every time?

Consider that the CSV file is updated frequently.

open source gis - Viewing Personal Geodatabase using software other than ArcGIS?

Is it possible to open a Personal Geodatabase file in software other than ArcGIS?

i.e are there any open source software available just for viewing the data?

wms - How to get getGetFeatureInfoUrl from multiple layers in OpenLayers 3

I have multiple layers in a map with layer group now onclick I want to get the getGetFeatureInfoUrl for the different layers on which I am clicking but the problem is that I am getting only the Info for the last layer appended not the others.

Here is my code:

 var wms_layer;
    var map = new ol.Map({
        target: 'map',
        layers: [
        new ol.layer.Group({
            'title': 'Base maps',
            layers: [
            new ol.layer.Tile({
                title: 'Water color',
                type: 'base',

                visible: false,
                source: new ol.source.Stamen({
                    layer: 'watercolor'
                })
            }),
            new ol.layer.Tile({
                title: 'OSM',
                type: 'base',
                visible: true,
                source: new ol.source.OSM({

                    attributions: [
                    'All maps © OpenCycleMap',
                    ol.source.OSM.ATTRIBUTION
                    ],
                    url: 'https://{a-c}.tile.thunderforest.com/cycle/{z}/{x}/{y}.png'
                })
            })
            ]
        }),
        new ol.layer.Group({

            title: 'Hisar Data Layers',
            layers: [
            new ol.layer.Tile({
                title: 'Habitations',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',
                    params: {
                        'LAYERS': 'HissarFullData:habitations',
                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',

                    },
                    serverType: 'geoserver'
                })
            }),
           wms_layer = new ol.layer.Tile({
                title: 'Hisar PWD',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',
                    params: {
                        'LAYERS': 'HissarFullData:pwd',

                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',
                    },
                    serverType: 'geoserver'
                })
            }),
            new ol.layer.Tile({
                title: 'Hisar Block',
                source: new ol.source.TileWMS({
                    url: 'http://192.168.2.2:8080/geoserver/HissarFullData/wms?service',

                    params: {
                        'LAYERS': 'HissarFullData:block_boundary',
                        'VERSION': '1.1.0',
                        'FORMAT': 'image/png',
                    },
                    serverType: 'geoserver'
                })
            })
            ]
        })

        ],
        view: new ol.View({
            center: ol.proj.transform([75.7217, 29.1492], 'EPSG:4326', 'EPSG:3857'),
            zoom: 10,
            minZoom: 5,
            maxZoom: 17
        })
    });
    var layerSwitcher = new ol.control.LayerSwitcher({
        tipLabel: 'Layers' // Optional label for button

    });
    map.addControl(layerSwitcher);
    var viewProjection = map.getView().getProjection();
    var viewResolution = map.getView().getResolution();
    map.on('click', function (evt) {
        var url = wms_layer.getSource().getGetFeatureInfoUrl(
            evt.coordinate, viewResolution, viewProjection,
            {
                'INFO_FORMAT': 'text/html',
            });

        if (url) {
            document.getElementById('info').innerHTML =
            '';
            $('#myModal').modal('show');
        }
    });

Whenever I click on the map I always getting info of the HissarFullData:block_boundary layer. I am using Layer switcher with this .

Needing ArcGIS Server to edit MS SQL Server database?

Here is what I have:

ArcGIS Desktop Advanced (10.2) & MS SQL Server 2008 R2

Here is my problem:

I need to be able to connect to the ms sql server and import, export, create, and edit spatial data. I can connect to my database and export features to it using the "Feature Class to Geodatabase" tool but can't edit the data once bringing it back into an ArcMap session. I have tried using the "Create Enterprise Geodatabase" tool available with the advanced desktop license but it asks for an Authorization File (for ArcGIS for Server I assume) that I don't have. I have looked at prices for ArcGIS for Server and it is not feasible option.

Are there alternatives to achieve what I need to using what I currently have licenses for?

Where does arcSDE fit into this?

Juno 3D open source software

Is there an open source software that I can use to collect data with the Juno 3D GPS?

I'm hesitant to buy a Terrasync or ArcPad license.

qgis 3 - QgsPrintLayout Setup from PyQGIS 3?

I am interested in making reproducible PDF figures in QGIS3, and therefore want to work from the console/PyQGIS.

Using some previous helpful advice from PDF Output with Python the problem is successfully reduced to setting up the QgsPrintLayout. I am not sure how to do this, below is my best attempt so far:

#setup the background classes for managing the project
canvas =iface.mapCanvas()
manager = project.layoutManager()
#new layout
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName('console')

itemMap = QgsLayoutItemMap(layout)
itemMap.setExtent(canvas.extent())
layout.addLayoutItem(itemMap)
#add the layout to the manager
manager.addLayout(layout)

The "console" layout is successfully added to QGIS 3 and I can open it to see if the content is satisfactory. The map does not render and its extent is nan,nan,nan,nan.

How do I setup the extent and size of the map object in the layout?

---

SOLVED. The following code produces a map without issues.

#get a reference to the layout manager
manager = project.layoutManager()
#make a new print layout object
layout = QgsPrintLayout(project)
#needs to call this according to API documentaiton
layout.initializeDefaults()
#cosmetic
layout.setName('console')
#add layout to manager
manager.addLayout(layout)

#create a map item to add
itemMap = QgsLayoutItemMap.create(layout)
#using ndawson's answer below, do this before setting extent
itemMap.attemptResize(QgsLayoutSize(6,4, QgsUnitTypes.LayoutInches))
#set an extent
itemMap.setExtent(canvas.extent())
#add the map to the layout
layout.addLayoutItem(itemMap)

Answer

You need to set a size for the map item (usually you do this before setting the extent):

itemMap.attemptResize(QgsLayoutSize(50, 60, QgsUnitTypes.LayoutMillimeters)) itemMap.attemptMove(QgsLayoutPoint(1,2,QgsUnitTypes.LayoutMillimeters)) itemMap.setExtent(...)

gdal - Error when clipping raster in QGIS: "'ascii' codec can't encode character"?

I have done it many times (clipping a raster with a polygon layer or by extension) but today I got an error for the first time and I don't know how to solve it. I am always using the same type of raster layers. When running the clipping tool I got the following message:

gdal_translate -of GTiff -ot Byte -projwin 410282.185971 4694950.41162 410783.882475 4694524.96815 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co ZLEVEL=6 -co TFW=YES "C:/Users/FERRAN ALA/Desktop/Justificació Fontanals/QGIS/Topo 1_5000/bt5mv20sd0f283077st1r031.sid" "C:/Users/FERRAN ALA/Desktop/Justificació Fontanals/QGIS/background/image.tif"

GDAL command output:

'ascii' codec can't encode character u'\xf3' in position 201: ordinal not in range(128) See log for more details

Does anybody know what that is?

java - Proj4J not precise for EPSG:3857 transformations

I try to convert some coordinates from "EPSG:31469" to "EPSG:3857" and back.

With the heavy framework GeoTools it seems to work pretty precise.

"EPSG:31469" - 5439627.33,          5661628.09
"EPSG:3857"  - 1573655.6648492748,  6636624.730959651
"EPSG:31469" - 5439627.330626475,   5661628.099510074

But I would like to use the lightweight library Proj4j (https://github.com/Proj4J/proj4j) which produces massive deviation:

"EPSG:31469" - 5439627.33,          5661628.09
"EPSG:3857"  - 1573659.9666159092   6603334.624358153

"EPSG:31469" - 5439626.576916553    5661585.491150079

Any ideas why this is the case?

My code:

import org.cts.crs.CRSException;
import org.osgeo.proj4j.CRSFactory;
import org.osgeo.proj4j.CoordinateTransform;
import org.osgeo.proj4j.CoordinateTransformFactory;
import org.osgeo.proj4j.ProjCoordinate;


public class Test
{

    public static void main(String[] args) throws CRSException
    {
        ProjCoordinate coord = new ProjCoordinate(5439627.33, 5661628.09);
        System.out.println(coord);

        ProjCoordinate target = new ProjCoordinate();
        CRSFactory crsFactory = new CRSFactory();

        CoordinateTransformFactory f = new CoordinateTransformFactory();
        CoordinateTransform t;

        t = f.createTransform(crsFactory.createFromName("EPSG:31469"),
            crsFactory.createFromName("EPSG:3857"));

        t.transform(coord, target);

        System.out.println(target);


        ProjCoordinate coord2 = new ProjCoordinate(target.x, target.y);

        t = f.createTransform(crsFactory.createFromName("EPSG:3857"),
            crsFactory.createFromName("EPSG:31469"));

        t.transform(coord2, target);

        System.out.println(target);
        System.out.println(target.x - coord.x);
        System.out.println(target.y - coord.y);


    }

}

Answer

Thanks to AndreJ! The fork of proj4j at https://github.com/dwins/proj4j works precisely on EPSG:3857

arcgis desktop - Adding PDF to backside printing for Data Driven Pages?

Is there an easy way to add a PDF to the backside of each map printed using data driven pages?

I know I could print one sided and then reinsert the printed pages upside down, and print the PDF on the reverse side. However, we may deliver this to our client digitally and would like it to be setup for them.

For example, I have data driven pages set up for 50 parcels. Each map printed will show a single parcel's data. I need a generic PDF to print on the backside of each of these 50 maps.

Answer

To do this I would use ArcPy with Data Driven Pages. There is some sample code in the help document entitled Creating a map book with facing pages which has all the necessary coding patterns but is a bit more complicated than what you want to do.

To create your PDF you just need to use arcpy.mapping to loop through your data driven pages one at a time, and then on each iteration to:

1. export the current page using exportToPDF and then appendPages


2. use appendPages again to put a copy of your generic PDF next

postgresql - Closest line to a point using PostGIS?

I have a table of roads and a table of points. I'd like for each point to have the ID of the closest road, as well as the distance to that road.

The line data has both the path type and the geometry type for its coordinates. The points have XY coordinates. What would be the best way to go about finding the closest road using PostGIS?

animation - Animating objects on canvas to change color by delay using PyQGIS?

On my map canvas I have various point objects spread through several layers. My plan is to 'animate' those objects by a loop from start to a predefined resolution.

By now I'm just testing the algorithm with the vertical positions of the existing objects and an exportCanvas function which saves the canvas as a single image for each step.

Because this method isn't very handy I want to have the output of the loop directly to the canvas - delayed by e.g. 0.5 seconds per step.

Above is an example of the merged jpg-files to a gif

The problem is that my loop doesn't work very well at all and freezes my QGIS for the whole process time (+delay).
I know that I need some kind of hyperthreading or callback from QGIS to my plugin but at this point I have no clue how to transfer this idea to useful code

Down below: the code so far

# Animation: Button clicked
def btnAnimationClicked(self):       
    rate = 24

    # Werte zum Färben der xyz
    for i in range(0,rate,1):
        self.caller(i, self.callback)

# Animation: Caller
def caller(self, val, func):

    func(val)

# Animation: Callback
def callback(self, i):
    yhi = 5699049523 + 400000 
    ylo = 5696625163
    delta = yhi - ylo
    rate = 24
    ste = delta / rate 


    ubound = ylo + (ste * (i+1))
    lbound = ylo + (ste * i)

    values = (
        ('hi', ubound + 1, yhi, 'red'),
        ('is', lbound, ubound, 'green'),
        ('lo', ylo, lbound - 1, 'yellow'),
    )           
    # HaKasten abgestuft färben     
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():

        if layer.name()[0:9] == "xyz_______":    
            ranges = []
            for label, lower, upper, color in values:
                symbol = QgsSymbolV2.defaultSymbol(layer.geometryType())
                symbol.setColor(QColor(color))
                rng = QgsRendererRangeV2(lower, upper, symbol, label)
                ranges.append(rng)
            attribut = 'y'
            renderer = QgsGraduatedSymbolRendererV2(attribut, ranges)
            layer.setRendererV2(renderer)    

    for layer in qgis.utils.iface.mapCanvas().layers():
        layer.triggerRepaint()
    #time.sleep(0.5) # optional
    self.exportCanvas() 

open source gis - QGIS: Is it possible to edit the datasource of a PostGIS layer?

After creating a layer with RT_SQL_Layer plugin, I want to change some details in the query. Is this possible, or do I need to create a new layer every time I want to change something in the query?

More specifically; I want to change the "Source for This Layer" as shown in example bellow:

Answer

You cannot change the query used in RT SQL Layer plugin. But you can try to only use a general query in RT and refine it in Layer Properties - General - Query builder. The part in Query builder can be easily adjusted while working with the same layer.

transportation - How can I calculate the average speed in various legs of a bus route?

I'm working on analyzing NYC MTA Bus traffic data.

There are numerous buses traveling along the route, emitting gps data at random locations.

I have about 30k records with bus_line, vehicle_id, longitude, latitude, bus_speed and timestamp. I refer to this dataset as dataset1

Calculating the speeds for an individual bus is trivial but I can't think of a simple way of averaging the speed in a particular leg along the route.

One idea I'm thinking about is to create another dataset (i.e., "dataset2") with longitude & latitude in 50 ft intervals. Then find a way to assign its points to the points from dataset1 based on proximity. At this point, it would be possible to find the average speed for each point because there will be multiple records with the same longitude and latitude.

The biggest challenge I see is in assigning the points from dataset2 to the points in dataset1: you need to sort the points exactly in the order they "appear" in the route. You can't simply order by longitude or latitude. You can sort conditionally such that if the bus is going west, you sort by longitude. If it's going north, you sort by latitude. However, I'd like something simpler if it exists. I'd appreciate any help!

Thanks.

arcgis desktop - Opening MXD saved in one version using earlier version

Working in an enterprise system. We need to determine if everyone needs to update to 10.6 at the same time for mxd cross-compatibility or can we update from 10.5 as needed.

snapping - Is there any option in QGIS to draw parallel lines that snap on the outside part of the line (adjacent lines) that are independent from scale?

I have to draw a map of a public bus transport network and i have a problem in the avenues where a lot of lines passes by. My objective is to make a visual map where the lines in those avenues are parallel one from the other and with no space between them, in order to see all of them at the same time, regardless the scale.

I've tried snapping options, and parallel drawing (CAD tools), but the issue is that when i change the scale, the distance between lines change and they collide or separate one from the other, messing everything.

So my question is:

• Is there any option to draw adjacent lines that stays adjacent independently of the scale of visualisation?

Note: I don't know if adjacent is used for what I mean, so here I copy an example of what I'm searching for:

---

This is what happens when I zoom in:

And when I zoom out:

EDIT: Here I share a situation of why the suggested solution wouldn't work for me. Three lines share the same street for a while, but then they separate into three different streets.

I can't use a single line with several symbolisation because all along the network, the lines split and join again (there are more than 15 bus lines).

The data contained in the attribute data doesn't help me, because it's only a layer full of lines, with no attribute other than some network information (number of buses/hour, passenger/hour, etc.). (QGIS 2.8.6 working on Windows 7)

qgis - Viewer not displaying the data in my standalone PyQGIS Application

I tried creating a standalone external application of QGIS using python and I followed the tutorial mentioned in the link http://download.osgeo.org/qgis/doc/workshops/foss4g2007_qgis0.9_workshop_en.pdf

After running the python file the viewer opens and allows the user to browse for Shapefiles. The problem is that the Shapefile cannot be viewed in the viewer. It raises an error "Invalid Shapefile" though I added a correct Shapefile everytime I loaded any shapefile.

My code as well as the error is attached. Please tell me where I am wrong.

def addLayer(self):
    file = QFileDialog.getOpenFileName(self, 
               "Open Shapefile", ".", "Shapefiles (*.shp)")

    fileInfo = QFileInfo(file)

    # Add the layer
    layer = QgsVectorLayer(file, fileInfo.fileName(), "ogr")

    if not layer.isValid():
        raise IOError("Invalid Shapefile")

    # Change the color of the layer to gray
    symbols = layer.renderer().symbols()

    symbol = symbols[0]
    symbol.setFillColor(QColor.fromRgb(192, 192, 192))

    # Add layer to the registry
    QgsMapLayerRegistry.instance().addMapLayer(layer)

    # Set extent to the extent of our layer
    self.canvas.setExtent(layer.extent())

    # Set up the map canvas layer set

    c1 = QgsMapCanvasLayer(layer)
    layers = [c1]
    self.canvas.setLayerSet(layers)

And the error window is

Answer

Non-valid layers are most of the times due to a wrong QGIS prefix definition.

Please try with:

qgis_prefix="C:\\Program Files\\QGIS Wiena\\apps\\qgis"    
QgsApplication.setPrefixPath(qgis_prefix, True) 

Which should go right before this line:

QgsApplication.initQgis()

Now your tested layers should be valid.

coordinate system - Manipulating Azimuthal Equidistant Projections in QGIS

I have a shapefile of world countries that is projected as Azimuthal Equidistant with Chicago, USA as the center. I would like to change the projection such that Edinburgh, UK is the center. One would think that this should be as simple as changing 2 numbers in the coordinate system parameters...

The problem that I seem to be experiencing is that QGIS makes it difficult if not impossible to view the complete proj4 code of a given shapefile's coordinate system parameters.

This is what I have tried so far:

• Added shapefile of world countries projected as Azimuthal Equidistant with Chicago as the center


• Opened properties, metadata tab, copy proj4 parameter code


• Went to settings, custom projection, paste in proj4 code, changed the lat/lon from that of Chicago to that of Edinburgh


• Saved shapefile of world countries, specified my new custom projection, added new file to map



• My new file looks like a map of Pangaea, with the continents all smushed together.

I think there must be other parameters aside from those shown in the proj4 code of the metadata tab. The only thing I'm changing in the proj4 code is the center of the map, but the change I'm seeing is much more than that. My observation suggests that more than one independent variable is changing.

I would be happy to send my shapefile to anyone who would like to play with it.

Thanks,

Daniel Wolf Environmental/Geospatial Enthusiast

gdal - How do I find out the resolution in (centimeters per pixel) of a TIF?

I tried using GDAL and got:

# gdalinfo nicholas_map_2.tif 

Driver: GTiff/GeoTIFF
Files: nicholas_map_2.tif
       nicholas_map_2.tfw
Size is 28799, 27875
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],

    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-122.325695999999994,37.876846000000000)
Pixel Size = (0.000000441883000,-0.000000350276000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL

Corner Coordinates:
Upper Left  (-122.3256960,  37.8768460) (122d19'32.51"W, 37d52'36.65"N)
Lower Left  (-122.3256960,  37.8670821) (122d19'32.51"W, 37d52' 1.50"N)
Upper Right (-122.3129702,  37.8768460) (122d18'46.69"W, 37d52'36.65"N)
Lower Right (-122.3129702,  37.8670821) (122d18'46.69"W, 37d52' 1.50"N)
Center      (-122.3193331,  37.8719640) (122d19' 9.60"W, 37d52'19.07"N)
Band 1 Block=28799x32 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=28799x32 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 

Band 3 Block=28799x32 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=28799x32 Type=Byte, ColorInterp=Alpha

I assume pixel size is what I'm looking for? But is that in degrees?

Querying related tables in QGIS

I have a parcel layer and a table of property sales. I have the two tables related (it is a one-to-many relationship) in a QGIS project so that when you click on a parcel you see all the related sales information. Is there a way to query that related sales table and have it select the appropriate parcels directly in QGIS?

I have come up with a round-about way where I create a view in the Spatialite DB that the data is stored in but it is asking my end user a bit too much to be able to change the views parameters every time he needs a new sales map.

Answer

I have the same issue. If find an alternative solution based on a simple selection... (Duplicate issue i've made)

From the parent table, I retrieve a list of ID's (PK) and in the child table, I set a filter based on my list of ID's to my FK.

Important notice, I do not use the relationship, only a selection based on PK FK. There is an important performance issue on large tables.

Here a snippet.

li_id = []
layers = QgsMapLayerRegistry.instance().mapLayers()

for name, layer in layers.iteritems():
     if name.lower().startswith('NAME OF MY PARENT TABLE'):

        features = layer.selectedFeatures()
        for f in features:
            li_id.append(str(f.attribute('F_ID')))
cLayer = iface.mapCanvas().currentLayer() # (CHILD TABLE SELECTED IN TOC)
if len(li_id) ==1:
    cWhere = """"F_ID" = '{0}'""".format(li_id[0])
else:
    cWhere=  """"F_ID" in {0}""".format(tuple(li_id))



expr  = QgsExpression(cWhere)
it = cLayer.getFeatures( QgsFeatureRequest( expr ) )
ids = [i.id() for i in it]
cLayer.setSelectedFeatures( ids )
iface.showAttributeTable(iface.activeLayer())

postgresql - Load PostGIS layer/table to QGIS canvas using Python

I want to load a PostGIS layer to QGIS canvas. I am using the following code but it is not loading the layer in canvas.

sql = "(SELECT * FROM table WHERE id = 20)"
uri = QgsDataSourceURI()
uri.setConnection("host", "5432", "databasename", "user", "password")

uri.setDataSource("public","table", "the_geom", sql,"gid")
print uri.uri()
vlayer = QgsVectorLayer(uri.uri(), "test", "postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)
canvas.refresh()

But the layer is not loaded. I don't know what I am doing wrong here.

kml - Calculated distance doesn't match Google Earth

I have a google earth route that measures 190.1km in Google Earth.

I then pulled down the .kml file, extracted the points and measured them using Vincenty's formula (ellipsoid earth) and got 190.2.

BUT, if I measure the SAME points along the SAME route to a different distance along that route, I'm way off. Sometimes as much as 2km.

For example, If I want to know the lat and lon of the point closest to 100km, I measure the distance between all the points along the route until the measurement is as close to 100km as it's going to get. When I compare that point to the route in Google Earth, it's around 98km.

Any ideas of why this is the case?

gdal - Import error:No module named _gdal

I have installed python 2.5.1 and FWTools2.4.7 in my laptop(Windows 7) when i run my pythonscript in FWTools i get following error Import error:No module named _gdal can any one suggest how to solve this problem?

openlayers 2 - WMS GetFeatureInfo request returns nothing with GeoServer

I'm having trouble with a select feature example in GeoServer and OpenLayers. For some reason the request returns nothing.

Here's the stripped down version that produces this:

I have a small WMS layer(~100 parcels) with polygons, and address data that I am overlaying Google Maps with a projection. It displays correctly.

map = map(...); //standard map creation
gmap = ...; // standard street map creation
wmslayer = new OpenLayers.Layer.WMS(
            "Parcels",
            "http://192.168.0.205/geoserver/wms", 
            {'layers': 'test:parcel', 'format':'image/png', 'transparent':'true'},
            {'opacity': 0.5, 'isBaseLayer': false, 'visibility' : false, 'zoomMin' : 15}
);

Then I create the click control...

infoControls = 
           {
            click: new OpenLayers.Control.WMSGetFeatureInfo
            ({
                url: 'http://192.168.0.205/geoserver/wms', 
                title: 'Get parcel data',
                layers: [wmslayer],
                queryVisible: true

            })
        };

Added the layers, add the control and activate it

map.addLayers([gmap, wmslayer]);  
infoControls[0].events.register("getfeatureinfo", this, showInfo);
map.addControl(infoControls[0]);
infoControls.click.activate();

//show info callback

function showInfo(evt) {    
    console.log(evt);
    alert(evt.text);
}

But when I click on the map inside I parcel, it gives the alert is blank, and so is the features in Firebug (from console.log(evt)).

The GET request in Firebug is red.

BBOX    -10845148.716676,5165071.048543,-10843335.725912,5166002.624825
FEATURE_COUNT   10
FORMAT  image/png

HEIGHT  390
INFO_FORMAT text/html
LAYERS  test:parcel
QUERY_LAYERS    test:parcel
REQUEST GetFeatureInfo
SERVICE WMS
SRS EPSG:900913
STYLES  
VERSION 1.1.1
WIDTH   759

X   582
Y   246

Edit: The database is CRS is EPSG:4327 and projecting it to Google Mercator... Do I need to project the BBOX of the request back to EPSG:4327 before I send the WMS request? How can I do this?

Here is the GeoServer log of the request

2011-07-22 14:49:01,301 INFO [geoserver.wms] - 
Request: getFeatureInfo
    GetMapRequest = 
GetMap Request
 version: 1.1.1

 output format: image/png
 width height: 759,390
 bbox: ReferencedEnvelope[-1.0845637448332E7 : -1.0842011466804E7, 5164278.964868 : 5166142.117432]
 layers: test:parcel
 styles: polygon
    QueryLayers = [org.geoserver.wms.MapLayerInfo@a837896f]
    XPixel = 385
    YPixel = 144
    FeatureCount = 10
    InfoFormat = text/html

    Exceptions = application/vnd.ogc.se_xml
    Version = 1.1.1
    Request = GetFeatureInfo
    BaseUrl = http://192.168.0.205:8080/geoserver/
    Get = false
    RawKvp = {INFO_FORMAT=text/html, BBOX=-10845637.448332,5164278.964868,-10842011.466804,5166142.117432, QUERY_LAYERS=test:parcel, SERVICE=WMS, HEIGHT=390, REQUEST=GetFeatureInfo, STYLES=, WIDTH=759, FEATURE_COUNT=10, VERSION=1.1.1, FORMAT=image/png, LAYERS=test:parcel, Y=144, X=385, SRS=EPSG:900913}
    RequestCharset = null

Answer

So mostly for future readers of this question - when you are doing getFeatureInfo requests to a different server (which includes a difference in the port number) you need a proxy - see http://trac.osgeo.org/openlayers/wiki/FrequentlyAskedQuestions#ProxyHost for more discussion.

Burning stream network into DEM layer using ArcGIS Desktop?

I am working on a DEM that in an area where there is little relief.

I would like to 'burn' the river network into the DEM so I can calculate flow accumulation and flow length accurately.

I am using ArcGIS Desktop 10.

Replacing ArcPy with ArcGIS API for Python?

Today I found out about the ArcGIS API for Python. I was wondering if you could use this to (eventually) replace ArcPy?

Recently I have been doing analyses that require cloud based solutions. When using AWS, Docker and Jupyter it is very convenient to stick with a Unix based OS such as Ubuntu. However, the Unix alternatives for ArcGIS have their own limitations as well such as the always hard to install python bindings for the GDAL library. ArcGIS for Python operates cross platform.

Perhaps there are the things I should and should not expect from the ArcGIS API for Python.

Perhaps Arcpy is going to disappear altogether.

API in action

gdal - DXF to georeferenced SHP

I know how to set projection and georeference bounds of scanned TIFF file, by using gdal, but I don't know how to do the same for vector file.

For example I georeferenced TIFF file by using this command:

gdal_translate -of GTiff -a_ullr x1 y1 x2 y2 -a_srs "7041.prj" input.tif output.tif

where 7041.prj is Gauss Kruger projection file.

I used CorelDraw to do the tracing, and both raster and DXF output share same bounds - vector file has boundbox same as raster image used for tracing. I wanted they overlay automatically when loaded in GIS.

For vector file, with ogr2ogr I can convert DXF to SHP and set projection, but I can't figure how to set georeference bounds as ogr2ogr doesn't accept -a_ullr argument, probably for some reason.

So my question is can I somehow set bounds for resulting SHP file?

Answer

Now it's possible by using latest gdal revision, released couple of minutes ago. ogr2ogr has new switch -gcp:

Every once in a while, we have the need to take non-georeferenced data (e.g. DXF or other CAD files) and draw it on top of a map.

One way to handle this would be to add the ability to specify GCPs in ogr2ogr that would be used to define a transformation that would then be applied to all vector coordinates during the translation. (The gdaltransform program kind of does this for small sets of coordinates, not for whole files)

more info and examples: http://trac.osgeo.org/gdal/ticket/4604

qgis - How to generate multiple maps where one highlighted attribute changes?

I have a layer that contains counties within the state of Washington.

I need to build one map (jpg) for each county. For that specific map all counties will show as light grey, but the specific county will be green, and the file name will be the county name.

Answer

What you want is creating the Atlas.

1. In print composer, go to Atas generation tab and check Generate atlas.


2. As covarage layer set your counties layer, and as output set your name of counties attribute (somtheing like "counties_names").


3. Go to layer properties and set rule-based symbology, create empty rule with grey color and in this rule nest a rule like $id = @atlas_featureid with green color (see picture below).

4. Go back to print composer and start preview of atlas (in atlas panel or menu Atlas --> Preview Atlas. Now you can browse atlas features by Next feature and Previous feture and so on.


5. If you wnat add other items like text fields, scalebar etc. (You can also control text fields by expression, so you can get for each atlas feature uniqe text fields values, for example counties names.)


6. Export via atlas panel or menu Atlas --> Export Atlas as Images...

Possible output (one of exported images):

pyqgis - Writing a python processing script with QGIS 3.0

Following the update to QGIS 3.0, it has become very difficult to find any information concerning the writing of processing scripts in QGIS 3.0.

@Underdark (see here ) has provided a basis for the skeleton. This code also seems to have been added in QGIS, when writing a new script from template (QGIS 3.0.2).

However, I couldn't find any way to help Python newbies like me to understand how to change that code, especially for the input and output layers.

My goal is to write a script taking 2 raster layers and a double as input, outputting two layers.

What would be the changes required to the example code to allow that?

For QGIS 2.x I would have used the following syntax :

##Layer1=raster
##Layer2=raster 
##myDouble=Double
##OutLayer1=output raster

##OutLayer2=output raster   

From what I understand, the changes have to be made in the following procedure, but I'm not sure what to put in place.

def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterFeatureSource(
        self.INPUT,
        self.tr("Input layer"),
        [QgsProcessing.TypeVectorAnyGeometry]))
    self.addParameter(QgsProcessingParameterFeatureSink(
        self.OUTPUT,

        self.tr("Output layer"),
        QgsProcessing.TypeVectorAnyGeometry))

On May 16th, the QGIS Python API documentation was released. However it is still unclear for me how to use it here. (Which might very well be a lack of Python knowledge)

Answer

With the transition from QGIS2.x to QGIS3.x the whole processing framework has been reworked and large parts of it run now as C++ classes that you can interact with using Python. Unfortunately the simple parameter syntax for data/dataset IO are no longer valid. The new parameter structure is much more orientated after the builtin (Python-) Processing algorithms that you find preinstalled in the toolbox.

As I see, you already followed the description of the new algorithm structure by @underdark. But in order to adjust this structure for your requirements (raster layers, double input, etc.) you have to change the code at multiple locations in the script. I have coded a rough example with a short explanation for you (just an algorithm skeleton based on @underdarks example):

from qgis.PyQt.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsProcessing, QgsProcessingAlgorithm, 
QgsProcessingParameterRasterLayer,QgsProcessingParameterNumber, 

QgsProcessingParameterRasterDestination)

class RasterAlg(QgsProcessingAlgorithm):
    INPUT_RASTER_A = 'INPUT_RASTER_A'
    INPUT_RASTER_B = 'INPUT_RASTER_B'
    INPUT_DOUBLE = 'INPUT_DOUBLE'
    OUTPUT_RASTER_A = 'OUTPUT_RASTER_A'
    OUTPUT_RASTER_B = 'OUTPUT_RASTER_B'

    def __init__(self):

        super().__init__()

    def name(self):
        return "RasterAlg"

    def tr(self, text):
        return QCoreApplication.translate("RasterAlg", text)

    def displayName(self):
        return self.tr("RasterAlg script")


    def group(self):
        return self.tr("RasterAlgs")

    def groupId(self):
        return "RasterAlgs"

    def shortHelpString(self):
        return self.tr("RasterAlg script without logic")


    def helpUrl(self):
        return "https://qgis.org"

    def createInstance(self):
        return type(self)()

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterRasterLayer(
            self.INPUT_RASTER_A,
            self.tr("Input Raster A"), None, False))

        self.addParameter(QgsProcessingParameterRasterLayer(
            self.INPUT_RASTER_B,
            self.tr("Input Raster B"), None, False))
        self.addParameter(QgsProcessingParameterNumber(
            self.INPUT_DOUBLE, 
            self.tr("Input Double"), 
            QgsProcessingParameterNumber.Double,
            QVariant(1.0)))
        self.addParameter(QgsProcessingParameterRasterDestination(
            self.OUTPUT_RASTER_A,

            self.tr("Output Raster A"),
            None, False))
        self.addParameter(QgsProcessingParameterRasterDestination(
            self.OUTPUT_RASTER_B,
            self.tr("Output Raster B"),
            None, False))

    def processAlgorithm(self, parameters, context, feedback):
        raster_a = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER_A, context)
        raster_b = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER_B, context)

        double_val = self.parameterAsDouble(parameters, self.INPUT_DOUBLE,context)
        output_path_raster_a = self.parameterAsOutputLayer(parameters, self.OUTPUT_RASTER_A, context)
        output_path_raster_b = self.parameterAsOutputLayer(parameters, self.OUTPUT_RASTER_B, context)

        #DO SOME CALCULATION

        results = {}
        results[self.OUTPUT_RASTER_A] = output_path_raster_a
        results[self.OUTPUT_RASTER_B] = output_path_raster_b
        return results

Which steps are done?

1. Import all necessary classes.


2. Define the algorithm as a class inheriting from QgsProcessingAlgorithm.


3. First you have to declare the names of the input and output parameters as string variables (parameter names) of the algorithm class (ie. INPUT_RASTER_A = 'INPUT_RASTER_A') in order to reference your algorithm with the parameters provided by the processing framework.


4. Add the methods that wire your algorithm to the processing toolbox gui and provide helpstrings, etc.


5. Then you add the parameters of the processing framework. Those are defined as child classes of QgsProcessingParameterType - in the case of your algorithm: QgsProcessingParameterRasterLayer, QgsProcessingParameterNumber, and so on. You can consult the API entries (ie. QgsProcessingParameterRasterLayer) in order to pass the right arguments and construct the parameter objects.


6. Pass the parameters alongside context and feedback objects to the processAlgorithm() method where you obtain the input datasets from the parameters at runtime (in this case QgsRasterLayer objects by using the parameterAsRasterLayer() method, etc.).


7. Do your computation.



8. Add the outputs to the results dictionary and return them as result of calling processAlgorithm().

I hope I could give you some insights on how to design your python algorithms in QGIS3. Whenever you are stuck, it is always helpful to look at how existing algorithms of the processing framework handle the parameters. You can have a look at them here.

geoprocessing - arcpy.MakeFeatureLayer in-memory layer still exists when subsequent step fails during testing

When I used arcgisscripting to create in-memory output layer with gp.makefeaturelayer, I would put it in a try/except block and delete the gp object in the except piece when the makefeaturelayer would fail. Now when I am using arcpy, it looks like I need to close the idle window and open again in order to get the in-memory output layer to get deleted.

How could I have the output layer from my makefeaturelayer tool get deleted in a try/except block? Thanks.

Answer

arcpy.Delete_management(featureLayer)

Is there QGIS plugin to allow 3d visualisation of geological borehole data similar to functionality of Target for ArcGIS?

I'm very new to QGIS, I have started using it as my new company is not willing to pay for ArcGIS and Target at this stage so I'm hoping QGIS can substitute not just in the interim but for good.

What I am looking for is a way to load in my borehole database and create 3d visualizations of the data, make cross sections and if possible model seams. Is there any capability within the QGIS world to undertake this kind of work?

javascript - Leaflet mouse wheel zoom only after click on map

I'm working with the Leaflet JavaScript Library and attached a (working) map to my HTML Document. It is in the middle of the page, and when I'm scrolling down with my mouse wheel and arrive at the map, it automatically zooms into the map.

I want to scroll through the page without stopping at the map. Is there a way to activate the wheel zoom only after the first click on the map?

Answer

It's simple: create L.Map with scrollWheelZoom: false option, then add a listener:

map.once('focus', function() { map.scrollWheelZoom.enable(); });

If you need to toggle zooming:

map.on('click', function() {

  if (map.scrollWheelZoom.enabled()) {
    map.scrollWheelZoom.disable();
    }
    else {
    map.scrollWheelZoom.enable();
    }
  });

arcgis desktop - Adding automated values to attribute field?

I have parcel and buildings layers.

In their attribute tables I want to replace "NULL" values with numbers from 1 to 2000, in "Parcel ID", or "Building ID" field.

How can I do this in ArcGIS Desktop?

Answer

Here is a code block for the Field Calculator that will do what you require.

1. Open the attribute table


2. Select all of the records that contain a NULL value in the desired field


3. Open up the Field Calculator and insert this code in the appropriate sections.

Accumulative and sequential calculations
Calculate a sequential ID or number based on an interval.

Expression:
autoIncrement()

Expression Type:
PYTHON_9.3

Code Block:

rec=0
def autoIncrement():
 global rec
 pStart = 1 #adjust start value, if req'd 
 pInterval = 1 #adjust interval value, if req'd
 if (rec == 0): 
  rec = pStart 
 else: 
  rec = rec + pInterval 
 return rec

This code sequence is taken from the ArcGIS 10.0 Desktop Help: Calculate Field Examples

I would highly recommend that you browse through the help documentation in that file and in the other help sections for ArcGIS, as there are many useful examples and code samples to pull from.

Enter the code above, into the Field Calculator as shown in the screen shot below: You would enter the code into the sections as detailed below:

1. Choose Python as the Parser


2. Pre-Logic Script Code: Enter code listed under Code Block


3. Under the field name, in this case, RD20FULL =: Enter code listed under Expression


4. Press OK

arcgis javascript api - Convert spatial Geographic data from sql server 2008 to Json

Is there any example to convert spatial data (point,line,polygon) stored as geographic in sql server 2008 into Json and display data using Arcgis Javascript Api 3.0.I have tried converting geographic data into Json given Geometry to Json but as my data stored in geographic data format it failed.Thank you in advance.

Answer

In theory just replace all SQL commands having geometry:: cast to geography casts (assuming that you 2008 R2 MSSQL server) Edit: looked that example code . Change Microsoft.SqlServer.Types.SqlGeometry to Microsoft.SqlServer.Types.SqlGeography OR Find SQL command which select code and force it to cast your geography objects to geometry

css - Using OpenLayers/jQuery without CDN

I was wondering what are the necessary files for hosting a webgis server with OpenLayers3 without using a CDN?

I am testing on a localhost without internet access. It doesn't work to merely copy the direct js (and css) file pointed to by a CDN such as

http://openlayers.org/en/v3.0.0/build/ol.js

to the local webserver as .

Also, the same question holds for jquery and bootstrap. I had

But I wasn't sure how to get a complete local version for those.

Answer

Download the v3.0.0 distribution from https://github.com/openlayers/ol3/releases , then copy the js file(s) inside the build/ subdirectory to your webserver.

convert - Shapefile to MSSQL with ogr2ogr fails to make connection

I am using ogr2ogr from OSGeo4W on my Windows 7 box to import a shapefile into a SQL Server database on one of my servers. I was able to successfully complete the conversion using shp2sqlserver, but I can't figure out why it fails when I try the same thing with ogr2ogr directly.

Here's the syntax of my ogr2ogr command:

ogr2ogr -f MSSQLSpatial "MSSQL:server=[server];database=[dbname];User Id=[uid];Password=[pwd]" C:\Path\to\file\shapefile.shp

When I run this, it fails first with an error saying "Cannot open database [dbname] requested by the login. The login failed." Next is another error that the "MSSQL Spatial driver doesn't currently support database creation. Please create database with the Microsoft SQL Server Client Tools."

I have checked and double-checked that the user id and password given works properly, that there is already a database with the name given (I substituted [dbname] here to speak in generalities), and that the database is owned by the user.

As mentioned, I was able to do this with shp2sqlserver so I don't think it's my database setup that is the problem. Is there a problem with my connection string? Is there some peculiarity within ogr2ogr that I'm not aware of?

I am using a Windows 7 machine connecting to SQL Server 2008 R2. Gdal is version 1.9.1.

Answer

The ODBC documentation suggests keywords UID and PWD instead of User ID and Password.

ArcGIS 10.3: How to merge lines divided by junction points?

I have one shapefile with train stations (points) from a region in the UK and another shapefile with railroad tracks (lines) between these stations. Problem is, in the latter I don't have a unique line from one station to the next, but multiple connected segments. The only thing that I was able to do was to add another shapefile with junction points, to just visualize where the segment endpoints are.

How can I merge these multiple lines together, so that from station A to station B there is only one line? This is intended to be repeated for any pair of connected stations.

I would also be keen on using Python to do this, making sure there is a way for me to tell the script to merge all the lines between two stations, ignoring the junction points.

I have already tried both Unsplit lines and Dissolve, but the resulting lines are longer than they should be, e.g. they incorporate one too many segments, thus going past their intended endpoint. Unsplit lines seems to be the closest match, but I'd need to be able to restrict the merging so that it is only performed when junction points stand in the way (see image below), thus stopping at station points. If someone came up with a manual procedure, that'd be ok as well, as my dataset is not massive and I could handle it.

For purposes of clarity, see the image below. Currently, I have a number of segments and junction points in between two consecutive stations. The desired output is, instead, a unique object (line) connecting the two stations.

EDIT

Dissolve does not work as the segments composing the A-to-B route do not share any common field value.

Answer

This would be a two or three step process. First dissolve all lines into a single feature (do not indicate a dissolve field). Next, split lines at points using your station feature class as the point input. You may then need to perform a multipart to singlepart as a final step.